RETIRED A STARS AND THEIR COMPANIONS. VI. A PAIR OF INTERACTING EXOPLANET PAIRS AROUND THE SUBGIANTS 24 SEXTANIS AND HD 200964^*

We obtained initial-epoch observations of 24 Sex at Lick Observatory in 2005 February, and since then we have obtained 50 RV measurements. After noticing time-correlated RV variations, we began additional monitoring at Keck Observatory in 2008 December, where we have obtained 24 additional RV measurements. The RVs from both observatories are listed in Table 2, along with the Julian Dates (JD) of observation and the internal measurement uncertainties (without jitter). Figure 1 shows the RV time series from both observatories, where the error bars represent the quadrature sum of the internal errors and 5 m s⁻¹ of jitter.

Zoom In Zoom Out Reset image size

Table 4. Orbital Parameters from Non-interacting Two-planet Solution

Parameter	24 Sex ba	24 Sex c	HD 200964 b	HD 200964 b
Period (days)	455.2(3.2)	910(21)	630.6(9.3)	829(21)
Tpb (JD)	2454758(30)	2454941(30)	2454916(30)	2455029(130)
Eccentricity	0.184(0.029)	0.412(0.064)	0.111(0.030)	0.113(0.076)
K (m s−1)	33.2(1.6)	23.5(2.9)	35.2(2.7)	22.1(2.3)
ω (deg)	227(20)	172.0(9)	223(20)	301(50)
MPsin i(MJup)	1.6(0.2)	1.4(0.2)	1.9(0.2)	1.3(0.2)
a (AU)	1.41(0.03)	2.22(0.06)	1.71(0.04)	2.03(0.06)
Lick rms (m s−1)	7.7	⋅⋅⋅	7.6	⋅⋅⋅
Keck rms (m s−1)	4.8	⋅⋅⋅	5.3	⋅⋅⋅
Jitter (m s−1)	5.0	⋅⋅⋅	5.0	⋅⋅⋅
$\sqrt{\chi ^2_\nu }$	1.14	⋅⋅⋅	1.15	⋅⋅⋅
Nobs Lick	50	⋅⋅⋅	61	⋅⋅⋅
Nobs Keck	24	⋅⋅⋅	35	⋅⋅⋅

Notes. ^aHD 90043. ^bTime of periastron passage.

Download table as:  ASCIITypeset image

Bowler et al. (2010) reported evidence of a two-planet system around 24 Sex, but the data at the time could not provide a unique orbital solution. An additional season of observations has provided stronger constraints on the possible orbits of the two planets. We fitted a model consisting of two non-interacting planets and a 1.54 M_☉ star orbiting their mutual center of mass. We find that two Keplerians provide an acceptable fit to the data with an rms scatter of 6.8 m s⁻¹ and a reduced $\sqrt{\chi ^2_\nu }$ =1.14. The inner planet has a period of P = 455.2 ± 3.2 days, velocity semiamplitude K = 33.2 ± 1.6 m s⁻¹, and eccentricity e = 0.184 ± 0.029. The outer planet has P = 910 ± 21 days, K = 23.5 ± 2.9 m s⁻¹, and e = 0.412 ± 0.064. Together with our stellar mass estimate these spectroscopic orbital parameters yield semimajor axes {a_b, a_c} = {1.41, 2.22} AU and minimum planet masses {M_b, M_c}sin i = {1.6, 1.4} M_Jup. A more detailed, dynamical analysis presented in Section 6 revises this two-Keplerian solution under the constraint of long-term stability.

We began monitoring HD 200964 at Lick Observatory in 2007 July and have obtained 61 RV measurements. Time-correlated variations in the star's RVs prompted additional monitoring at Keck Observatory where we have obtained 35 measurements since 2007 October. The RVs from both observatories are listed in Table 3, along with the JD of observation and the internal measurement uncertainties (without jitter). Figure 2 shows the RV time series from both observatories, where the error bars represent the quadrature sum of the internal errors and 5 m s⁻¹ of jitter.

Zoom In Zoom Out Reset image size

As is the case for 24 Sex, Bowler et al. (2010) reported evidence of a two-planet system around HD 200964, and an additional season of observations has provided stronger constraints on the possible orbits of the planets in the system. We find that a two-Keplerian model provides an acceptable fit to the data with an rms scatter of 6.8 m s⁻¹ and a reduced $\sqrt{\chi ^2_\nu }$ =1.14. The inner planet has a period of P = 630.6 ± 9.3 days, velocity semiamplitude K = 35.2 ± 2.7 m s⁻¹, and eccentricity e = 0.111 ± 0.030. The outer planet has P = 829 ± 21 days, K = 22.1 ± 2.3 m s⁻¹, and e = 0.113 ± 0.076. Together with our stellar mass estimate these spectroscopic orbital parameters yield semimajor axes {a_b, a_c} = {1.71, 2.03} AU and minimum planet masses {M_b, M_c}sin i = {1.9, 1.3} M_Jup. A more in-depth dynamical analysis presented in Section 6 revises this two-Keplerian solution under the constraint of long-term stability.

We analyze the RV observations using a Bayesian framework following Ford (2005, 2006a). We assume priors that are uniform in logarithmic intervals of orbital period, eccentricity, argument of pericenter, mean anomaly at epoch, and the velocity zero point. For the velocity amplitude (K) and stellar jitter (σ_j), we adopt a prior of the form p(x) = (x + x_o)⁻¹[log(1 + x/x_o)]⁻¹, with K_o = σ_j,o = 1 m s⁻¹. For a discussion of priors, see Ford & Gregory (2007). The likelihood for RV terms assumes that each RV observation (v_i) is independent and normally distributed about the true RV with a variance of σ²_i + σ²_j, where σ_i is the published measurement uncertainty and σ_j is a jitter parameter that accounts for additional scatter due to stellar variability, instrumental errors, and/or inaccuracies in the model (i.e., neglecting planet–planet interactions or additional, low-amplitude planet signals). It is important to bear in mind that by adding the variance due to jitter to the measurement errors we are making the simplifying assumption that jitter is normally distributed about our model, which is certainly not entirely true as many sources of jitter are systematic rather than purely random.

In our initial phase of analysis, we use an MCMC method based upon Keplerian orbit fitting to calculate a sample from the posterior distribution (Ford 2006a). We calculate multiple Markov chains, each with ∼2 × 10⁸ states. We discard the first half of the chains and calculate Gelman–Rubin test statistics for each model parameter and several ancillary variables. We find no indications of non-convergence. Thus, we randomly choose a subsample (∼ 16,000 samples) from the posterior distribution for further investigation.

Next, we use this subsample as the basis for a much more computationally demanding analysis that uses fully self-consistent n-body integrations to account for planet–planet interactions when modeling the RV observations. We again perform a Bayesian analysis, but replace the standard MCMC algorithm with a Differential Evolution Markov Chain Monte Carlo (DEMCMC) algorithm (ter Braak 2006; Veras & Ford 2009, 2010). In the DEMCMC algorithm, each state of the Markov chain is an ensemble of orbital solutions. The candidate transition probability function is based on the orbital parameters in the current ensemble, allowing the DEMCMC algorithm to sample more efficiently from high-dimensional parameter spaces that have strong correlations between model parameters. More details of this DEMCMC algorithm and associated tests of its accuracy will be presented in a forthcoming paper (R. P. Nelson et al. 2011, in preparation)

The priors for model parameters are the same as those of the Keplerian MCMC simulations. The initial conditions of the n-body simulations are calculated by converting between Keplerian and Cartesian coordinates. In this paper, we consider only coplanar, edge-on two-planet systems, and hence the planetary masses resulting from this n-body analysis (see Table 6) will be minimum masses. We note that relaxing the assumption of coplanarity would (1) be significantly more computationally intensive, and (2) potentially lead to a much greater proportion of unstable systems being generated, as the planetary masses used in the integrations would be increased. Previous experience suggests that inclinations greater than >45° frequently lead to highly unstable results, while inclinations <30° can often be rather similar to uninclined systems. As such, we emphasize that our results should not be interpreted as giving the precise planetary solutions (masses, period ratios, etc.), but rather as being indicative that stable coplanar solutions exist, and that further detailed investigation is warranted.

For the n-body integrations, we use a time-symmetric fourth-order Hermite integrator that has been optimized for planetary systems (Kokubo et al. 1998). We extract the RV of the star (in the barycentric frame) at each of the observation times for comparison to RV data. During the DEMCMC analysis, we also impose the constraint of short-term (100 years) orbital stability. We check whether the planetary semimajor axes remain within a factor of 50% of their starting value, and that no close-approaches occur within 0.1× the semimajor axis during the 100 year n-body integration. Any systems failing these tests are rejected as unstable (regardless of the quality of the fit to RV data). Thus, the DEMCMC simulations avoid orbital solutions that are violently unstable. In our DEMCMC simulations, this process is repeated for 10,000 generations, each of which contains 16,000 systems, for a total of ∼10⁸ n-body integrations in each DEMCMC simulation.

Since the DEMCMC simulations only require stability for ∼100 years, the orbital solutions in the final generation may or may not be stable for longer timescales. Since nearly all of these systems are strongly interacting, we take this final generation (16,000 systems) and demand that they also be stable (according to the same criteria above) over the course of a 10⁷ year integration, performed using the hybrid Bulirsch-Stoer/Symplectic integrator Mercury (Chambers 1999). Only the orbital solutions which are stable over the course of this long-term n-body integration are regarded as being acceptable solutions. More details on the dynamical analysis performed and the results obtained will be presented in a companion paper (M. Payne et al. 2011, in preparation).

We find that the n-body DEMCMC routine results for 24 Sex are concentrated around solutions with {P_b, P_c} ≈ {450, 900} days, i.e., close to the 2:1 resonance, and that they span a significantly smaller range of parameter space than do the double-Keplerian fits. The n-body RV fitting aspect of the DEMCMC routine acts to shrink the parameter space, whereas the stability requirement in the routine has the effect of shifting the DEMCMC solutions toward the 2:1 period ratio mark, primarily by favoring lower periods for the inner planet. We illustrate in Figure 4 the difference in the endpoints of the two analyses, showing a scatter plot of the planetary periods at the end of both the Keplerian analysis and the n-body analysis.

Zoom In Zoom Out Reset image size

We have applied the DEMCMC method described in Section 6.1 multiple times using different values for the initial ensemble of orbital solutions. Each time, these simulations reached qualitatively similar results to those of Figure 4, but quantitatively there are signs that the method has not fully converged. Therefore, we do not interpret the results as a precise estimate of the posterior probability distribution. Instead, we use these results to demonstrate that there exist orbital solutions that are both stable and consistent with the Doppler observations. Most importantly, all our simulations indicate that the posterior distribution is concentrated at solutions with a ratio of orbital periods between ∼1.9 and ∼2.1 and very close to the 2:1 MMR. Thus, we conclude that the current observations strongly favor orbital solutions in (or at least near) the 2:1 MMR.

Next, we performed long-term stability tests for the 16,000 orbital solutions in the final generation of the n-body DEMCMC analysis. We find that ∼90% of the orbital solutions are unstable over the course of 10⁷ years of integration. We consider the remaining ∼1000 (of the 16,000 orbital solutions) that are stable for 10⁷ years to be plausible orbital solutions given both the RV data for 24 Sex and the requirement of long-term dynamical stability.

We show in Figure 4 that the stable systems remain strongly clustered around the 2:1 period commensurability region. Taking all of these long-term-stable systems into account, we find that the data indicate that the inner planet has a period P_b= 452.8^+2.1_−4.5 days, semimajor axis a_b= 1.333^+0.004_−0.009 AU, eccentricity e_b= 0.09^+0.14_−0.06, and mass M_bsin i= 1.99^+0.26_−0.38 M_Jup, while the outer planet has a period P_c= 883.0⁺³²₋₁₄ days, semimajor axis a_c= 2.08^+0.05_−0.02 AU, eccentricity e_c= 0.29^+0.16_−0.09, and mass M_csin i= 0.86^+0.35_−0.22 M_Jup. We detail all the fitted parameters from this analysis in Table 6.

When we examine specific systems in detail, we find that the pericenter of the outer planet overlaps the location of the pericenter of the inner planet, meaning that the planets therefore experience detectable gravitational interactions. Individual systems typically undergo short-term eccentricity and semimajor axis oscillations with amplitudes that can be at least as large as the spread in median values quoted in Table 6. As an example, one stable system was found to exhibit planetary eccentricity oscillations with a period ${\sim} 800\, \rm years$ and amplitude 0.02 < e_b < 0.20 and 0.15 < e_c < 0.36, while the semimajor axes oscillated with a period ∼25 years and amplitudes 1.330 < a_b < 1.338 AU and 2.05 < a_c < 2.12 AU.

We note that the PDFs from the DEMCMC analysis are bimodal. To investigate the cause of the bimodality, we performed a similar n-body+DEMCMC analysis without the requirement of short-term stability. The resulting PDFs have a single broader mode consistent with the output of the Keplerian MCMC analysis. We conclude that the bimodal nature of the PDF from the dynamical analysis is most likely the result of demanding short-term stability.

The overlap of the pericenters suggests that an MMR is needed to stabilize the system over long timescales by preventing close encounters. However, an analysis of the resonant angles θ_2,1 = 2λ_c − 1λ_b − 1ϖ_2,1 (where λ is the mean longitude and ϖ is the longitude of pericenter) suggests that, for some of the long-term-stable systems, the planets circulate, rather than librate, i.e., the systems are observed to have angular ranges for θ_2,1 ∼ 360°. We therefore caution that the current state of the observations and dynamical analysis cannot confirm whether the system is truly in a resonant state, and hence further work will be required. A more detailed investigation (M. Payne et al. 2011, in preparation) will probe the nature of these dynamical interactions in 24 Sex and HD200964 in greater detail.

Finally, we note that, as discussed in Section 6.1, jitter is allowed to vary during the DEMCMC procedure. As such we find a best-fit value for the jitter in the same manner as we do for the various planetary orbital parameters. From the 24 Sex analysis, we find a large jitter value of 9.9^+2.9_−1.2 m s⁻¹. This is substantially larger than the value of 5 m s⁻¹ assumed in Section 4, indicating that 24 Sex has intrinsic RV variability at the high end of the observed distribution, or that there are other unmodeled sources of variability such as an additional planet in the system. However, a periodogram of the residuals shows no significant power above the noise. Using the approach of Bowler et al. (2010), we can rule out additional companions with masses greater than 0.3 M_Jup out to 1 AU and "hot" planets with masses M_Psin i <0.1 M_Jup within 0.1 AU. Additional monitoring at Keck and a more in-depth dynamical analysis will help clarify the situation.

The results of our n-body DEMCMC simulations for HD 200964 are again more tightly confined than the results of the double-Keplerian analysis. As an example of this we illustrate in Figure 5 the difference in the endpoints of the two analyses, showing a scatter plot of the planetary periods at the end of both the Keplerian analysis and the n-body analysis, as well as histograms for the same data. Again, we note that the shift to an n-body analysis acts to concentrate the preferred region into a smaller area, while the stability requirements act to shift the favored parameters closer toward the 4:3 period commensurability.

Zoom In Zoom Out Reset image size

As for the case of 24 Sex, we have applied our DEMCMC method multiple times using numerous randomized initial conditions. Again, we find qualitatively similar results from all sets of the DEMCMC analysis, but with some indication that the results have not truly converged. As such, we interpret these results as demonstrative that stable systems exist which can fit the RV data, and moreover, the extremely narrow range of period ratios favored by the analysis (∼1.32 to ∼1.36) shows that the current observations strongly favor orbital solutions in or close to the 4:3 MMR.

Next, we test the long-term orbital stability of the orbital solutions identified by the DEMCMC analysis. For HD 200964, >90% of the systems are clearly unstable during a 10⁷ year integration (i.e., they experience a collision or a change in semimajor axes of more than 50%.) We find that ∼1000 of the 16,000 simulations remain stable for 10⁷ years. We interpret these as plausible orbital solutions consistent with the RV observations of HD 200964.

As illustrated in Figure 5, the stable systems occupy a region of parameter space corresponding to a region of parameter space near the 4:3 period ratio. Taking only these long-term-stable systems into account, we find that the inner planet has a period P_b= 613.8^+1.3_−1.4 days, semimajor axis a_b= 1.601^0.002_−0.002 AU, eccentricity e_b= 0.04^+0.04_−0.02, and mass M_bsin i= 1.85^+0.14_−0.08 M_Jup, while the outer planet has a period P_c= 825.0^+5.1_−3.1 days, semimajor axis a_c= 1.95^+0.008_−0.005 AU, eccentricity e_c= 0.181^+0.024_−0.017, and mass M_csin i= 0.90^+0.12_−0.06 M_Jup. We detail all the fitted parameters from this analysis in Table 6.

We find that most of the stable planetary orbits overlap, producing strongly interacting systems, resulting in significant oscillations in the semimajor axes and the eccentricities of both planets. An example stable solution exhibits observable eccentricity oscillations (0.03 < e_b < 0.1 and 0.13 < e_c < 0.18) on an approximately 250 year timescale and semimajor axis oscillations (1.57 < a_b < 1.58 AU and 1.9 < a_c < 1.93 AU) on an approximately 60 year timescale. As in the 24 Sex system, some of the long-term-stable systems examined here for HD 200964 seem to be circulating rather librating and we are conducting a more detailed examination to ascertain the fraction of system solutions which actually librate. We find a jitter value of 8.23^+0.38_−0.88 m s⁻¹, which is larger than the empirical jitter estimate of Johnson et al. (2010).