Two Jovian Planets around the Giant Star HD 202696: A Growing Population of Packed Massive Planetary Pairs around Massive Stars?

The HIRES RV time series of HD 202696 consists of 42 precise Doppler measurements taken between 2007 July and 2014 September. With a median RV precision of 1.37 m s⁻¹ and an end-to-end amplitude of 127 m s⁻¹, these data clearly show systematic motions, resembling two Keplerians caused by two massive orbiting planets. Figure 3 shows the precise RV data of HD 202696, together with our best two-planet Keplerian model and its uncertainties, while Figure 4 shows the two fairly well-sampled Keplerian signals phase-folded at their best-fit periods (see Section 3.2 for a quantitative analysis of the data). Note that the data uncertainties displayed in Figures 3 and 4 already contain the best-fit RV jitter estimate of ∼8 m s⁻¹, which was automatically added in quadrature to the error budget during the fitting. This RV jitter estimate is typical for early K giants such as HD 202696 and is likely due to solar-like oscillations, which have periods of less than a day and are thus not resolved in the long-term RV time series. The scaling relations from Kjeldsen & Bedding (2011) predict 3.7 m s⁻¹ RV jitter for HD 202969, while the typical RV jitter value for a star with a B − V color of 1.003 mag in the Lick sample is 10–20 m s⁻¹ (see Figure 3 in Trifonov et al. 2014). Since the stars in the Lick sample are, on average, slightly more evolved than HD 202696, it might be reasonable to expect HD 202696 to have slightly smaller RV jitter than most other stars in this sample.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is worth noting that with the release of the HIRES database, Butler et al. (2017) reported the presence of at least three periods in the Keck data of HD 202696. They identified two significant signals with periods of 522.3 ± 6.4 and 946.0 ± 19.0 days, classified as planetary "candidates," and an additional period of 214.2 ± 4.8 days, which required further confirmation.

As an independent test, we derive the generalized Lomb–Scargle periodogram (GLS; Zechmeister & Kürster 2009) of the HIRES RV data, as well as that of the S- and H-index activity indicators (Figure 5) also provided by Butler et al. (2017). These time series are tabulated in Table 2. For reference, in all panels of Figure 5, the range of possible stellar rotation frequencies of HD 202696 is indicated by the red shaded area, with the red dashed line at the most likely value (1/${P}_{\mathrm{rot}.}={0.00615}_{-0.00177}^{+0.00418}$ day⁻¹; see Section 2). The horizontal dotted, dot-dashed, and dashed lines correspond to false-alarm probability (FAP) levels of 10%, 1%, and 0.1%, respectively. We use an FAP <0.1% as our significance threshold. The top panel of Figure 5 reveals a dominant signal around 540 days (blue dashed line) in the HIRES RVs, which is significant. The second panel of Figure 5 shows the periodogram of the RV data after the ∼540 day signal has been subtracted; it indicates a second significant signal with a period of ∼970 days (green dashed line). The middle panel of Figure 5 shows that no period with an FAP <0.1% can be detected after the two most significant signals have been consecutively subtracted from the data.

Zoom In Zoom Out Reset image size

The ∼214 day signal previously reported by Butler et al. (2017) is also visible in the residuals of the data, but with 1% < FAP < 10%, which we consider insignificant. Moreover, the 214 day signal falls within the red shaded area of possible stellar rotation frequencies and thus is not considered further in our analysis. No significant periodicities can be identified in the two panels of Figure 5, which show the periodograms of the H- and S-index activity indicators, respectively.

Additionally, we inspected the V-band measurements of HD 202696 from the All-Sky Automated Survey (ASAS; Pojmanski et al. 2005). We find only 48 high-quality ASAS V-band measurements (graded A and B with $\overline{{\sigma }_{{\rm{V}}}}=0.037$ mag), mostly taken between 2003 May and 2007 October, i.e., prior to the HIRES Doppler observations of HD 202696. Thus, the sparse ASAS photometry is only relevant to check for evident intrinsic photometric variability, which could be associated with the RV signals. No significant periodicity is observed in any of the five available ASAS apertures. The bottom panel of Figure 5 shows the GLS periodogram of the averaged ASAS photometry measurements, indicating that HD 202696 is most likely a photometrically stable star.

The two significant signals at ∼540 and ∼970 days are consistent with those reported by Butler et al. (2017) as planet candidates. We find that these signals are not consistent with the range of possible stellar rotation frequencies of HD 202696. Also, the S- and H-index activity indicators of the HIRES data and ASAS photometry do not show any significant periodic structure or correlation that could be associated with the RV signals. Thus, we consider the ∼540 and ∼970 day signals as robust planetary detections.

In order to determine the orbital parameters of the HD 202696 b and c planetary companions, we adopt a maximum-likelihood estimator (MLE) scheme coupled with a Nelder–Mead algorithm (also known as a downhill simplex algorithm; see Nelder & Mead 1965; Press et al. 1992). In our scheme, we minimize the negative logarithm of the model's likelihood function ($-\mathrm{ln}{ \mathcal L }$), while we simultaneously optimize the planetary RV semi-amplitudes K_b,c, periods P_b,c, eccentricities e_b,c, arguments of periastron ω_b,c, mean anomalies M_b,c, and the RV zero-point offset. Prior information on the period, phase, and amplitude for each of the two-planet candidates is taken from the GLS periodogram test (see Figure 5), which provides good starting parameters for the MLE algorithm. All parameters are derived for a uniform reference epoch of 2,450,000 (BJD). Following Baluev (2009), we also include the RV jitter⁹ as an additional parameter in the modeling.

We apply two models to the RV data: a standard unperturbed two-planet Keplerian model and a more accurate self-consistent N-body dynamical model, which models the gravitational interactions between the planets by integrating the equations of motion using the Gragg–Bulirsch–Stoer integration method (Press et al. 1992). We cross-check the obtained best-fit results for consistency with the Systemic Console 2.2 package (Meschiari et al. 2009) and find excellent agreement between our MLE code¹⁰ and Systemic.

For parameter distribution analysis and uncertainty estimates, we couple our MLE fitting algorithm with a Markov chain Monte Carlo (MCMC) sampling scheme using the emcee sampler (Foreman-Mackey et al. 2013). For all parameters, we adopt flat priors (i.e., equal probability of occurrence), and we run emcee from the best fit obtained by the MLE. We select the 68.3% confidence levels of the posterior MCMC parameter distribution as 1σ parameter uncertainties.

The best-fit parameters of our two-planet Keplerian and dynamical models, together with their 1σ MCMC uncertainties, are tabulated in Table 3. We first fit the Keplerian model, and as a next step, we adopt the best two-planet Keplerian parameters as an input for the more sophisticated dynamical model. For consistency with the unperturbed Keplerian frame and in order to work with minimum dynamical masses, we assume an edge-on and coplanar configuration for the HD 202696 system (i.e., i_b,c = 90° and ΔΩ = 0°). We find that the models are practically equivalent, suggesting a pair of planets with equal masses of the order of ∼2.0 M_jup and a period ratio of ${P}_{\mathrm{rat}.}\approx 1.83$ near the 11:6 MMR (for details, see Table 3).

Table 2.  HIRES Doppler Measurements for HD 202696 from Butler et al. (2017) and RV Residuals of the Best-fit Models from Table 3

Epoch (BJD)	RV (m s−1)	${\sigma }_{\mathrm{RV}}$ (m s−1)	S Index	H Index	o−cKep. (m s−1)	o−cN-body (m s−1)
2,454,287.927	2.88	1.96	0.1091	0.0312	−2.94	−3.57
2,454,399.822	−63.90	1.25	0.1072	0.0307	−2.53	−4.50
2,454,634.038	−20.49	1.14	0.0899	0.0312	−5.05	−5.41
2,454,674.900	2.57	1.39	0.1066	0.0308	−0.79	−0.85
2,454,718.006	10.64	1.49	0.0868	0.0308	−4.21	−4.00
2,454,777.857	11.66	1.19	0.1132	0.0309	−4.38	−3.84
2,454,778.819	13.33	1.39	0.1065	0.0308	−2.61	−2.07
2,454,805.740	10.92	1.40	0.1138	0.0308	−1.13	−0.41
2,454,955.107	−5.77	1.32	0.1055	0.0308	16.03	16.72
2,454,964.128	−15.69	1.23	0.1075	0.0310	7.25	7.84
2,454,984.066	−23.92	1.38	0.1065	0.0311	0.58	0.93
2,454,985.007	−26.07	1.40	0.1114	0.0309	−1.53	−1.20
2,455,014.963	−39.57	1.31	0.1072	0.0308	−15.66	−15.70
2,455,015.957	−41.36	1.31	0.1061	0.0310	−17.54	−17.59
2,455,019.021	−25.42	1.30	0.1143	0.0309	−1.89	−1.97
2,455,043.881	−18.51	1.40	0.1074	0.0307	1.18	0.87
2,455,076.033	−9.97	1.37	0.1129	0.0306	1.10	0.74
2,455,084.036	0.00	1.61	0.1251	0.0306	8.42	8.08
2,455,106.909	5.94	1.47	0.1086	0.0307	6.10	5.92
2,455,133.918	17.65	1.40	0.0936	0.0305	7.88	7.96
2,455,169.772	21.46	1.40	0.1145	0.0307	1.14	1.44
2,455,170.691	17.93	1.43	0.1069	0.0308	−2.59	−2.29
2,455,187.695	20.55	1.39	0.1011	0.0308	−2.85	−2.57
2,455,198.700	21.99	1.42	0.1075	0.0309	−2.35	−2.15
2,455,290.146	2.47	1.33	0.1010	0.0309	0.17	−1.29
2,455,373.126	−33.29	1.23	0.1137	0.0309	7.89	7.22
2,455,396.059	−52.38	1.36	0.1097	0.0309	−1.22	−1.26
2,455,436.743	−69.28	1.27	0.0975	0.0307	−5.03	−4.34
2,455,455.739	−73.42	1.34	0.1085	0.0307	−5.25	−4.49
2,455,521.793	−53.00	1.37	0.1062	0.0308	16.87	16.67
2,455,698.124	−10.58	1.35	0.1006	0.0309	−16.91	−17.80
2,455,770.042	53.71	1.56	0.1099	0.0309	28.09	27.74
2,455,787.946	20.04	9.66	0.0659	0.0285	−6.39	−6.71
2,455,841.869	26.89	1.26	0.0882	0.0308	5.86	5.65
2,455,880.836	8.95	1.23	0.1292	0.0312	−3.05	−3.06
2,455,903.717	−1.19	1.25	0.1071	0.0309	−6.85	−6.68
2,455,931.693	2.89	1.29	0.0945	0.0307	5.29	5.69
2,456,098.123	−28.24	1.31	0.1084	0.0309	−8.23	−7.59
2,456,154.011	−2.81	1.46	0.1084	0.0308	4.67	5.64
2,456,194.986	3.05	1.31	0.1098	0.0307	2.08	3.36
2,456,522.092	−57.44	1.43	0.1149	0.0309	1.73	0.45
2,456,894.081	19.32	1.23	0.1113	0.0309	−5.09	−5.03

Download table as:  ASCIITypeset image

Both the Keplerian and dynamical fits reveal a moderate best-fit value of e_c for the eccentricity of the outer planet but with a large uncertainty toward nearly circular orbits. This indicates that e_c is likely poorly constrained, and perhaps this moderate eccentricity value is a result of overfitting. Indeed, a two-planet Keplerian model with circular orbits (i.e., e_b,c, ω_b,c = 0) has $-\mathrm{ln}{ \mathcal L }=150.286$, meaning that the difference between the circular and noncircular Keplerian models is only ${\rm{\Delta }}\mathrm{ln}{ \mathcal L }$ ≈ 1.5, which is insignificant.¹¹ Thus, we conclude that neither the Keplerian nor the dynamical model applied to the current RV data can place tight constraints on the orbital eccentricities.

The inclusion of the planetary inclinations i_b,c and ΔΩ (i.e., ${\rm{\Delta }}i$) in the dynamical modeling was also not justified. We find that all dynamical fits with coplanar inclinations in the range of 5° < i_b,c < 90° are practically indistinguishable from the coplanar and edge-on configuration (i_b,c = 90° and ΔΩ = 0°). Fits with i_b,c < 5° are of poorer quality due to the increased dynamical masses of the companions and the strong gravitational interactions, which are not consistent with the data anymore. Therefore, any further attempt to constrain the planetary inclinations as coplanar or mutually inclined leads to inconclusive results.

Overall, our orbital parameter estimates show that taking the planetary gravitational interactions into account in the fitting does not yield a significant improvement over the Keplerian model. Yet we decide to rely on the dynamical modeling, since this scheme provides physical constraints on the derived orbital parameters and naturally penalizes highly unstable (unrealistic) orbital solutions.

The long-term stability and dynamical properties of the HD 202696 system are tested using the SyMBA N-body symplectic integrator (Duncan et al. 1998), which was modified to work in Jacobi coordinates (e.g., Lee & Peale 2003). We integrate each MCMC sample for a maximum of 1 Myr with a time step of 4 days. In case of planet–planet close encounters, however, SyMBA automatically reduces the time step to ensure an accurate simulation with high orbital resolution. SyMBA also checks for planet–planet or planet–star collisions or planetary ejections and interrupts the integration if they occur. A planet is considered lost and the system unstable if, at any time: (i) the mutual planet–planet separation is below the sum of their physical radii (assuming Jupiter mean density), i.e., the planets undergo collision; (ii) the star–planet separation exceeds two times the initial semimajor axis of the outermost planet (${r}_{\max }\gt 2\times {a}_{{\rm{c}}\mathrm{in}.}$), which we define as planetary ejection; or (iii) the star–planet separation is below the physical stellar radius (R ≈ 0.03 au), which we consider a collision with the star. All of these events are associated with large planetary eccentricities leading to crossing orbits, close planetary encounters, rapid exchange of energy and angular momentum, and eventually instability. Therefore, these somewhat arbitrary stability criteria are efficient to detect unstable configurations and save CPU time.

Figure 6 shows the posterior MCMC distribution of the fitted parameters with a dynamical modeling scheme whose orbital configuration is edge-on, coplanar, and stable for at least 1 Myr. The histogram panels in Figure 6 provide a comparison between the probability density distribution of the overall MCMC samples (black) and the stable samples (red) for each fitted parameter. The two-dimensional parameter distribution panels represent all possible parameter correlations with respect to the best dynamical fit from Table 2, whose position is marked with blue lines. The black 2D contours are constructed from the overall MCMC samples (gray) and indicate the 68.3%, 95.5%, and 99.7% confidence levels (i.e., 1σ, 2σ, and 3σ). For clarity, in Figure 6, the stable samples are overplotted in red.

Zoom In Zoom Out Reset image size

We find that ∼46% of the MCMC samples are stable. The histograms in Figure 6 clearly show that the distributions of all MCMC samples and their stable counterparts are in mutual agreement, peaking near the position of the best dynamical fit. The main exception is the distribution of e_c, which is bimodal with a stronger peak suggesting nearly circular orbits and a smaller peak at somewhat moderate eccentricities from where the best dynamical fit originates. Here P_b and ω_b also show bimodal distributions, which are related to the bimodal distribution of e_c. As we discussed in Section 3.2, however, the best-fit eccentricities are statistically insignificant, and a simpler two-planet model with circular orbits is equally good in terms of $-\mathrm{ln}{ \mathcal L }$. This indicates that the best-fit solution is likely a global $-\mathrm{ln}{ \mathcal L }$ minimum, which lies in a smaller and isolated orbital parameter space with moderate e_c. This minimum and its surroundings, however, are highly unstable and therefore unlikely to represent the HD 202696 system configuration. Thus, running an MCMC test reveals a larger, stable, and statistically confident orbital parameter space toward e_b,c → 0, where the stable fraction of the sample distribution is ∼75%. From the long-term dynamical analysis of the posterior MCMC samples, we can place a boundary at planetary eccentricities e_b,c ≤ 0.1, while larger eccentricities are very unlikely. The overall posterior distribution of ω_b,c and M_b,c are consistent with the best-fit estimate, but the stability test cannot constrain these parameters. Because of the low eccentricities of the stable solutions, these parameters are ambiguous and can be found anywhere between 0° and 360°. However, the orbital mean longitudes λ_b,c = ω_b,c + M_b,c are rather well constrained to ∼±50° (see Table 3).

Given the stability constraints from our MCMC test, the most realistic configuration of the HD 202696 system is represented by the parameter distribution of the stable fits. In Table 4, we present our final orbital estimates of the two-planet system assuming coplanar and edge-on (i_b,c = 90°, ΔΩ = 0°) configurations. These parameters are adopted from the maximum value of the stable parameter probability density function, while the uncertainties are the 0.16 and 0.84 quantiles of the distribution. Our new orbital estimate is consistent with circular orbits ${e}_{{\rm{b}}}={0.011}_{-0.011}^{+0.078}$ and ${e}_{{\rm{c}}}={0.029}_{-0.012}^{+0.065}$ but otherwise consistent with the best two-planet dynamical fit shown in Table 3.

Table 3.  Best Keplerian and Dynamical Fits for the HD 202696 System

Two-planet Keplerian Model
Parameter	HD 202696 b	HD 202696 c
Semi-amplitude K [m s−1]	${31.7}_{-2.1}^{+4.8}$	${28.1}_{-5.3}^{+1.1}$
Period P [days]	${521.0}_{-7.3}^{+6.2}$	${956.1}_{-29.9}^{+22.4}$
Eccentricity e	${0.056}_{-0.040}^{+0.066}$	${0.261}_{-0.202}^{+0.060}$
Arg. of periastron ω [deg]	${259.2}_{-89.2}^{+106.0}$	${129.1}_{-33.6}^{+142.7}$
Mean anomaly M0 [deg]	${69.7}_{-115.8}^{+85.8}$	${152.0}_{-74.3}^{+112.8}$
Mean longitude λ [deg]	${328.8}_{-44.9}^{+50.4}$	${281.1}_{-58.9}^{+50.5}$
Semimajor axis a [au]	${1.573}_{-0.015}^{+0.012}$	${2.357}_{-0.049}^{+0.037}$
Minimum mass $m\sin i$ $[{M}_{\mathrm{Jup}}]$	${1.930}_{-0.132}^{+0.285}$	${2.028}_{-0.354}^{+0.073}$
RV offset [m s−1]	$-{15.32}_{-1.15}^{+2.59}$
RV jitter [m s−1]	${8.20}_{-0.46}^{+3.07}$
$-\mathrm{ln}{ \mathcal L }$	148.799
Two-planet Dynamical Model
Parameter	HD 202696 b	HD 202696c
Semi-amplitude K [m s−1]	${32.5}_{-1.6}^{+3.1}$	${26.6}_{-3.5}^{+1.1}$
Period P [days]	${528.2}_{-12.5}^{+1.8}$	${958.0}_{-26.4}^{+2.2}$
Eccentricity e	${0.026}_{-0.006}^{+0.100}$	${0.265}_{-0.229}^{+0.005}$
Arg. of periastron ω [deg]	${289.7}_{-91.8}^{+35.1}$	${140.0}_{-44.9}^{+88.7}$
Mean anomaly M0 [deg]	${70.2}_{-40.8}^{+67.9}$	${155.5}_{-52.2}^{+70.41}$
Mean longitude λ [deg]	${359.9}_{-63.0}^{+8.2}$	${259.5}_{-68.5}^{+1.0}$
Semimajor axis a [au]	${1.587}_{-0.025}^{+0.004}$	${2.360}_{-0.040}^{+0.003}$
Min. dyn. mass m $[{M}_{\mathrm{Jup}}]$	${1.991}_{-0.011}^{+0.174}$	${1.917}_{-0.212}^{+0.116}$
Inclination i [deg]	90.0 (fixed)	90.0 (fixed)
Node ${\rm{\Delta }}{\rm{\Omega }}$ [deg]	0.0 (fixed)
RV offset [m s−1]	$-{15.12}_{-1.44}^{+1.95}$
RV jitter [m s−1]	${8.31}_{-0.47}^{+3.08}$
$-\mathrm{ln}{ \mathcal L }$	148.973

Note. The parameters are valid for BJD = 2,450,000 and are derived using our stellar mass estimate of M = 1.91 M_⊙. The uncertainties of ${a}_{{\rm{b}},{\rm{c}}}$, ${m}_{{\rm{b}},{\rm{c}}}\sin i$ (Keplerian model), and ${m}_{{\rm{b}},{\rm{c}}}$ (dynamical model) do not take into account the small uncertainty of the stellar mass.

Download table as:  ASCIITypeset image

Table 4.  Mode of the Dynamically Stable MCMC Distribution for the HD 202696 System (Peak of the Red Histograms from Figure 6)

Adopted Two-planet Coplanar Edge-on Configuration Based on Stability
Parameter	HD 202696 b	HD 202696c
Semi-amplitude K [m s−1]	${34.1}_{-2.9}^{+2.3}$	${25.3}_{-3.3}^{+2.1}$
Period P [days]	${517.8}_{-3.9}^{+8.9}$	${946.6}_{-20.9}^{+20.7}$
Eccentricity e	${0.011}_{-0.011}^{+0.078}$	${0.028}_{-0.012}^{+0.065}$
Arg. of periastron ω [deg]	${201.6}_{-98.5}^{+97.8}$	${136.8}_{-62.7}^{+143.0}$
Mean anomaly M0 [deg]	${79.2}_{-28.1}^{+198.6}$	${180.0}_{-108.9}^{+98.9}$
Mean longitude λ [deg]	${298.8}_{-14.6}^{+70.0}$	${262.8}_{-53.3}^{+37.8}$
Semimajor axis a [au]	${1.566}_{-0.007}^{+0.016}$	${2.342}_{-0.035}^{+0.034}$
Dyn. mass m $[{M}_{\mathrm{Jup}}]$	${1.996}_{-0.100}^{+0.220}$	${1.864}_{-0.227}^{+0.177}$
Inclination i [deg]	90.0 (fixed)	90.0 (fixed)
Node ${\rm{\Delta }}{\rm{\Omega }}$ [deg]	0.0 (fixed)
RV offset [m s−1]	$-{14.42}_{-2.22}^{+1.70}$
RV jitter [m s−1]	${9.36}_{-0.57}^{+2.03}$

Note. Same comments as in Table 3.

Download table as:  ASCIITypeset image

Integrating this configuration for 100 Myr reveals a stable planetary system with a mean period ratio $\overline{{P}_{\mathrm{rat}.}}$ ≈ 1.83 and orbital eccentricities oscillating in opposite phase in the range from 0 to 0.05 for e_b and 0 to 0.04 for e_c. This configuration librates mostly in an anti-aligned geometry with a secular apsidal angle ${\rm{\Delta }}\omega ={\omega }_{{\rm{b}}}-{\omega }_{{\rm{c}}}$ ∼ 180° exhibiting a large libration amplitude of ±80°.

The long-term stability of the orbital estimate shown in Table 4 boosts the confidence that the most likely orbital configuration of the HD 202696 system has been identified.

While the dynamical modeling of the HIRES Doppler measurements of HD 202696 was unable to qualitatively constrain the masses of HD 202696 b and c, we also use the long-term stability of fits with forced inclination as an efficient way to limit the range of the possible planetary masses further. Figure 7 shows the results from this test. We inspect the range of inclinations from i = 90° to 5°. To ensure enhanced stability, we adopt the planetary eccentricity estimates from Table 4, which, for each adopted i, remained fixed in the fit. The small ${\rm{\Delta }}\mathrm{ln}{ \mathcal L }$ of the dynamical models in the range of i = 90°–15° shows that these models are statistically equivalent, but the fit quality rapidly deteriorates for i < 15° (not shown in Figure 7). This is likely due to the significantly increased planetary masses, which may alter the system stability even at lower eccentricities.

Zoom In Zoom Out Reset image size

Inclination i = 20° seems to be a stability boundary of the coplanar case, leaving the companion masses in the planetary regime with a maximum of about three times the minimum planetary masses of the coplanar edge-on case. However, since orbital inclinations near i = 90° are statistically more likely, we limit our dynamical analysis of the HD 202696 system to the coplanar and edge-on configuration.

Overall dynamical analysis of the stable samples reveals that the mean period ratio over 1 Myr of integrations is $\overline{{P}_{\mathrm{rat}.}}={1.83}_{-0.05}^{+0.04}$, which is consistent with the period ratio from which the system is initially integrated, ${P}_{\mathrm{rat}.\mathrm{in}.}={1.82}_{-0.03}^{+0.03}$. This implies that the stable MCMC orbital configurations remain regular and well separated at any given time of the stability test. From these stable samples, ∼29% are found with an aligned apsidal libration (${\rm{\Delta }}\omega \sim 0^\circ$), ∼34% in an anti-aligned apsidal libration (${\rm{\Delta }}\omega \sim 180^\circ ;$ similar to the orbital evolution of the long-term stable configuration presented in Table 4), and ∼37% in apsidal circulation (Δω in the range 0°–360°). In fact, the latter case usually involves e_b decreasing from near its maximum value to near zero and e_c increasing from near zero to near its maximum value at ${\rm{\Delta }}\omega \approx 90^\circ$, and vice versa at ${\rm{\Delta }}\omega \approx 270^\circ$, with episodes of rapid circulation when e_b or e_c is nearly zero (see right panel of Figure 9). This type of evolution is consistent with a phase space dominated by large libration islands about Δω = 0° and 180° and having only a narrow region of circulation (see, e.g., Lee & Peale 2003).

Figure 8 shows the distribution of the mean period ratio $\overline{{P}_{\mathrm{rat}.}}$ and the libration amplitudes of Δω for the stable dynamical samples with aligned apsidal libration Δω ∼ 0° (left panel) and anti-aligned apsidal libration Δω ∼ 180° (right panel). From Figure 8, it is clear that both alignment and anti-alignment have similar libration amplitudes of Δω peaking around 70° and with a sharp border at 90°, beyond which the system is involved in apsidal circulation. An interesting pileup in the posterior distribution of $\overline{{P}_{\mathrm{rat}.}}$ can be seen around 1.80, 1.83, and 1.86, which correspond to the 9:5, 11:6, and 13:7 MMRs. We inspected these possible high-order MMRs for librating resonance angles,¹² but we did not detect any. Despite the apsidal libration of the Δω in ∼60% of the configurations, which is suggestive of resonant behavior, the system seems to be out of any MMR. In the Δω ∼ 180° case, a small number of system configurations are consistent with a 2:1 period ratio, but these samples have statistically poorer $-\mathrm{ln}{ \mathcal L }$ and are thus unlikely to represent the HD 202696 system.

Zoom In Zoom Out Reset image size

Given the large fraction of nonresonant stable configurations, we concluded that the system HD 202696 is most likely dominated by secular interactions at low eccentricities. Figure 9 shows the orbital evolution of three confident stable configurations that have a mean period ratio of $\overline{{P}_{\mathrm{rat}.}}$ ∼ 1.83 exhibiting Δω ∼ 0°, Δω ∼ 180°, and Δω circulating, respectively. The various stable configurations show libration amplitudes in their eccentricities between 0 and 0.15 but with mean values consistent with the posterior eccentricity distribution. The secular timescales of these oscillations strongly depend on the initial orbital geometries and planetary masses, ranging mostly between ∼200 and ∼1000 yr.

Zoom In Zoom Out Reset image size

The HD 202696 system is another example of a growing number of multiplanetary systems around evolved intermediate-mass stars with semimajor axes close to the critical stability boundary. Including HD 202696, there are now seven tightly packed planetary systems with period ratios smaller than or about 2:1.

Earlier examples are the 24 Sex and HD 200964 systems (Johnson et al. 2011). A 1.5 M_⊙ subgiant, 24 Sex is orbited by two Jovian planets with minimum dynamical masses of m_b  ∼ 2 and m_c ∼ 0.9 M_Jup and periods of P_b ≈ 450 and P_c ≈ 883 days. The best-fit period ratio for the 24 Sex system is slightly below 2:1, but the planets are almost certainly involved in a 2:1 MMR (Johnson et al. 2011; Wittenmyer et al. 2012). The HD 200964 system is similar to 24 Sex in terms of stellar mass and minimum planetary masses, but it has a much more compact orbital configuration, with periods of P_b ≈ 610 and P_c ≈  830 days. The HD 200964 system is only stable if the orbits are in a 4:3 MMR (Wittenmyer et al. 2012; Tadeu dos Santos et al. 2015).

The 1.5 M_⊙ subgiant star HD 5319 is another 4:3 MMR candidate host (Robinson et al. 2007; Giguere et al. 2015). With minimum planetary masses of ${m}_{{\rm{b}}}\sin i$ ∼ 1.8 and ${m}_{{\rm{c}}}\sin i$ ∼ 1.2 M_Jup and periods of P_b  ≈ 640 and P_c  ≈ 890 days, this system is dynamically very challenging. Giguere et al. (2015) showed that, although most of the possible configurations of HD 5319 are unstable, stable solutions at the exact 4:3 MMR do exist. Therefore, the 4:3 MMR seems like the most reasonable explanation of the coherent Doppler signal of HD 5319.

A 1.7 M_⊙ slightly evolved star, HD 33844 was reported by Wittenmyer et al. (2016) to host two Jovian-mass planets with minimum dynamical masses m_b ∼ 1.9 and m_c ∼ 1.7 M_Jup and periods of P_b  ≈  550 and P_c ≈ 920 days. These authors concluded that the planetary pair is most likely trapped in the 5:3 MMR, but stable solutions consistent with 8:5 and 12:7 MMR were also found within the 1σ uncertainty in the semimajor axes.

A 1.8 M_⊙ giant, HD 47366 has two nearly equal-mass planetary companions (Sato et al. 2016). With ${m}_{{\rm{b}}}\sin i$ ∼1.8 and m_c sin i ∼ 1.9 M_Jup and periods of P_b ≈ 363 and P_c ≈ 685 days, the system is dynamically challenging. Within the uncertainties of the orbital parameters given, the initial period ratio is close to 15:8, where the system could be stable if the planets are on circular orbits. Sato et al. (2016) also proposed that HD 47366 could be locked in a retrograde 2:1 MMR. Tightly packed retrograde configurations are, however, difficult to reconcile with current planet formation scenarios. Marshall et al. (2019) showed that alternatives to these configurations for HD 47366 exist, including stable prograde solutions in the 2:1 MMR.

The η Ceti system is composed of a giant star orbited by a planet pair with minimum dynamical masses of m_b ∼ 2 and m_c ∼ 3 M_Jup and nominal best-fit periods of P_b ≈ 400 and P_c ≈ 750 days. Trifonov et al. (2014) found two possible stability regimes for η Ceti: either it is a strongly interacting two-planet system locked in an anti-aligned 2:1 MMR or, with lower probability, it could be a low-eccentricity two-planet system dominated by secular interactions close to a high-order MMR like 9:5, 11:6, 13:7, or 15:8, resembling to some extent the HD 202696 system.

Figure 10 shows all known RV multiple-planet systems around intermediate-mass stars with masses larger than 1.3 M_⊙.¹³ The stellar host mass is shown as a function of the logarithm of the system's period ratio. A total of 21 planetary¹⁴ systems are shown in Figure 10. The seven planetary systems found between the 2:1 and 4:3 first-order MMRs (and briefly described above) represent 33% of all known massive multiple-planet systems around intermediate-mass stars. This fraction is rather large given the broad range of the observed planetary period ratios. The planet pair occurrence rate at low period ratios, of course, could be a result of a selection effect, since these systems are likely to be discovered first by precise Doppler surveys. However, a more intriguing possibility is that the increasing number of these tight planetary pairs indicates the existence of a massive multiple-planet system population, some of which have successfully broken the 2:1 MMR during the planet migration phase.

Zoom In Zoom Out Reset image size

It should be noted that a strongly interacting 2:1 MMR system could, in principle, exhibit prominent oscillations of its semimajor axis (and period ratio). The secular oscillation timescales are usually longer than the temporal baseline of the Doppler measurements and therefore normally only marginally detected in the best-fit modeling. An example is the η Ceti system, where best-fit dynamical modeling of the osculating orbits suggests P_rat. ∼ 1.86, but the stable solutions of the system have $\overline{{P}_{\mathrm{rat}.}}$ = 2:1, with an MMR libration amplitude of ${\rm{\Delta }}\overline{{P}_{\mathrm{rat}.}}\sim 0.2$ and timescales of ∼500 yr (Trifonov et al. 2014). In this context, the 24 Sex, HD 47366, and even HD 202696 systems could in fact be true 2:1 MMR systems with oscillating period ratios, observed slightly below the 2:1 period ratio at the present epoch. More precise Doppler data for these systems taken in the future, with an extended temporal baseline, could reveal the true orbital configurations of these systems; we encourage such extended RV campaigns. Still, the HD 33844 system is likely stable at or near the 5:3 MMR, while the HD 5319 and HD 200964 planetary systems can only be stable if they are involved in the 4:3 MMR. The latter systems are dynamically not consistent with the 2:1 MMR and thus indicate that planet pairs could break the 2:1 MMR and settle at other stable orbits in lower-order MMRs.

The HD 202696 system, like the majority of the known Jovian-mass planetary systems, is currently located within the so-called "ice line," where it is unlikely to form gas giants. This can be easily seen from Figure 11, which shows the planetary systems from Figure 10 versus their planetary semimajor axes. The dashed blue line denotes the approximate ice-line radius, beyond which it is commonly accepted that planetary embryos are able to accumulate icy material and grow massive enough to start gas accretion from the protoplanetary disk to eventually become Jovian-mass planets. While accreting disk material, massive planets undergo inward migration toward the stellar host before the disk dissipates, and planets end up at their observed semimajor axes. We note that the ice line plotted in Figure 11 is calculated for the protoplanetary phase assuming only stellar irradiation, where the disk temperature would be T_irr = 150 K (Ida et al. 2016), but adopting a stellar mass–luminosity relation of the form L = M^3.5, which is reasonable for main-sequence stars. We are aware that this is a crude approximation, but it is sufficient to demonstrate the generally "warm orbits" of the discovered Jovian-mass planets orbiting intermediate stars.

Zoom In Zoom Out Reset image size

Figure 11, in connection with current planet formation theories, suggests that the massive planet pairs have formed further out than observed today and have migrated inward together. If the migration is convergent, there is a high probability of planets being trapped in an orbital resonance, with a final stop at the strong 2:1 MMR (e.g., Lee & Peale 2002; Beaugé et al. 2003). Perhaps due to not yet well understood dynamical processes between the disk and the planets, pairs of planets may have broken the 2:1 MMR and migrated further in.

André & Papaloizou (2016) studied systems similar to 24 Sex and η Ceti and were able to reproduce the 2:1 MMR in these systems. However, in their simulations, there were also cases in which the planets passed through the 2:1 MMR, which made them undergo capture into the destabilizing 5:3 resonance. André & Papaloizou (2016) concluded that finding systems in the 5:3 MMR is unlikely.

In this context, the 4:3 MMR systems are also very challenging from the formation perspective. Rein et al. (2012) failed to reproduce the 4:3 MMR of the HD 200964 system using the standard formation scenarios, such as convergent migration, scattering, and in situ formation. On the other hand, Tadeu dos Santos et al. (2015) were able to closely reproduce the 4:3 resonant dynamics of HD 200964 by including interactions between type I and type II migration, planetary growth, and stellar evolution. As pointed out by Tadeu dos Santos et al. (2015), however, the 4:3 MMR capture of the HD 200964 system occurred only for a very narrow range of initial conditions.

In general, massive planets that have reached the gap opening mass migrate via type II migration, which follows the viscous disk evolution (Lin & Papaloizou 1986). In this case, a two-planet pair would migrate on the same timescale, and it is very hard for the outer planet to catch up with the inner planet, let alone to break free from a resonance trapping via planet–disk interactions.

However, if the outer planet is less massive and does not open a full gap, it can undergo type III migration (Masset & Papaloizou 2003). This migration mode is even faster than type I migration and can easily allow outer planets to catch up with inner, more massive planets and be trapped in resonance (Masset & Snellgrove 2001).

This situation has been studied in the framework of the Grand–Tack scenario for Jupiter and Saturn (Walsh et al. 2011; Pierens et al. 2014). In particular, Pierens et al. (2014) found that the inward-migrating Saturn could, under certain disk conditions, have broken free from the 2:1 resonance and migrated closer to Jupiter, where trapping could have occurred, e.g., in the 3:2 resonance.

While the planets are still embedded in the gas disk, they will continue to accrete gas and grow. The outer planet can feed from a larger gas reservoir than the inner one because it has direct access to the disk's gas supply with the full disk accretion rate. The inner planet can only accrete the gas that passed the outer planet, which is typically only up to 20% of the disk's accretion rate (Lubow & D'Angelo 2006). Eventually, the outer planet can outgrow the inner planet if the disk lifetime is long enough. During their final growth, the dynamical interactions between the giant planets change, which could lead to small orbital disturbances breaking the final resonance.

A possible escape from the 2:1 MMR could also be a resonance overstability, as in the context of Goldreich & Schlichting (2014). In this case, convergent planetary migration with strongly damped eccentricities may only lead to a temporal capture at the 2:1 MMR. However, these scenarios require further investigation, which is beyond the scope of this paper.