A Pair of Warm Giant Planets near the 2:1 Mean Motion Resonance around the K-dwarf Star TOI-2202^*

2.1.1. TESS

TOI-2202 was observed in five out of the 13 Sectors of the first year of the TESS primary mission. Observations were performed with the 30 minute cadence mode in Sectors 1, 2, 6, 9, and 13 between 2018 and 2019 July. TOI-2202 b was identified in the light curves extracted from the TESS Full Frame Images (FFIs). This was done using a pipeline called tesseract ³² (F. Rojas et al. 2021, in preparation). Pipeline tesseract receives the TIC ID 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. In the context of the WINE collaboration, we generated light curves for all Sectors of the first year of the TESS primary mission for stars in the TICv8 catalog brighter than T = 13 mag. Transiting candidates are identified by running the transitleastsquares package (TLS; Hippke & Heller 2019) on each light curve, and also by searching for individual negative deviations from the median flux. This latter processing allowed us to identify long-period planets (e.g., Schlecker et al. 2020a) and single transiters (Gill et al. 2020). TOI-2202 was identified as a candidate in individual TESS light curves using both detection methods, finding a periodicity of the transits of ∼12 days. Figure 1 shows the target pixel file (TPF) image of TOI-2202 constructed from the TESS image frames and Gaia DR2 data. We investigated if the transit signal was coming from neighboring stars by generating light curves for each of the pixels in a region 20 pixels wide around the target star, and also by analyzing the time series associated with the background flux. We found that the signal indeed originated close ($\lt 2^{\prime}$) to the target star. Nonetheless, as is shown in Figure 1, TOI-2202 has a slightly fainter (T = 12.7) neighboring star located at ∼24'' from it (TIC358107518) that was not possible to fully reject as the source of the observed transits. Our ground-based photometric and Doppler follow-up confirmed that the source of the transit signal was indeed TOI-2202 (see Sections 2.1.2 and 2.2).

Zoom In Zoom Out Reset image size

Since tesseract does not correct for contamination in the TESS apertures from nearby stars, we also apply a dilution correction for the contamination of TIC 358107518. To estimate the dilution, we use the R−P Gaia DR3 fluxes of TOI-2202 and TIC 358107518, and Equation (6) in Espinoza et al. (2019). For TOI-2202 the mean R−P flux is 113,856 electron s⁻¹, while the mean R−P flux of TIC 358107518 is 66,547 electron s⁻¹, and we obtain a dilution factor of 0.63, which was applied to the light curves.

The top panel of Figure 2 shows the dilution-corrected, median-normalized, TESS FFI light curves of TOI-2202. All of the TESS 30 minute cadence FFI data taken in Sectors 1, 2, 6, 9, and 13 suffered from notable systematics, which could be attributed to stellar activity, instrumental effects, or a combination of both. The transit events with a depth of approximately 1%, which was consistent with being produced by a Jupiter-sized planet with a period of about 11.9 days, can be also easily identified on the combined TESS light curve. The bottom panel of Figure 2 shows a GLS periodogram of the TESS light curve with prominent peaks in the range between 20 and 30 days, which come from the light-curve systematics, and which are generally in line with the most likely rotational period range of TOI-2202 (see Section 3).

Zoom In Zoom Out Reset image size

2.1.2. CHAT

2.2.1. FEROS

2.2.2. HARPS

2.2.3. PFS

For data and orbital analysis, we employed the Exo-Striker exoplanet toolbox ³⁴ (Trifonov 2019). The Exo-Striker provides easy access to a large variety of public tools for exoplanet data analysis, such as a generalized Lomb–Scargle periodogram (GLS; Zechmeister & Kürster 2009), a maximum likelihood periodogram (MLP; Baluev 2008, 2009; Zechmeister et al. 2019), transit photometry de-trending via wotan (Hippke et al. 2019), and a transit period search via the transitleastsquares package (TLS; Hippke & Heller 2019), which we used in this work for the RV and transit photometry signal analysis. The Exo-Striker is able to model data with a standard Keplerian model or with a more complex dynamical model in the case of gravitationally interacting multiple-planet systems detected in RVs or transit photometry data. The Exo-Striker works in the Jacobi coordinate system, which is a natural frame for the orbital parameterization of multiple-planet systems (Lee & Peale 2003). The modeling can be performed by "best-fit" optimization schemes (i.e., Levenberg–Maquardt, Nelder–Mead, Newton, etc.), sampling schemes such as an affine-invariant ensemble Markov Chain Monte Carlo (MCMC) sampler (Goodman & Weare 2010) via the emcee package (Foreman-Mackey et al. 2013), or the nested sampling technique (Skilling 2004) via the dynesty sampler (Speagle 2020).

The RV modeling schemes are intrinsic ³⁵ to the Exo-Striker, while currently, for transit light-curve models, the Exo-Striker employs the BAsic Transit Model cAlculatioN package (batman; Kreidberg 2015). In addition, dynamical modeling of TTV data is done with a Python wrapper of the TTV-fast package (Deck et al. 2014). When needed, RV and transit data can be additionally modeled with Gaussian process (GP) regression models using the celerite package (Foreman-Mackey et al. 2017a), which is also included in the Exo-Striker.

We inspected the TESS FFI light curves of TOI-2202 derived with tesseract. For our preliminary transit light-curve analysis, we further de-trended each of the tesseract sector light curves with a robust (iterative) Matérn GP kernel, with the aim to capture the systematic variation of the light curves (see Hippke et al. 2019). For CHAT data, we applied additional transit photometry de-trending as a function of the airmass at the epoch of observations. Our de-trending scheme resulted in nearly flat, normalized, TESS and CHAT light curves, which we inspected for transits using the TLS algorithm. Figure 3 shows the constructed TLS power spectra of the available transit data of TOI-2202. A significant peak occurs at a period of P_b = 11.91 days, followed by peaks at 5.96 days and 23.83 days, etc., which are simply low-order harmonics of the actual transit signal.

Zoom In Zoom Out Reset image size

However, a single-planet transit model with a period and phase adopted from the TLS failed to produce a good fit to the available TESS and CHAT data. Performing a transit light-curve model optimization by adjusting the orbital period P, eccentricity e and argument of periastron ω (or $e\sin \omega$, $e\cos \omega$), inclination i, time of inferior transit conjunction t₀, the semimajor axis relative to the stellar radius a/R_⋆, and planetary radius relative to the stellar radius R/R_⋆ did not help either. We found that the light-curve transit signals exhibit strong deviations in the expected individual time-of-transits assuming a constant period, suggesting strong transit timing variations (TTVs). To extract the TTVs, we performed a separate one-planet fit to the TESS and CHAT light curves, assuming a circular orbit (i.e., $e\sin \omega$, $e\cos {\omega }_{b}$ = 0), but with variable mid-transit times as fitting parameters. The rest of the transit parameters in this modeling scheme across each individual TESS and CHAT light curve were shared, with exception of the limb-darkening (LD) coefficients of TESS and CHAT, for which we adopted separate quadratic LD models, and nuisance parameters such as the transit light-curve relative photometric offset and additional data jitter. ³⁶ To perform an adequate parameter search, we ran a nested sampling scheme with 1000 "live-points," focused on the posterior convergence instead of Bayesian evidence (see, Speagle 2020, for details). We adopted the 68.3% confidence levels of the nested sampling posterior probability parameter distribution as 1σ parameter uncertainties.

The extracted mid-transit time estimates and their precise uncertainties yielded very strong periodic TTVs in the TESS and CHAT data, whose amplitude around the mean period was ∼1.2 hr, covering one full TTV superperiod. Table 2 lists the precisely extracted individual transit times and their errors. No significant TLS power is detected in the TESS photometry residuals, meaning that only one planetary companion of TOI-2202 is detectable on the TESS FFI light curves. These results suggest the presence of an additional nontransiting companion in the TOI-2202 systems that is close and massive enough to perturb the Jovian-like transiting planet TOI-2202 b.

For a period search in the precise RVs and activity indices data of TOI-2202, we computed maximum likelihood periodograms, which calculate the log-likelihood ($\mathrm{ln}{ \mathcal L }$) power by optimizing for each test frequency. The MLP algorithm allows for multiple data sets, each with an additive offset and jitter parameters (Baluev 2009; Zechmeister et al. 2019), which makes it more suitable for period analysis of multi-instrument data. We adopted the significance thresholds of the ${\rm{\Delta }}\mathrm{ln}{ \mathcal L }$ improvements, which correspond to false-alarm probabilities (FAPs) of 10%, 1%, and 0.1%.

Figure 4 shows the MLP periodograms of the combined FEROS and HARPS RVs and activity time series, separately. The PFS data consist of only four RVs and these were not included in the MLP analysis. The MLP periodogram of the combined FEROS and HARPS data shows a strong power at 11.91 days, which is the period of the known transiting planet candidate. The RV residuals, after removing the 11.9 days signal, however, do not indicate the presence of additional significant periodicity in the RV data. No significant activity periodicity is evident in the FEROS and HARPS activity indices. The only exception is the HARPS Na i D activity index data, which show marginally significant ${\rm{\Delta }}\mathrm{ln}{ \mathcal L }$ at lower frequencies that, however, do not have a counterpart in the HARPS RVs. Thus, we concluded that the MLP analysis of the activity data alone does not suggest that TOI-2202 is an active star.

Zoom In Zoom Out Reset image size

Further RV analysis in this work was performed jointly with just the available TTVs (Section 4.4) and transit photometry data (Section 4.5) of TOI-2202, using dynamical modeling.

Logically, our next step was to study the TTVs obtained from the TESS and CHAT light-curve transit events. TTVs contain essential information of the system's dynamics (Agol et al. 2005), and adequate modeling could reveal the most likely orbital parameters and planetary masses. This is under the condition that at least one complete TTV superperiod (i.e., one cycle of the TTV signal) is covered (Lithwick et al. 2012). We detected just over one full TTV superperiod of the transiting TOI-2202 b, with an estimated period of approximately 357 days. Although this is a promising lead, following Lithwick et al. (2012), it is evident that the observed superperiod and the TTVs amplitude can be explained with almost all possible period ratios in the first- or second-order commensurability, in which TOI-2202 c could be inner or outer planet. Therefore, we examined wide ranges of the period, mass, and eccentricity of the nontransiting TOI-2202 c.

First, we performed TTV fitting with the Exo-Striker by adopting a self-consistent dynamical model on the extracted TTVs. The fitted parameters for each planet were the dynamical planetary mass m, orbital period P, eccentricity e, argument of periastron ω, and mean-anomaly M₀, which in this work are always valid for the epoch of the first transit event t₀ = 2458327.103 (BJD). For this test we assumed a coplanar, edge-on, and prograde two-planet system (i.e., i_b , i_c = 90° and Δi = 0°), whereas for the dynamical mass of TOI-2202, we adopted our best estimate of 0.823 M_⊙. The time step in the dynamical model was set to 0.1 days to assure a precise orbital resolution of the inner transiting planet TOI-2202 b.

We ran a nested sampling scheme, which allowed us to efficiently explore the complex parameter space of osculating orbital parameters and study the parameter posteriors. In our nested sampling scheme with dynesty, we ran 100 "live-points" per fitted parameter, focused on the Bayesian log-evidence convergence, using a "static" nested sampler (see Speagle 2020 for details). For all parameters, we adopted uniform priors, which define an equal probability of occurrence within the experimentally chosen parameter ranges. For TOI-2202 b, we selected prior parameter range estimates taken from our single-planet parameter analysis from Section 4.2, which assure a transit occurrence near t₀ = 2458327.103, whereas for the perturber in the system (to become TOI-2202 c), we explored a wide parameter space of eccentricities, masses, and periods. For instance, we explored ${P}_{c}\in { \mathcal U }$(20.0, 40.0) days, ${e}_{c}\in { \mathcal U }$(0.0, 0.2), and ${m}_{c}\in { \mathcal U }$(0.01, 0.6) M_Jup; the rest of the priors shall not be discussed here for the sake of brevity. However, the prior ranges can be visually assessed in Figure A1, which shows the resultant posterior probability distribution from the TTV analysis. Figure A1 indicates that the posteriors are multi-modal, suggesting that more than one orbital solution could explain the data.

We found that many pairs of planetary masses and eccentricities and orbital periods of TOI-2202 c provide a plausible explanation of the extracted TTVs of TOI-2202 b. The ambiguity in eccentricity versus dynamical planetary mass was already observed in another K-dwarf TESS system consistent with two warm Jovian-mass planets, TOI-216 (Dawson et al. 2019). Based only on the TTVs for this system, Dawson et al. (2019) were unable to firmly decide whether the system is an eccentric pair of a Saturn-mass planet accompanied by a Neptune-mass planet or is composed of a Jovian-mass planet in a 2:1 MMR with a sub-Saturn-mass planet, where both planets are consistent with more circular orbits. Only after securing a large number of RV measurements and expanding the TTV baseline of TOI-216 did Dawson et al. (2021) confirm the latter configuration. From the posterior probability distribution, we found that TTVs of TOI-2202 b are most likely induced by an exterior sub-Saturn, whose orbital period is close to the first-order eccentricity-type 2:1 MMR with the transiting planet. Such a planetary configuration resembles the solution for TOI-216. Nevertheless, as can be seen from Figure A1, at this stage we cannot completely rule out that the TOI-2202 system is more eccentric and resides in the second-order MMR in the 3:1 period ratio commensurability.

The precise RV measurements that we obtained for TOI-2202 could further constrain the orbital eccentricity and planetary masses, and break the ambiguity. Therefore, as a next step, we included the RV data in the analysis by performing a joint TTV+RV nested sampling analysis, repeating the steps and prior ranges listed above. The RV inclusion in the modeling leads to a more complex model, which now fits the RV semi-amplitude K parameter for each planet, constraining the planetary masses and the RV data offsets and RV jitter parameters for HARPS, PFS, and FEROS, adding six more free parameters. For these, we defined experimentally defined uniform priors of RV off. $\in { \mathcal U }$(–140.0, –70.0) m s⁻¹, and log-uniform (Jeffreys) priors of RV jitter $\in { \mathcal J }$(0.0, 50.0) m s⁻¹.

Figure A2 shows the results from this analysis. The posterior probability distribution of the parameters from this test is definitive, suggesting a Jovian–Saturn pair close to the 2:1 MMR. Yet, some fraction of the samples are samples are consistent with configurations at the 3:1 period ratio commensurability. Figure 5 shows the TTV+RV dynamical fit with maximum $-\mathrm{ln}{ \mathcal L }$ from the nested sampling posteriors, which is consistent with a pair of massive planets close to the 2:1 MMR. For completeness, we examined the alternative best-fit model near the 3:1 commensurability. Figure 6 shows the competing ∼3:1 period ratio best-fit solution. The quality of these fits and overall posterior probability median values and 1σ uncertainties are listed in Table 3. From Figure 5 and Figure 6 it is very clear that a pair close to the 2:1 MMR or the 3:1 MMR can explain the TTV signal. However, the only reasonable solution to the RV data is shown in Figure 5, i.e., the system in the 2:1 commensurability. We note that in both cases shown in Figures 5 and 6 the RV scatter is large, which we attribute to the relatively low S/N of the spectra that we achieved for this rather faint star, but the RV jitter parameter estimates of the HARPS, FEROS, and PFS data are far more reasonable in the 2:1 period ratio case. A visual inspection of Figures 5 and 6 shows that the RV data follows the two-planet dynamical model adequately. The orbital solution close to the 2:1 MMR has $-\mathrm{ln}{ \mathcal L }$ = −136.45, which is statistically more significant than the one close to the 3:1 MMR, which has $-\mathrm{ln}{ \mathcal L }$ = −159.34 (i.e., ${\rm{\Delta }}\mathrm{ln}{ \mathcal L }$ = 22.9). Therefore, we conclude that the combination of TTVs from TESS and CHAT, and the RV data from FEROS, HARPS, and PFS, firmly point to a massive pair of Jovian planets with periods of P_b = ${11.9101}_{-0.0009}^{+0.0009}$ days, and P_c = ${24.754}_{-0.008}^{+0.007}$ days, eccentricities of e_b = ${0.078}_{-0.010}^{+0.014}$ and e_c = ${0.011}_{-0.006}^{+0.009}$, and dynamical masses of m_b = ${0.927}_{-0.036}^{+0.035}$ M_Jup and m_c = ${0.191}_{-0.030}^{+0.032}$ M_Jup.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Exo-Striker Posteriors Probability and Maximum $-\mathrm{ln}{ \mathcal L }$ Orbital Parameters Estimates of the Two-planet System TOI-2202

	∼2:1 MMR Fit				∼3:1 MMR Fit
	Median and 1σ	Max. $-\mathrm{ln}{ \mathcal L }$			Median and 1σ	Max. $-\mathrm{ln}{ \mathcal L }$
Parameter	Planet b	Planet c	Planet b	Planet c	Planet b	Planet c	Planet b	Planet c
K (m s−1)	${94.2}_{-3.6}^{+3.4}$	${15.2}_{-2.4}^{+2.6}$	95.3	19.6	${83.1}_{-7.7}^{+8.0}$	${22.6}_{-2.7}^{+2.6}$	86.9	19.5
P (days)	${11.9108}_{-0.0009}^{+0.0009}$	${24.7545}_{-0.0078}^{+0.0073}$	11.9103	24.7557	${11.9145}_{-0.0010}^{+0.0009}$	${34.6091}_{-0.0125}^{+0.0151}$	11.9164	34.5803
e	${0.0782}_{-0.0097}^{+0.0135}$	${0.0110}_{-0.0064}^{+0.0094}$	0.0672	0.0104	${0.1795}_{-0.0172}^{+0.0184}$	${0.1318}_{-0.0170}^{+0.0210}$	0.2069	0.1308
ω (deg)	${26.9}_{-5.6}^{+5.2}$	${55.7}_{-14.8}^{+16.9}$	44.4	103.24	${109.8}_{-4.5}^{+5.9}$	${198.9}_{-7.8}^{+7.1}$	117.9	205.8
M0 (deg)	${71.2}_{-5.5}^{+6.5}$	${145.8}_{-18.7}^{+13.3}$	51.2	101.7	${332.3}_{-8.3}^{+6.2}$	${202.8}_{-4.7}^{+5.6}$	319.1	200.3
a (au)	${0.0956}_{-0.0016}^{+0.0015}$	${0.1558}_{-0.0026}^{+0.0025}$	0.0957	0.1558	${0.0957}_{-0.0016}^{+0.0015}$	${0.1948}_{-0.0032}^{+0.0031}$	0.0957	0.1947
m (MJup)	${0.927}_{-0.048}^{+0.047}$	${0.191}_{-0.030}^{+0.033}$	0.939	0.246	${0.806}_{-0.081}^{+0.081}$	${0.316}_{-0.039}^{+0.037}$	0.840	0.271
RVoff FEROS (m s−1)	−126.5${}_{-7.1}^{+7.4}$		−125.9		−128.4${}_{-7.3}^{+8.9}$		−123.5
RVoff PFS (m s−1)	−76.3${}_{-3.7}^{+3.6}$		−79.5		−65.5${}_{-11.7}^{+9.1}$		−67.0
RVoff HARPS (m s−1)	−109.7${}_{-5.1}^{+4.7}$		−110.5		−111.6${}_{-7.2}^{+6.2}$		−116.6
RVjit FEROS (m s−1)	${30.7}_{-5.5}^{+7.3}$		29.0		${44.3}_{-6.5}^{+8.4}$		45.2
RVjit PFS (m s−1)	${5.9}_{-5.7}^{+7.0}$		5.4		${20.6}_{-6.7}^{+12.3}$		12.3
RVjit HARPS (m s−1)	${19.4}_{-3.8}^{+4.6}$		17.4		${31.1}_{-5.1}^{+7.1}$		32.0
$-\mathrm{ln}{ \mathcal L }$			−136.448				−159.336

Note. Two possible orbital configurations could explain the observed data: a system close to the first-order 2:1 MMR and a system close to the second-order 3:1 MMR, with a significant statistical preference to the former. The orbital elements are in the Jacobi frame and are valid for epoch BJD = 2458327.103. Only coplanar and edge-on systems (i_b , i_c = 90° and Δi = 0°) are assumed. The joint TTV+RV dynamical model accepts a fixed value of the stellar mass (0.823 M_⊙); however, the derived planetary posterior parameters of a and m are calculated taking into account the stellar mass uncertainty according to the floor uncertainties predicted by Tayar et al. (2020), listed in Table 1.

Download table as:  ASCIITypeset image

The transit light curves and the RV data of TOI-2202 contain the dynamical signature of the gravitationally interacting planets in the system. Therefore, we performed an alternative, more complex, orbital fitting with respect to the analysis presented in Section 4.4. We adopted a self-consistent photodynamical model that fits directly the transit light curves and RVs in an attempt to extract more accurate estimates of the planetary orbital and physical parameters. This comes at the cost of significantly more CPU time.

We chose to apply the photodynamical model on the raw photometric light curves of TOI-2202, which we simultaneously de-trended during the orbital fitting. We included linear models fit to the CHAT data, which simultaneously de-trend the light curves against airmass at the time of observation, and a GP regression model fit to the TESS transit light curves, which aims to capture the evident stellar activity signals seen in the TESS data (see Section 2.1 and Figure 2). The transit GP model parameters we shared with a GP model applied to the RV time series. The inclusion of a complex transit+RV GP model component to the already complex photodynamical model is well justified, since the TESS light curves exhibit periodicity near the estimated orbital period of TOI-2202 c, which could possibly affect its RV signature and, thereafter, our mass and eccentricity estimates of the planets. This is similar to the case of the GJ 143 system (TOI-186, Dragomir et al. 2019; Trifonov et al. 2019a), for which the TESS light curves and spectroscopy data are consistent with a stellar activity periodicity that is very close to the orbital period of the transiting Neptune-mass planet GJ 143 b (Gan et al. 2021). Although we did not detect significant periodic signals in the spectroscopic activity indices and RVs near the orbital frequency TOI-2202 c, the large RV jitter observed in all three RV data sets motivated us to adopt a "common" GP regression model. For this purpose, we adopted the rotational GP regression kernel as formulated by Foreman-Mackey et al. (2017b):

$$\begin{eqnarray}&&k(\tau )=\displaystyle \frac{B}{2+C}{e}^{-\tau /L}\left[\cos \left(\displaystyle \frac{2\pi \tau }{{P}_{\mathrm{rot}}}\right)+(1+C)\right],\end{eqnarray} \tag{ 1 }$$

where P_rot is a proxy for the rotation period of the star, L is the coherence timescale (e.g., the life-time of stellar spots), τ is the time lag between two consecutive data points, and C is a balance parameter for the periodic and the nonperiodic parts of the GP kernel. These parameters were set common to the RV and transit light-curve parts of the model. The parameter B defines the GP co-variance amplitude, and thus, naturally, the RV and the transit models have separate amplitude parameters.

Similar to that described in Section 4.4, we performed a nested sampling test, which we used to estimate the fitted parameter posteriors and confidence intervals, with a numerical time step in the dynamical model set to 0.1 days. We fitted the TESS and the CHAT light curves, adopting a central body mass of 0.823 M_⊙ and free orbital parameters, which for the two planets are the planetary orbital period P, eccentricity e and argument of periastron ω, and orbital inclination i. The dynamical modeling scheme within the Exo-Striker fits the osculating orbital parameters for a given epoch (i.e., nonstatic, perturbed orbital parameters), and thus, it requires the mean-anomaly M₀ parameterization instead of the time of inferior transit conjunction t₀, which is commonly used in Keplerian models when fitting transit light curves. For each sampling iteration, the Exo-Striker computes the perturbed orbital elements and transit times, which allows for the precise modeling of the light curves. Since only TOI-2202 b transits, we only fit the light curves with the planetary scaled semimajor axis, a_b /R_⋆, and radius R_b /R_⋆, whereas for TOI-2202 c these are unconstrained. However, we allowed i_c to vary in the fitting to account for mutual orbital inclinations, and therefore non-coplanar orbital dynamics. Thus, we allowed TOI-2202 c to transit in our orbital analysis, although such models are naturally penalized by a poorer $-\mathrm{ln}{ \mathcal L }$. The fitted orbital parameters for the RV model are shared with those of the transit model, while the RV signal semi-amplitude K parameters for each planet constrain the planetary masses in the dynamical model. For the TESS and the CHAT light curves, we adopted different quadratic limb-darkening models and varied the quadratic limb-darkening parameters u₁ and u₂ for each instrument. We also varied the flux offset and jitter parameter of each transit light curve, and the offset and jitter parameters of each RV data set.

Figure 7 shows the transit component of the photodynamical model constructed together with a GP model for TESS and linear models for CHAT, and constrained by the RV data. Figure 8 shows the RV component of the photodynamical model fitted to the FEROS, the HARPS, and the PFS data, including the GP model component. The final estimates we derived from the joint photodynamical model posterior probability distribution are planetary periods P_b = ${11.9101}_{-0.0036}^{+0.0022}$ days and P_c = ${24.674}_{-0.034}^{+0.026}$ days, eccentricities e_b = ${0.042}_{-0.008}^{+0.025}$ and e_c = ${0.062}_{-0.021}^{+0.026}$, and dynamical masses m_b = ${0.978}_{-0.059}^{+0.063}$ M_Jup and m_c = ${0.369}_{-0.084}^{+0.103}$ M_Jup. The mutual inclination we constrained to Δi = ${6.56}_{-2.10}^{+1.92}$ deg. We note that the posterior distribution for the mass TOI-2202 c planet is bimodal, with peaks at ∼0.30 M_Jup and ∼0.47 M_Jup. We took the median of this bimodal distribution, meaning that it is possible that the mass uncertainties of TOI-2202 c are slightly underestimated. Also, it shows that a full photodynamical model can also suffer from mass-eccentricity ambiguities. The full set of fitted and derived orbital and physical planetary parameter estimates of TOI-2202 b and c are listed in Table 4. The nuisance parameter estimates from the photodynamical nested sampling analysis are listed in Table A4. ³⁷

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Nested Sampling Priors, Posteriors, and Maximum $-\mathrm{ln}{ \mathcal L }$ Orbital Parameters of the Two-planet System TOI-2202

	Median and 1σ		Max. $-\mathrm{ln}{ \mathcal L }$		Adopted Priors
Parameter	Planet b	Planet c	Planet b	Planet c	Planet b	Planet c
K (m s−1)	${99.2}_{-5.2}^{+5.6}$	${29.3}_{-6.6}^{+8.3}$	92.9	38.3	${ \mathcal U }$(90.0, 110.00)	${ \mathcal U }$(20.0, 40.0)
P (days)	${11.9101}_{-0.0036}^{+0.0022}$	${24.6744}_{-0.0339}^{+0.0258}$	11.9123	24.6797	${ \mathcal U }$(11.90, 11.92)	${ \mathcal U }$(24.60, 24.80)
e	${0.0420}_{-0.0075}^{+0.0255}$	${0.0622}_{-0.0211}^{+0.0452}$	0.0672	0.0104	${ \mathcal U }$(0.00, 0.2)	${ \mathcal U }$(0.00, 0.2)
ω (deg)	${84.1}_{-16.4}^{+9.8}$	${320.8}_{-12.5}^{+129.7}$	86.0	322.6	${ \mathcal U }$(0.0, 360.00)	${ \mathcal U }$(0.0, 360.00)
M0 (deg)	${7.6}_{-8.5}^{+26.4}$	${267.2}_{-269.9}^{+9.6}$	4.3	265.8	${ \mathcal U }$(0.0, 360.00)	${ \mathcal U }$(0.0, 360.00)
λ (deg)	${90.4}_{-0.8}^{+14.6}$	${226.4}_{-15.5}^{+7.1}$	90.3	228.4	(derived)	(derived)
i (deg)	${88.4}_{-3.3}^{+0.6}$	${84.7}_{-2.9}^{+2.4}$	88.9	87.2	${ \mathcal U }$(80.0, 90.0)	${ \mathcal U }$(80.0, 90.0)
Ω (deg)	0.0	${5.0}_{-2.8}^{+2.3}$	0.0	5.9	(fixed)	${ \mathcal U }$(0.0, 10.0)
a/R⋆	${23.08}_{-2.18}^{+2.90}$	...	23.45	...	${ \mathcal U }$(20.0, 35.00)	...
R/R⋆	${0.1261}_{-0.0065}^{+0.068}$	...	0.1255	...	${ \mathcal U }$(0.01, 0.25)	...
Δi (deg)	${6.56}_{-2.10}^{+1.92}$	...	6.14	...	(derived)	...
a (au)	${0.09564}_{-0.00161}^{+0.00156}$	${0.15544}_{-0.00263}^{+0.00255}$	0.09569	0.15554	(derived)	(derived)
m (MJup)	${0.978}_{-0.0588}^{+0.0630}$	${0.369}_{-0.0836}^{+0.103}$	0.917	0.482	(derived)	(derived)
R (RJup)	${1.01}_{-0.080}^{+0.522}$	...	0.992	...	(derived)	...

Note. Additionally, these parameters were estimated including a Gaussian processes regression kernel, which was used to model the stellar rotation effects and was common to the TESS light curves and RV data. CHAT light curves were simultaneously modeled with linear regression models to account for airmass optical effects. These and other nuisance parameter estimates are listed in Table A4. The orbital elements are in the Jacobi frame and are valid for epoch BJD = 2458327.103. The joint dynamical model accepts a fixed value for the stellar mass (0.823 M_⊙); however, the derived planetary posterior parameters of a, m, and R are calculated taking into account the stellar parameter uncertainties (see the Note in Table 2). The median value of m_c comes from a bimodal distribution (see Figure 9).

Download table as:  ASCIITypeset image

We performed a long-term stability and dynamical analysis of the TOI-2202 system using a custom version of the Wisdom–Holman N-body algorithm (Wisdom & Holman 1991), which directly adopts and integrates the Jacobi orbital elements from the Exo-Striker. Because of the relatively short orbital periods of the TOI-2202 planets, we chose a very small integration time step equal to 0.02 days, which was necessary for the accurate numerical calculation analysis of the dynamical properties of the system. However, this short time step leads to severe numerical overhead, which limited our ability to test the system over the system's estimated age. For example, an N-body integration of the TOI-2202 system for only 10 Myr takes approximately one day on a modern CPU. An educated guess, however, suggests that longer integration times are likely not necessary. The Exo-Striker provides instant Angular Momentum Deficiency (AMD, Laskar & Petit 2017) monitoring, which indicates that given the estimated small orbital eccentricities, semimajor axes, mutual orbital inclination, and planetary masses, the TOI-2202 system is AMD-stable. For more details on the calculation of the AMD-stability criteria, we refer to Laskar & Petit (2017). Additionally, given the estimated semimajor axes and masses in the system, we can calculate the mutual Hill radii of the planets as

$$\begin{eqnarray}&&{R}_{\mathrm{Hill},{\rm{m}}}\approx \sqrt[3]{\displaystyle \frac{({m}_{b}+{m}_{c})}{3{M}_{\star }}}\displaystyle \frac{({a}_{b}+{a}_{c})}{2}\approx 0.01\,\mathrm{au}\end{eqnarray} \tag{ 2 }$$

and since a_c –a_b = 0.06 au, then the planets have about 6 R_Hill,m, which is above the ∼3.5 R_Hill,m threshold needed for the system to be considered Hill-stable (see, Gladman 1993). In terms of AMD and Hill stability, the TOI-2202 planetary system must be generally stable despite the close planetary orbits. Nevertheless, the AMD and Hill stability criteria do not account for the system's dynamics near mean motion resonances, thus they can only be used as a proxy for long-term stability. Therefore, as a compromise, our long-term stability simulations of the TOI-2202 system were performed for a maximum of 1 Myr by numerically integrating 10,000 randomly chosen samples of the achieved parameter posteriors from the photodynamical modeling scheme. We ran our stability test on a modern 40 CPU Intel Xeon based workstation, which took about four weeks to complete, and which we find to be reasonable for the N-body dynamical analysis in this work.

Given the period ratio of the system close to the 2:1 commensurability, we inspected the first-order MMR angles θ₁ and θ₂ , which are defined as:

$$\begin{eqnarray}&&{\theta }_{1}={\lambda }_{b}-2{\lambda }_{c}+{\varpi }_{b},\ {\theta }_{2}={\lambda }_{b}-2{\lambda }_{c}+{\varpi }_{c},\end{eqnarray} \tag{ 3 }$$

where ϖ_b,c = Ω_b,c + ω_b,c are the planetary longitudes of periatron and λ_b,c = M_0b,c + ϖ_b,c are the mean longitudes. We also monitored for libration of the secular apsidal angle Δω, which is defined as:

$$\begin{eqnarray}&&{\rm{\Delta }}\omega ={\theta }_{1}-{\theta }_{2}={\varpi }_{b}-{\varpi }_{c},\end{eqnarray} \tag{ 4 }$$

which indicates if the dynamics of the system is dominated by secular interactions, exhibiting apsidal libration in alignment (Δω librating around 0°), anti-alignment (Δω librating around 180°), or an asymmetric libration.

The results from our long-term stability analysis indicate that all examined 10,000 samples are stable for 1 Myr with very similar dynamical behavior. Figure 9 shows the derived posteriors of the dynamical properties of the studied 10,000 samples. The distribution of dynamical parameters reveals low-eccentricity evolution, but despite being close to the 2:1 MMR, the TOI-2202 pair seems to reside outside of the low-order eccentricity-type 2:1 MMR. The mean period ratio evolution is osculating around 2.07, while we did not detect libration of the resonance angles θ₁, θ₂, and the apsidal alignment angle Δω. The posterior distributions of θ₁, θ₂, and Δω are consistent with circulation, with libration amplitudes between 0° and 360°.

Zoom In Zoom Out Reset image size

Figure 10 shows a 200 yr extent of the dynamical evolution of the photodynamical fit with a maximum $-\mathrm{ln}{ \mathcal L }$ from the posteriors probability samples (see Table 4), which is representative of the overall dynamics of the tested posterior samples. Figure 10 shows the evolution of the semimajor axes a_b and a_c, eccentricities e_b and e_c, the period ratio, the apsidal alignment angle Δϖ, and the characteristic 2:1 MMR angles θ₁ and θ₂. No resonance angle libration is observed. For the same fit, Figure 11 shows the trajectory evolution of different combinations of e_b , e_c , and sine and cosine functions of the 2:1 MMR resonance angles θ₁ and θ₂. There is no observed libration around a fixed point in the trajectory evolution that could suggest a 2:1 MMR. Similar trajectory evolution is observed in the majority of the posterior samples that are within the 1σ credible interval.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The eccentricities of both planets are low and when that is the case, the resonant and near-resonant dynamics can be studied analytically following Nesvorný & Vokrouhlický (2016). There are three variables to consider, each of them being a combination of orbital elements. Constant δ is an orbital invariant that defines the position of the system relative to the 2:1 commensurability, the resonant angle ψ is a combination of θ₁ and θ₂, and variable Ψ is a combination of planetary masses, semimajor axes, and eccentricities (see Nesvorný & Vokrouhlický 2016 for details). Figure 12 shows the position of TOI-2202 in the context of the dynamical variables δ, ψ, and Ψ. The resonant librations of ψ can only happen for δ > 0.945. The best-fit and median values of TOI-2202, listed in Table 4, lead to δ ≃ −0.532 and −0.77, respectively, and the system is therefore firmly outside the libration region. For TOI-2202, Ψ and ψ follow a deformed circle that is slightly offset from the origin, which means that ψ circulates and Ψ oscillates (0.4 ≲ Ψ ≲ 1.1).

Zoom In Zoom Out Reset image size

Analytic expressions can be used to relate the TTVs measured for TOI-2202 to the system's architecture. Given that the two orbits are not librating in the 2:1 resonance, here we use analytic TTV expressions given in Agol & Deck (2016). Equation (6) in that paper gives the nonresonant TTVs to the first order in orbital eccentricities. We only consider the highest amplitude term for the inner planet:

$$\begin{eqnarray}&&\delta {t}_{1}=\displaystyle \frac{{P}_{b}}{2\pi }\displaystyle \frac{{m}_{c}}{{M}_{* }}{f}_{1,1}^{(-2)}{e}_{c}\sin {\theta }_{2},\end{eqnarray} \tag{ 5 }$$

where ${f}_{1,1}^{(-2)}\simeq 0.05$. The expected amplitude of TTVs is therefore A_TTV ≃ 0.05 days, which compares well with the measured amplitude (see Figure 5). Other terms, including TTV chopping during planet conjunctions (Nesvorný & Vokrouhlický 2014; Deck & Agol 2015), are responsible for the slight deviations of measured TTVs from a perfect sinusoid. These terms are important to break degeneracies in the TTV inversion. The expected TTV period is the period of the θ₂ angle. Neglecting the long-term evolution of ϖ_c , the TTV period can be approximated by the superperiod:

$$\begin{eqnarray}&&{P}_{\mathrm{TTV}}={\left[\displaystyle \frac{1}{{P}_{b}}-\displaystyle \frac{2}{{P}_{c}}\right]}^{-1}.\end{eqnarray} \tag{ 6 }$$

Adopting the best-fit orbital periods from Table 4, we obtain P_TTV ≃ 340 days, which is a good match to Figure 5.