Characterizing the Orbital and Dynamical State of the HD 82943 Planetary System With Keck Radial Velocity Data1footnote 11footnote 1 Based on observations obtained at the W. M. Keck Observatory, which is operated jointly by the University of California and the California Institute of Technology.

Characterizing the Orbital and Dynamical State of the HD 82943 Planetary System With Keck Radial Velocity Data1


We present an updated analysis of radial velocity data of the HD 82943 planetary system based on 10 years of measurements obtained with the Keck telescope. Previous studies have shown that the HD 82943 system has two planets that are likely in 2:1 mean-motion resonance (MMR), with the orbital periods about 220 and 440 days (Lee et al. 2006). However, alternative fits that are qualitatively different have also been suggested, with two planets in a 1:1 resonance (Goździewski & Konacki 2006) or three planets in a Laplace 4:2:1 resonance (Beaugé et al. 2008). Here we use minimization combined with parameter grid search to investigate the orbital parameters and dynamical states of the qualitatively different types of fits, and we compare the results to those obtained with the differential evolution Markov chain Monte Carlo method. Our results support the coplanar 2:1 MMR configuration for the HD 82943 system, and show no evidence for either the 1:1 or 3-planet Laplace resonance fits. The inclination of the system with respect to the sky plane is well constrained at degrees, and the system contains two planets with masses of about 4.78 and 4.80 (where is the mass of Jupiter) and orbital periods of about 219 and 442 days for the inner and outer planet, respectively. The best fit is dynamically stable with both eccentricity-type resonant angles and librating around 0.

celestial mechanics – planetary systems – stars: individual (HD 82943)

1 Introduction

To date, more than 125 extrasolar multiple-planet systems have been confirmed, in which there are nearly 40 systems that are suspected to contain planets in or near mean-motion resonance (MMR) (see, e.g., Wright et al. 2011a, Wright et al. 2011b14), hinting that MMRs play an important role in the orbital configurations of exoplanetary systems. Candidates from the Kepler transit survey show a significant fraction of adjacent planets with their period ratio in or near first order mean-motion resonance commensurabilities (Lissauer et al. 2011; Fabrycky et al. 2012). However, Veras & Ford (2012) have analytically shown that the vast majority of Kepler systems with two near-resonance transiting candidates cannot be in resonances. As such, the fact that the orbital period ratio is close to a low-order integer ratio does not necessarily indicate that a system is in resonance, and that a more detailed examination of the system dynamics is required to rule in or out this possibility.

The first pair of exoplanets discovered to be in mean-motion resonance was the GJ 876 planetary system (Marcy et al. 2001). By fitting radial velocity (RV) measurements, it is well established that the pair of planets is in a deep 2:1 MMR with both lowest order, eccentricity-type MMR angles,


and the secular apsidal resonance angle,


librating about with small libration amplitudes (Laughlin & Chambers 2001; Rivera & Lissauer 2001; Laughlin et al. 2005). Here, is the mean longitude, is the longitude of periapse, and subscripts 1 and 2 represent the inner and outer planets, respectively. An updated RV analysis by Rivera et al. (2010) showed that the GJ 876 planetary system contains an additional outer planet, and it is in a Laplace resonance with the previously known 2:1 resonant pair. The orbital configuration of the resonant pair of planets b and c in GJ 876 system differs from that of the Galilean satellites Io and Europa, which is also in 2:1 MMR but with librating around and both and librating around . The anti-aligned corotational configuration ensures that Io and Europa remain stable in the 2:1 MMR because of their small eccentricities (Lee & Peale 2002), whereas moderate and large eccentricities allow a wide variety of other stable 2:1 MMR configurations, including the GJ 876 configuration (Lee & Peale 2002; Beaugé et al. 2003; Ferraz-Mello et al. 2003; Lee 2004; Beaugé et al. 2006). The 2:1 MMR configuration of the GJ 876 system can be established by convergent orbital migration caused by disk-planet interactions, and thus provided constraints on the migration processes within, and the physical environment constituting, the protoplanetary disk. -body simulations with forced inward orbital migration of the outer planet for the GJ 876 system show that either significant eccentricity damping during migration or quick dispersal of the disk after resonance capture is necessary to explain the observed eccentricity values (Lee & Peale 2002). Although earlier hydrodynamic simulations did not show significant eccentricity damping from disk-planet interactions (Papaloizou 2003; Kley et al. 2004, 2005; Moorhead & Adams 2008), a recent study for the GJ 876 system by Crida et al. (2008) shows that a stronger eccentricity damping from disk-planet interactions can be produced if the inner disk and its interactions with the inner planet are considered.

It is important to increase the sample of well-determined resonant systems, as the GJ 876 system has shown interesting properties in terms of orbital dynamics as well as providing constraints on the evolution of planetary system. This paper focuses on characterizing the orbital and dynamical state of the suspected 2:1 MMR system HD 82943. The detection of the HD 82943 planetary system was announced in European Southern Observatory (ESO) press releases15, with the first planet discovered in 2000 and the second planet in 2001. Based on the orbital parameters posted in the Geneva group web page, Ji et al. (2003) claimed that the fit tends to be unstable unless the system is locked in 2:1 MMR. Mayor et al. (2004) published a double-Keplerian fit of the HD 82943 system based on 142 CORALIE RV measurements, and claimed that this system is in 2:1 MMR with orbital periods of about 220 and 440 days, and planetary mass ratio of . However, the actual data set was never published.

By direct -body integrations, Ferraz-Mello et al. (2005) found that the solution of the HD 82943 system found by Mayor et al. (2004) was unstable on the order of yr. The instability of this solution did not change when the integration was started at different initial times instead of starting only at the first observational time epoch of Mayor et al. (2004), but with the other Keplerian orbital elements the same. Using the CORALIE RV data extracted from the graphs in Mayor et al. (2004), Ferraz-Mello et al. (2005) found a best fit with rms = based on a double-Keplerian model, and explored the parameter space around the best fit. The best fit was unstable, but the rms as a function of the primary parameters was shallow around the minimum and there were many stable good fits giving slightly higher rms. In a statistical sense the real solution may correspond to one of the many stable good fits.

With 23 additional new RV data from the Keck Observatory, Lee et al. (2006) analyzed the combination of the CORALIE data set (which was extracted from graphs in Mayor et al. 2004) and the Keck data set, using the double-Keplerian model combined with parameter grid search. The best fit had and rms = , but it was also dynamically unstable. Fits in the parameter grid as a function of the eccentricity and argument of periapse of the outer planet were systematically explored. The minimum was shallow around the best fit, with good fits that have just slightly higher being stable. Dynamical stability exploration of fits in the parameter grid showed that all fits that are stable are in 2:1 MMR, assuming that the planets are in coplanar orbits. The results suggested that the HD 82943 system was almost certainly in 2:1 MMR.

Using the same data sets, Goździewski & Konacki (2006) undertook a dynamical fitting analysis using a hybrid algorithm which introduces an orbital instability penalty factor into the minimization. The edge-on coplanar 2:1 MMR best fit with approximately the same as the best fit found by Lee et al. (2006) was unstable, but two islands of stable 2:1 MMR fits near the best fit in parameter space were found. Using a genetic search algorithm, qualitatively different fits associated with the 1:1 eccentric resonance configuration were found. The best stable fit with a 1:1 resonance is highly mutually inclined, and just slightly worse than the best fit of 2:1 MMR in terms of . This raises the interesting possibility that similar radial velocities could be produced by 2:1 MMR and 1:1 eccentric resonance orbits for the HD 82943 system.

Based on the same data sets, Beaugé et al. (2008) analyzed the sensitivity of the fits with respect to data set by the so-called Jackknife method. The parameters of the best fit are sensitive to a few data points in the Keck data set, suggesting a plausible inconsistency between the CORALIE and Keck data set. A new type of fit corresponding to a coplanar 3-planet orbital configuration which included an additional outer planet was suggested. Beaugé et al. (2008) claimed that the quality of the edge-on coplanar 3-planet dynamical best fit was a significant improvement accounting for the increased number of fitting parameters. The best-fit 3-planet configuration was unstable, but interestingly, a fit that corresponds to a stable Laplace resonance configuration with was suggested.

Two major problems remain in the orbital characterization of the HD 82943 planetary system. First, the minimum of and rms of the 2:1 MMR best fit was very shallow in parameter space and almost all of the best fits found to date have been dynamically unstable. Second, based on the combined RV data of CORALIE and Keck observations, three qualitatively different type of fits, i.e., 2:1 MMR fits, 3-planet fits and 1:1 resonance fits, with almost equally good quality have been suggested, complicating the situation on this system. The purpose of this analysis of the HD 82943 system is to see whether we can discriminate between these three different type of fits, and whether our -body dynamical fit can constrain the inclinations and masses of the planets with the longer time span of the Keck data. In our analysis, we only use the Keck data instead of combining with the CORALIE data for fitting for the following reasons: the CORALIE data were never published except in graphical form, and the analysis by Beaugé et al. (2008) implied an inconsistency between the CORALIE and Keck RV measurements. Although the Keck data is more accurate and has a longer time span, the large number of CORALIE data points means that they could strongly influence the properties of the fits if they are included.

This paper is organized as follows: In §2, we introduce the stellar properties of HD 82943 and the available RV data from observations at Keck Observatory. In §3, we introduce our -body dynamical fitting method based the framework of minimization. In §4, we present our fitting results which are organized for discriminating three qualitatively different type of fits described above. We show that the coplanar 2:1 MMR configuration is preferred, and the inclination of the system with respect to the sky plane is well constrained at about . In §5, distributions and uncertainties of the best-fit parameters are presented. In §6, we compare results from minimization combined with grid search to those from the Markov chain Monte Carlo analysis. Finally, our conclusions are summarized in §7.

2 Stellar Characteristics and Radial Velocity Measurements

The stellar properties of the G0 star HD 82943 were summarized by Mayor et al. (2004), Fischer & Valenti (2005) and Takeda et al. (2007). A stellar mass of for HD 82943 was suggested by Laws et al. (2003) and Fischer & Valenti (2005) independently, and a smaller stellar mass of was determined by Santos et al. (2000). A median mass of was adopted by Mayor et al. (2004), Lee et al. (2006) and Goździewski & Konacki (2006). Here we adopt the slightly lower stellar mass of from the more recent spectroscopic analysis by Takeda et al. (2007). The age of the star HD 82943 is about 3 Gyr (Mayor et al. 2004; Takeda et al. 2007). Thus any good fit that interprets the configuration of the HD 82943 planetary system should be constrained to be dynamically stable on timescales of the order of yr. Uncertainties in the stellar radial velocity may arise from the stellar jitter caused by acoustic -modes, turbulent convection, star spots, and flows controlled by surface magnetic fields. We adopt a stellar jitter of estimated by Lee et al. (2006) for HD 82943, and the jitter is quadratically added to the internal uncertainties of the RV data sets for results in §4 and §5. The jitter is treated as an unknown parameter in the Bayesian analysis of §6.

We began RV measurements for HD 82943 in 2001 with the HIRES echelle spectrograph on the Keck \@slowromancapi@ telescope. We have to date obtained 64 RV measurements spanning about 10 yr, and the measurements are listed in Table 1. This is a significant improvement on the 23 Keck measurements spanning 3.8 yr used by Lee et al. (2006). For the 23 data points previously published by Lee et al. (2006), the RV in Table 1 are slightly different and the uncertainties are smaller, because we have analyzed all spectra uniformly, using an improved Doppler analysis code and a new template taken in 2009 with higher resolution and higher signal-to-noise ratio. The split in the data in Table 1 corresponds to the 2004 upgrade of the Keck-HIRES CCD System. The median of the internal velocity uncertainties for the two sets of Keck data are about 1.5 and 1.2 , respectively. Because the Keck velocities were measured with two different instrumental configurations, we model the two different Keck data sets as having two unknown velocity offsets.

3 Methodology

More than 8 orbital periods of the outer planet of the HD 82943 planetary system have been observed, and the mutual gravitational interactions between the planets are expected to significantly influence their orbital configuration. Therefore, it is important to adopt a self-consistent dynamical analysis for the RV data.

The Levenberg-Marquardt (LM) method (Press et al. 1992) can efficiently search for the local minimum given an initial guess of a set of fitting parameters. A key ingredient in the LM method is the partial derivative of the radial velocity function with respect to each of the fitting parameters , where is the radial velocity function and is the set of fitting parameters. A method for calculating these derivatives in dynamical fitting has been developed from the variational equations that describe the difference of two adjacent orbital motions and its evolution with time. A detailed description of this method and its application to dynamical fitting for multiple planetary systems is given by Pál (2010). Compared to the numerical derivatives method which has been widely used in dynamical RV fitting, the advantage of the variational equation method is the better accuracy in calculating the derivatives and the lower computational cost when the number of planets is not large.

For coplanar orbits, the fitting parameters are the zero-point offset velocity for each RV data set, the inclination of the system’s orbital plane relative to the sky plane and the initial osculating Keplerian orbital elements for each planet, i.e., (, , , , ). For mutually inclined orbits, the fitting parameters are the zero-point offset velocity for each RV data set and the initial osculating Keplerian orbital elements for each planet, i.e., (, , , , , , ). Here, is the amplitude of the radial velocity, is the orbital period, is the orbital eccentricity, is the argument of periapse, is the mean anomaly at the first observational epoch (which is BJD 2452006.913 for the current data set), is the orbital inclination relative to the sky plane and is the ascending node. Since RV fitting can only determine the difference in the ascending nodes of orbits, not the absolute ascending node for each orbit measured from arbitrary reference direction in the sky plane, it is convenient to fix the initial ascending node of the first planet to 0, and the ascending nodes of the other planets are simply treated as the mutual ascending nodes relative to the first planet at the first observational epoch. Note that during the -body integration in our dynamical fitting, all the orbital parameters are allowed to evolve, so the values of all (including ) evolve with time. In some situations variants of the fitting parameters are preferred. For example, when the eccentricities are small, we prefer (, , ) rather than (, , ). The observed radial velocity of the central body is


where is the total mass of the system, is planetary mass, are velocity components for planet in astrocentric coordinates, is the number of planets, is the number of data sets, the unit vector is directed towards the observer and here we set it to be .

The equations of planetary motion in astrocentric coordinates are


Here, is the mass of the central body, is the gravitational constant, for planet , is the position vector and is the velocity vector. Let us denote the left hand sides of Eqs. (5) – (7) as variables , and the right hand sides as functions . The variations denote arbitrary deviations from the original variables , and the evolution of these deviations with time can be calculated by the linearized variational equations


where is the number of parameters. Because of the linear property, the relation between initial variations and variations at arbitrary time can be written in the form


The matrix measures the ratio of the deviations at an arbitrary time, , to the initial deviations, and the initial matrix (at ) is a unit matrix, i.e., . Combining Eqs. (8) and (9), the variational equations can be written as a set of linear differential equations in matrix form


Explicit expressions of the derivatives in Eq. (10) can be found in Baluev (2011). The radial velocity in a dynamical model is the stellar velocity caused by the planets. The partial derivatives of the radial velocity at time with respect to initial orbital parameters are in fact related to gradients of orbital motions with respect to initial orbital parameters. These gradients are calculated from the matrix in the variational equations. Thus, the partial derivatives of the radial velocity can be calculated using the chain rule


Explicit expression for the gradients of for coplanar orbits with variables can be found in Pál (2010).

The initial orbital parameters in hierarchical multiple-planet systems should be interpreted as Jacobi coordinates. The Jacobi coordinates better emulate the assumptions used to fit the parameters to the data than either barycentric coordinates or astrocentric coordinates (Lissauer & Rivera 2001; Lee & Peale 2003). The fitting parameters in Jacobi coordinates are denoted as , and their definition can be found in Lee et al. (2006). In this case, Eq. (11) should be replaced by


where are the variables in Jacobi coordinates at the first observed epoch. The transformation from variables to can be found in Saha & Tremaine (1994), and it can be written in matrix form


Thus, the component in Eq. (12) is substituted with matrix .

The variational equations (Eq. [10]) are integrated simultaneously with the equations of motions (Eqs. [5]-[7]) using the Bulirsch-Stoer integrator, and the RV values in the model and their partial derivatives with respect to the fitting parameters (Eqs. [11] and [12]) are obtained in every observational time. There are several versions of the Bulirsch-Stoer integrator, and we adopt the one implemented in the SWIFT package (Levison & Duncan 1994). Because the LM method is a local minimization method, we use the systematic grid search techniques (Lee et al. 2006) to ensure the search for global best fit.

4 Fitting Results

The HD 82943 planetary system has been studied for years, and previous RV fitting results provide good initial guesses for minimization in our study. However, we must keep in mind that because our RV data sets are different from previous ones (see the discussion in section 1), other solutions may be found. We first fit with Keplerian model, and the best-fit Keplerian model is taken as an initial guess for dynamical fitting.

4.1 2:1 MMR Fits

Coplanar Edge-on Fits

We first conduct a brief Keplerian analysis with the initial parameters adopted from the global best fit (Fit \@slowromancapi@) found by Lee et al. (2006). The reaches a minimum, similar to that of Fit \@slowromancapi@ in the upper left corner of the grid shown in Fig. 3 of Lee et al. (2006). We then explore fits in the grid starting from this minimum. There is a second minimum in the lower left corner of the grid lower in . Based on the lower minimum, we allow all parameters to vary and obtain the best fit with and . The parameters of the best fit is listed in the first part of Table 2. The and rms are significantly improved compared to Fit \@slowromancapi@ of Lee et al. (2006) (), and the Keplerian best fit is dynamically stable for at least yr according to a direct -body integration.

The Keplerian best fit is used as an initial guess for coplanar edge-on dynamical fitting. The LM algorithm converges to a minimum with and . Starting from different initial guesses, we do not find any other local minimum, so this is likely the global best fit for the coplanar edge-on dynamical case. The parameters of the best fit are listed in the second part of Table 2 and the RV curve and residuals are shown in the left panel of Fig. 1. The dynamical fitting is an improvement over the Keplerian fitting, as the of the coplanar edge-on dynamical best fit is less than the of the Keplerian best fit by about 8.8. The eccentricity of the outer planet () is smaller than the best fit () found by Lee et al. (2006). The edge-on best fit is dynamically stable for at least yr, and it is in a 2:1 MMR with both resonance angles and librating around . The evolution of semimajor axes and , eccentricities and and resonance angles and are shown in the right panel of Fig. 1. Notice that although is relatively small in the first observed epoch, it fluctuates because of the large libration amplitudes of the resonance angles. The maximum of is about 0.23 and the minimum is about 0.07.

Parameter grid search is necessary to examine whether there are any other minima. The errors given by the covariance matrix in the LM algorithm suggest that , and are the least certain parameters. The mean anomaly is a time-dependent parameter and it does not directly correlate to the orbital spatial configuration, so it is normally not used for grid searching. We search for local best fits with and fixed in the grid, because when becomes sufficiently small, and are more natural for representing the orbital configuration. After some experiments we decided the range of and to be (0.1 - 0.3) and (0.15 - 0.15) respectively. The left panel of Fig. 2 shows the contours (1.165 - 1.80) for the edge-on fits in grid. The best fit is the only minimum in , with = 1.16 for 10 adjustable parameters. The smooth contours show that the in the parameter space smoothly converge to the minimum. We also search fits in other grids. The grid is the most relevant grid to MMR configuration, and it is shown in the right panel of Fig. 2. Similar to the grid, there is only one minimum in the parameter space.

We use the symplectic integrator SyMBA (Duncan et al. 1998) to perform a direct -body integration with a maximum time of 50,000 yr for each local best-fit model in the grids. The initial semimajor axes of the inner and outer planet are about 0.75 AU and 1.18 AU, respectively. For the system to be considered stable, two criteria have been set: (1) the maximum distance of a planet to the star must be less than 3 AU (5 AU for 3-planet models in §4.2) and the minimum distance must be larger than 0.075 AU; (2) the distance between the planets must be larger than the sum of their physical size assuming their mean density of 1 g . The minimum and maximum distances from the star adopted in criterion (1) are sufficiently different from the initial semimajor axes of the planets that the system is clearly unstable if these limits are exceeded. The dynamical properties of the fits in the grid are shown in Fig. 3. In the left panel, the thick dashed lines are contours of survival time of the integrations, with the thickest line corresponding to 50,000 yr (the maximum integration time) and the thinnest line corresponding to 2,000 yrs. The blank region is the stable region in which all fits are stable. All stable fits are in 2:1 MMR with librating around . Most stable fits have librating around but some fits close to the dynamical stability boundary have circulating. The thin solid black and thin dashed gray (magenta in the color version) curves represent the contours of libration amplitudes for and , respectively. In the right panel, the thick dashed lines are the stability boundary, and the thin curves are contours corresponding to the , and contours of confidence levels based on of 2.3, 6.17, 11.8 (or of 0.043, 0.114, 0.219) larger than the minimum for 2 parameters (Press et al. 1992). The star point represents the fit with the smallest libration amplitudes of both resonant angles (), and it is outside the confidence region. The best-fit model is far away from the stability boundary, with libration amplitudes about and for and , respectively (see also Fig. 1). A large fraction of the confidence regions are in the stable region, and a small fraction of fits in confidence region are unstable. We do similar analysis in the grid, and the results are shown in Fig. 4, with the same contents as the grid. The confidence regions are narrowed by the stability boundaries compared to grid, but almost all fits in the confidence region are stable.

Coplanar Inclined Fits

Dynamical fitting can be sensitive to the true masses of planets, so it may be possible to determine the inclinations of the planetary orbits. In this subsection, we assume that the planets are on coplanar orbits. We allow the system’s orbital inclination relative to the sky plane to vary together with the other orbital parameters. and fitting parameters of the best fit as a function of inclination are shown in Fig. 5 and Fig. 6. In Fig. 5, the open circles are values, and the dashed lines represent 1, 2 and 3 confidence levels, which are of 1.0, 4.0 and 9.0 larger than the minimum . The changes slowly until the inclination drops below about . The decreasing stops at about , then increases with decreasing inclination. The minimum of at about inclination is lower than the value of for the edge-on best-fit model, which is a close to improvement over the edge-on best fit. The parameters of the best fit also change with the inclination as shown in Fig. 6, and there are inflection points near inclination for some fitting parameters. In Bayesian inference, an isotropic distribution of orbit normals implies a prior of for the inclination and an effective likelihood function , where . The equivalent quantity () is shown as a function of inclination by the filled triangles in the upper panel of Fig. 5. Although the prior reduces the effective likelihood when the orbit is highly inclined from the line of sight, the effective likelihood function also reaches a maximum at about inclination, similar to the original likelihood function.

The best fit is at near inclination, and its fitting parameters are listed in the third part of Table 2. The of 1.08 is close to 1, and the rms of is consistent with the estimated stellar jitter . The masses of the inner and outer planets are 4.78 and 4.80 respectively, where is the mass of Jupiter. The RV curve is shown in the left panel of Fig. 7. The best fit is dynamically stable for at least yr and is in the 2:1 MMR with both resonance angles, and , librating around , as shown in the right panel of Fig. 7. The major difference from the edge-on best fit is that the eccentricity of the outer planet is larger and the masses of the two planets are almost equal. As shown in the right panel of Fig. 7, the librating behavior of and is similar to that of the edge-on best fit, but the libration period is shorter and the average value of is slightly larger than those of the edge-on best fit.

Similar to the edge-on fits in § 4.1.1, we conduct a grid search around the coplanar inclined best fit allowing the inclination to vary. Fig. 8 shows the contours for the and grids. Similar to edge-on fits, only one minimum is found for each grid, and the inclinations of the fits near the minimum are not far from . However, the contours show discontinuities, which are not present in the edge-on fits. The discontinuities in and the free fitting parameters are illustrated in the left panel of Fig. 9 by extracting and inclination (one of the free fitting parameters) along the arrow (from to ) in the left panel of Fig. 8. The reason for the discontinuities is the appearance of a second local minimum with respect to the free fitting parameters when . The right panel of Fig. 9 shows as a function of inclination with and different fixed values taken along the arrow in the left panel of Fig. 8. For , there is a single minimum with at an inclination that increases with increasing . For , a second minimum with appears around . Because the starting condition of minimization is small inclination along the arrow, the fit is trapped in the minimum around when . The discontinuities in the grid shown in the right panel of Fig. 8 can be explained similarly.

Fig. 10 shows the dynamical properties of fits in the grid, with similar contents as Fig. 3 of edge-on fits. Similar to the edge-on case, the stable region is the region where there are libration amplitude contours (the thin solid black and thin dashed gray curves). The minimum is also far from the dynamical stability boundary, with libration amplitudes about and for and , respectively (see also Fig. 7). Most fits in the confidence region and almost all fits in the confidence region are in the stable region, different from the edge-on fits where there is still a small fraction of fits in the confidence region being unstable. All fits in the stable region are in 2:1 MMR with both and librating around , and the fit with the smallest libration amplitudes of both resonant angles () is near the confidence level. Fig. 11 shows the dynamical properties of fits in the grid. The dynamical properties and the statistics of this grid are similar to those of the grid.

Mutually Inclined Fits

Very few extrasolar planetary systems have successfully had their mutual inclinations measured by radial velocity alone. The most familiar case would be the GJ 876 planetary system, but even in this case the conclusions are not yet consistent. Based on the combination of RV and astrometry data, Bean & Seifahrt (2009) showed that the mutual inclination between the planets GJ 876 b and c is . Based on pure RV data, Correia et al. (2010) showed that the mutual inclination is less than , and an updated analysis by Baluev (2011) limited it by .

Since we are able to constrain the inclination of the HD 82943 system if we assume coplanar orbits, it is interesting to try to constrain the mutual inclination between the two planets in the system. We use the best-fit coplanar model of 2:1 MMR as an initial guess. After convergence of the LM method, reached 1.10 and the inclinations of the planets relative to the sky plane are about and , respectively. The becomes better, but the becomes worse after 2 more fitting parameters are introduced. Moreover, the fit is unstable on the order of 10,000 yr. A systematic parameter grid search show that we cannot constrain the mutual inclination. Fig. 12 shows the results of the grid search coupled with dynamical analysis. In the left panel, the thin curves are the contours with the thin dotted lines representing the minimum (near the lower left corner of the grid) and the 1 and confidence levels. The thick dashed line is the dynamical stability boundary, and the region in the middle of the grid is the stable region. The minimum located in the lower left corner of the grid is dynamically unstable and is shallow in parameter space. The right panel is an expansion of the lower left part of the grid. The solid curves represent contours and the thin dashed lines represent the contours of mutual inclination of and . Note that the mutual inclination is determined by and the mutual longitude of ascending node . The contours cross the region within and mutual inclination, showing that fits of different mutual inclination have the same quality and thus cannot be constrained by statistics. The stability test cannot constrain the mutual inclination either. In the grid searching, another local minimum was found with of 0.96. However this orbital configuration is extremely unstable because the system is highly mutually inclined with one of the “planets” having the mass of a brown dwarf. Based on the current data, the coplanar inclined best fit is already adequate, and it seems difficult to find better fits in mutually inclined configuration of 2:1 MMR.

Summary for 2:1 Resonance Fits

In summary, the coplanar inclined best fit results in significant improvement in statistics ( of the inclined best fit compared to of the edge-on best fit and of Keplerian best fit). Almost all fits inside the confidence region are stable, and all stable fits have both resonant angles and librating around . The stable 2:1 MMR configuration is robust for the HD 82943 system because all good fits are stable and in the 2:1 MMR. The fits with the smallest libration amplitudes of both resonant angles () are about from the best fit (with ) in the parameter grids, suggesting that the system does not favor small-libration-amplitudes configuration. On the other hand, we cannot solve for the mutual inclination for the HD 82943 system.

4.2 3-Planet Fits

The Laplace resonance configuration is well known as the double 2:1 MMR among the Galilean satellites Io, Europa, and Ganymede. In extrasolar planetary systems, the Laplace resonance may also play an important role in various MMR configurations. For example, the GJ 876 system (Rivera et al. 2010) and the HR8799 system (Marois et al. 2010) are suspected to contain planets in Laplace resonance. The existence of a third outer planet in a Laplace resonance with the two existing planets for the HD 82943 system was suggested by Beaugé et al. (2008). It is interesting to examine the viability of the 3-planet (Laplace resonance) fits with the new Keck data.

First, we input the residuals from the 2:1 MMR best fits of coplanar edge-on and inclined orbits to the Lomb-Scargle periodogram function (e.g., Press et al. 1992) implemented in the Systemic Console (Meschiari et al. 2009). The power spectra as a function of period are shown as the left and middle panels of Fig. 13. There is a peak at around 1100 days for both power spectra of edge-on and inclined orbits. However, both peaks have false alarm probability (FAP) larger than 10%, i.e., the peak is lower than the 10% FAP line which is the lowest dashed lines in the figure. Conventionally, the statement of having a new planet should be based on having a periodogram power spectrum peak at least higher than the line corresponding to (e.g., Marcy et al. 2005), which is the middle dashed line in the figure. The right panel of Fig. 13 shows the power spectrum of the window function which evaluates the periodicity contributed to the data from the choice of observational epochs. There is a peak at about 1100 days, which corresponds to an observational period of about 3 years. The coincidence for the three analysis in Fig. 13 all showing peaks at about 1100 days hints that the periodic signals in the residuals are partially due to the structured systematic noise.16 To verify the validity of the analytic FAPs calculated by Systemic, we have conducted a separate false alarm analysis using a complementary bootstrapping approach, similar to that of Wright et al. (2007) and Marcy et al. (2005). We randomly redrew the velocity residuals to the edge-on and inclined 2:1 resonant cases (with replacement), maintaining the temporal spacing of the observations, 1000 times. In each case we calculated the height of the tallest peak in the periodogram and compared to the tallest peak in the periodogram of the unscrambled data. In our unrestricted analysis (periods from 1 day to 10,000 days), we find that peaks near 1 day and 1100 days are tall, which result from the window function of this data set. Our FAP for the peak around 1100 days is 3.3% for edge-on case and 8.5% for inclined case, which are lower than the analytic FAP but still not comparable to . We then restrict our analysis to periods between 2 and 900 days, to avoid the tall peaks near 1 day and 1100 days. We find in the inclined case, 26% of our synthetic data sets had a peak taller than the tallest seen in the actual data set, indicating a FAP value . In the edge-on case, we find FAP . Thus, there is no evidence for the existence of a third planet in the Keck data.

Nevertheless, we try to fit a 3-planet model and look for any configurations associated with the Laplace resonance. Assuming coplanar edge-on orbits, the residuals from the coplanar edge-on best fit are treated as a new data set. A Keplerian orbit is fitted to the residuals with a period of about 1100 days. Then we input the initial guess of the 3-planet model from the sum of the edge-on two-planet best fit and the Keplerian fit of the third planet to our dynamical fitting code, allowing all parameters (except the fixed inclination ) to vary. A local minimum is found with a of 0.59 and rms of . The parameters of this local best fit are listed in Table 3 and its RV curve and residuals are shown in the left panel of Fig. 14. The period of the third planet is 1077 days, close to the peak in the periodogram, and is 0.402. This fit is dynamically unstable in a few hundred years as shown in the right panel of Fig. 14. Because of the high eccentricity of the third planet, its orbit is easily perturbed by the two massive planets inside and thus becomes unstable. The large uncertainties of some fitting parameters of the third planet from the covariance matrix suggest that there may be other minima. Starting from this local best fit, we search in the grid, and found two other local minima with slightly higher . So the fit in Table 3 is likely the global best fit of the 3-planet coplanar edge-on configuration. Similar to the best fit, the other two local minima are dynamically unstable in a few hundred years. In fact, dynamical stability test in the grid shows that all fits with are unstable. In order to find fits associated with Laplace resonance, we search the grids starting separately from the three minima found in the grid. We do not find local minimum near the nominal Laplace configuration, i.e., . Fits that are close to the nominal Laplace resonance configuration have , which is much larger than that of the 3-planet best fit, and they are dynamically unstable. Finally, we allow the inclination to float assuming coplanar orbits, and then fit the data starting with the coplanar edge-on best fit. Similar to the coplanar 2:1 MMR fits, the shows a minimum () at about .

In summary for the 3-planet fits, the periodograms of the residuals from the 2:1 MMR best fits do not provide evidence for the existence of a third planet. Additionally, the of the best fit is significantly lower than 1.0, which together with the goodness of the 2:1 MMR best fit (), hints that the 3-planet model results in over fitting the current data (i.e., the 3-planet model has too many parameters and its fit to the data is “too good”17). Finally, dynamical exploration shows that all good 3-planet best fits are dynamically unstable, and there is no good fit corresponding to the nominal Laplace resonance configuration.

4.3 1:1 Resonance Fits

Laughlin & Chambers (2002) pointed out that a 1:1 eccentric resonance configuration could be found in extrasolar planetary systems. The eccentric 1:1 resonant configuration can be generated by initially placing a planet in a circular orbit and the other planet in a highly eccentric orbit, with the period ratio nearly 1.0. The system maintains a stable configuration with angular momentum exchange between the two planets (i.e., the eccentricities are oscillating). Goździewski & Konacki (2006) reported a group of fits with 1:1 resonance that fit the combined RV data of CORALIE and Keck observations of the HD 82943 system almost as well as the 2:1 MMR fits.

Coplanar Fits

First, we explore edge-on coplanar 1:1 resonance fits. We use Keplerian fitting to search for an initial guess. We skip the strongest periodic signal of about 220 days in the original data set and directly fit with a Keplerian orbit at about 440 days. As we already know that the orbits may have high eccentricities, we force the first Keplerian orbit to have a relatively high eccentricity and then check the periodogram of the residuals. A periodicity of about 450 days is identified, and then the second Keplerian orbit with about 450 days is fit. The Keplerian best fit near the 1:1 resonance is adopted as an initial guess for dynamical fitting. The LM method quickly converges to a fit with of 1.69 and rms of . The parameters of this local best fit are listed in Table 4 with both and large ( and ). Similar to previous cases, we explore parameter grids around this local best fit to see if there are other minima. The results show that the local best fit is the only minimum in nearby parameter space for coplanar edge-on orbits. The best fit is unstable after several hundred years and dynamical stability analysis in the grids also shows that all considered fits in coplanar edge-on orbits are unstable after a short time. Based on the edge-on best fit, we vary the inclination to explore inclined coplanar fits. Unlike the 2:1 MMR and 3-planet cases, the of coplanar 1:1 resonance fits do not improve when the orbits are allow to be inclined, and we do not find any stable fit when all inclinations are explored. The results are similar to that of Goździewski & Konacki (2006), who did not find any stable fit for coplanar 1:1 resonance orbits.

Mutually Inclined Fits

Next we allow the orbits to be mutually inclined and adopt the best-fit model with a coplanar edge-on 1:1 resonant configuration as an initial guess. The LM algorithm found a local minimum with of 1.51 and rms of . This fit is reported as fit (a), and the fitting parameters are listed in Table 5. The mutual inclination of the orbits in fit (a) is about , and the fit becomes unstable quickly in hundreds of years. A grid search based on fit (a) does not find any better fit in nearby parameter space and dynamical analysis does not find any stable fit in these grids.

However, the exploration of mutually inclined fits in Goździewski & Konacki (2006) shows that for large mutual inclinations there may be multiple minima present in the parameter space, and some of the fits with large mutual longitude of ascending node () are stable. We adopt the mutually inclined stable best fit of Goździewski & Konacki (2006) as the initial guess for another exploration. After we adjusted the mean anomaly of their best fit to our initial epoch, we find a local minimum with of 1.61 and rms of and it is reported as fit (b), with parameters listed in Table 5. Fit (b) is highly mutually inclined with the mutual inclination of about , and becomes unstable in a few hundred years. However, grid search based on fit (b) yields other minima. Here we show a representative grid which starts from the fit (b) and yield two other minima in the grid. The left panel of Fig. 15 shows the contours , with the arrows pointing to the locations of three (potential) minima. The minimum in the lower left grid is labeled as fit (c) in the figure, whose parameters are listed in Table 5. The of 1.59 of fit (c) is slightly less than fit (b), and Fig. 16 shows the RV curve and residuals of fit (c) as well as its dynamical evolution. Interestingly, fit (c) is dynamically stable with the mass of planet 2 being about 37 times that of Jupiter, and with a small libration amplitude of as shown in the right panel of Fig. 16. Finally, the fit (d) which is pointed out in the left upper corner of Fig. 15 is not actually located in the grid, but is recognized from the tendency of the contour directions. We take a fit in the region of the grid where arrow (d) is pointing, and allow all parameters to vary. A local minimum is then found and it is reported as fit (d), with its parameters listed in Table 5. Its of 1.43 is the lowest among all best fits of 1:1 resonance, however, the orbital configuration of fit (d) is retrograde with the mutual inclination of about , and is unstable in less than a hundred years. The right panel of Fig. 15 shows the dynamical analysis in the grid. The thick lines are the dynamical stability boundary, and the stable region is the region with thin dashed lines. The thin dashed lines represent the libration amplitude of . All stable fits are in 1:1 resonance with librating around . Surprisingly, fits in the lower left corner of the grid, whose masses are much larger than Jupiter mass, are in a stable 1:1 resonance, whereas fits at high inclinations are unstable.

In summary for the 1:1 resonance fits, five local best fits have been found for 1:1 resonance, and stable 1:1 resonance fits have been found by grid search. However, the lowest among all 1:1 fits is 1.43, which is still much larger than 1.08 from the 2:1 MMR best fit. Therefore, there is no evidence supporting the 1:1 resonance model for the HD 82943 system.

5 Error Estimation

Error estimation is important because it provides an evaluation of the uncertainties in the planetary masses and orbital parameters in the best-fit model. Based on the analysis of §4, the coplanar inclined 2:1 MMR best fit is adopted as the orbital solution for the HD 82943 planetary system. Here we analyze the errors and distribution of fitting parameters for this best fit based on the bootstrap method.

The prescription in Press et al. (1992) is adopted as our bootstrap method for fitting parameter distribution estimation18. The bootstrap method uses the actual data set containing data points to generate synthetic data sets also with data points. Here, each synthetic data set consists of 64 entries, and each data entry is chosen randomly from all 64 entries in the real data set (Table 1). Each entry includes the observational time, radial velocity and instrumental uncertainty. Because of the random process, the synthetic data set almost certainly contains duplicated data points, i.e., they have the same observational time, radial velocity and uncertainty. For convenience, in the procedure of generating synthetic data set, when an entry of the real data set is chosen more than once, a random number of absolute value days is added to the observational time for every duplicated data point in the synthetic data set. We generate and fit 5000 samples to estimate the distribution of fitting parameters and calculate the 68.3 percent confidence interval for the model parameters by:


where is the probability density function as a function of , and , are the lower and upper value of 68.3 percent confidence errors, respectively.

The probability density distribution of fitting parameters determined by the bootstrap method is illustrated in Fig. 17. In the figure, the solid curves are the probability density distributions from all fits in bootstrap samples and the dashed lines in vertical direction represent the best-fit parameters. As shown in the figure, most parameters are centrally peaked, and some of them are asymmetrically distributed. Some parameters’ distributions do not peak at the best-fit parameters. For example, the peak of the distributions of and are slightly shifted from the best-fit parameters, hinting that these fitting parameters may be sensitive to some data points in the original data set. When some of the sensitive data points are absent from the synthetic data sets, the fitting results are slightly shifted from that of the original data set. Interestingly, the distribution of has double peaks, with a smaller one near about 0.05, hinting a small probability for small initial orbital configurations. All fits from bootstrap are integrated for a maximum time of 50,000 yr in order to examine whether dynamical stability will provide any constraints on the parameter distributions. As shown in Fig. 17, the dotted curves are probability density distributions from only the stable fits in bootstrap samples. The probability density distributions from the stable fits are similar to those from all fits, but almost all distributions from stable fits are more centrally peaked than those from all fits, meaning that although dynamical stability does not prefer a significantly different distribution, it does constrain the parameters better. In particular, the distributions of show a difference for small eccentricities. The small peak at vanishes after dynamical stability constraints. Thus the possibility of orbital configurations with small initial is ruled out by dynamical stability test. Finally, all stable fits are in 2:1 MMR with at least librating about as shown in the lower right panel of Fig. 17, where there are a few cases with circulating. The distributions of the libration amplitudes of both angles peak near the values of the best fit, and the overall range favor moderate libration amplitudes for both and .

Finally, the uncertainties in the orbital parameters for the 2:1 MMR coplanar best fit are listed in Table 6, as determined by three methods: the covariance matrix , the constant method. The inclination of the orbits is well constrained at degrees. Uncertainties determined by the bootstrap method are suggested as the reported uncertainties for the best-fit parameters of this system. The intervals of error bars determined by the constant method are comparable to those obtained from the covariance matrix. The errors from bootstrap are the largest amongst the three in every orbital fitting parameters. Both constant and bootstrap methods show asymmetric errors.

6 Comparisons With Self-Consistent Bayesian Analysis

6.1 Bayesian (DEMCMC) Approach

We also analyzed the radial velocity measurements using a Bayesian framework following Ford (2005) and Ford (2006). We assume priors that are uniform in the logarithm of orbital period, eccentricity, argument of pericenter, mean anomaly at epoch, and the velocity zero-point. For the velocity amplitude () and jitter (), we adopted a prior of the form with m s, i.e., high values are penalized (for a discussion of priors, see Ford & Gregory 2007). The likelihood for radial velocity terms assumes that each radial velocity observation is independent and normally distributed about the true radial velocity with a variance of , where is the internal measurement uncertainty, and is the jitter parameter.

We used an MCMC method based upon Keplerian orbit fitting to calculate a sample from the posterior distribution (Ford 2006). We calculated multiple Markov chains, each with states, and discarded the first half of the chains. We calculated Gelman-Rubin test statistics for each model parameter and several ancillary variables and found no indications of non-convergence amongst the individual chains. Finding no indications of non-convergence, we randomly choose a subsample () from the posterior distribution for further investigation.

Following the Keplerian fitting procedure, we use the method described in Payne & Ford (2011), Johnson et al. (2011) and Wang et al. (2012), using the subsample as the basis for a much more computationally demanding analysis that uses fully self-consistent -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. The priors for the model parameters are the same as those of the MCMC simulations.

For the -body integrations, we use a time symmetric 4th order Hermite integrator that has been optimized for planetary systems (Kokubo et al. 1998). We extract the radial velocity 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 of their starting value, and that no close-approaches occur within 0.1 times the semimajor axis during the 100-year -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 systems, for a total of -body integrations in each DEMCMC simulation.

Due to the very high computational cost of running large number of -body integrations, we confine the majority of our DEMCMC investigations to the coplanar fixed-inclination regime (with or ). We leave DEMCMC analysis of fits that allow the inclination to vary for future work.

6.2 Results of MCMC RV Analysis

We now present the results of our application of the MCMC methodology described in §6.1 to the RV data sets listed in Table 1, and compare the results to those from the minimization with grid search and the bootstrap method.

2-Planet Keplerian Fits

We illustrate a sample of the Keplerian MCMC analysis by plotting in Fig. 18 the periods of the two planets, which result from analyzing the Keck data sets. We show contour plots for the planetary periods () for the frequentist Levenberg-Marquardt (thin, black, solid) and Bayesian MCMC (thick, black, dotted) approaches. The thin solid black contours display the different levels (the minimum and according to from inside out). The thick dotted black contours display different information, as the MCMC algorithm provides us with a final density of solutions, so we plot iso-density contours containing , and of solutions. The 2:1 period ratio is plotted as a gray dotted line. We find that the two approaches agree very well, in the sense that the size and shape of contours are very similar and they both find consistent best-fit solutions very close to the 2:1 period ratio.

The probability density distribution of fitting parameters resulting from the Bayesian MCMC approach is plotted in Fig. 19 as the red curves, and the vertical lines represent the best-fit values from minimization. The majority of the parameters display smooth Gaussian profiles, and most of the peaks nearly coincide with the best-fit values. It should be noted that the eccentricity of the outer planet is rather poorly constrained in this Keplerian analysis. Using the bootstrap method described in §5, we generate 5,000 synthetic data sets and fit them with a Keplerian model. The distribution from bootstrap is shown as the black curves in Fig. 19, and many of them show less constraints (e.g., ) than the distribution from the MCMC approach. A full table of the mean values and their errors from both MCMC and bootstrap is given in Table 7. All the uncertainties from bootstrap are larger than those from the MCMC approach.

2-Planet Coplanar Inclined Fits: -Body

Implementing the -body fitting procedures detailed in §6.1, we arrive at the parameter fits illustrated by the blue contours in Fig. 18 for , with similar contents as the Keplerian fits. We find in Fig. 18 that the inclusion of mutual interactions leads to the period of the inner planet (at the first observing epoch) shifting to a slightly shorter period (from days, to days), while the outer period shifts to slightly larger values (from days, to days). Similar to Keplerian fits, the LM and DEMCMC approaches agree very well for inclined -body fitting results. The contours from -body fits are larger than those from Keplerian fits, suggesting that the Keplerian fits constrain the parameters better than the inclined -body fits.

We also compare the results from DEMCMC to those from bootstrap, as illustrated in Fig. 20, where the vertical lines represent the best-fit values from minimization from a coplanar inclined dynamical model, the red solid curves are from DEMCMC, and the black solid curves are from bootstrap. Almost all distributions from both methods show smooth gaussian profiles, and all of them peak around the best-fit values. Distributions from DEMCMC and bootstrap are similar in terms of both their shapes and sizes. The mean values and their errors are listed in Table 7. The stellar jitter is treated as an unknown parameter in Bayesian analysis. The distribution of jitter (not plotted) peaks at about with uncertainty of , which is consistent with the estimated value used in §4 and §5.

We also provide in Table 7 the mean values and errors of fitting parameters from both DEMCMC and bootstrap for coplanar edge-on () models, which are not discussed above, as a reference for future study. Our DEMCMC Bayesian algorithms have been tested in ever greater detail and found to perform successfully, with the strong overlap in the result from the different methods giving increased confidence in the robustness of our conclusions.

6.3 3-Planet Coplanar Edge-on Fits

In §4.2 we found using the periodograms and LM analysis that there was little evidence for a third planet. However, for the sake of completeness, we now model the system with three planets using the same MCMC approach (Keplerian and -body) outlined in §6.1.

We find that the best fit solution (not shown) has the inner two planets essentially unaltered at and days, with the third planet at days (although this is extremely poorly constrained, with significant uncertainties days), i.e., results that are very similar to those found using the LM approach. Hence we are confident that there is no evidence of the 1:2:4 resonance (requiring days) suggested in Beaugé et al. (2008).

7 Summary and Discussion

We have analyzed the orbital and dynamical state of the HD 82943 planetary system by dynamically fitting 10 years of Keck RV measurements. Based on parameter grid search, fits around the best fits as a function of various pairs of parameters have been systematically explored. Three type of fits associated with qualitatively different orbital configurations, the 2:1 MMR, 3-planet and 1:1 resonance configurations, have been examined.

In terms of the 2:1 MMR fits, our Keplerian best fit has = 1.38, significantly better than previous results () which were based on only 3.8 years of Keck RV data combined with the lower quality CORALIE data. The dynamical best fit of the coplanar edge-on orbits has = 1.21 and rms of , and it is in a 2:1 MMR with both resonance angles and librating around . Grid search coupled with dynamical stability test in the and grids shows that the best fit is the only minimum. The best fit is deep in the stable region and all fits in the stable region are in 2:1 MMR. When the inclination of coplanar orbits is varied, the as a function of inclination clearly shows a deep minimum at about , with of 1.08 and rms of , which is close to confidently better than the edge-on best fit. The inclined best fit contains two planets of masses 4.78 and 4.80 , and it is in 2:1 MMR with both and librating around . Systematic search for fits allowing the inclination to vary in and grids shows that the best fit is also the only minimum. All good fits are in the stable region and all stable fits are in 2:1 MMR. The contours and dynamical properties of fits in the grids are similar to that of coplanar edge-on fits, except that the contours in the grids show discontinuities. Finally, the mutual inclination of 2:1 MMR fits cannot be constrained by either statistics or dynamical stability test. Compared to previous fitting results of the 2:1 MMR configuration based on the lower-quality CORALIE and shorter Keck RV data (Lee et al. 2006; Goździewski & Konacki 2006; Beaugé et al. 2008), our 2:1 MMR best-fit model improves significantly in both and the rms. More importantly, assuming coplanar configuration, the inclination relative to the sky plane is well contained at about , and the system is stable in a 2:1 MMR configuration.

The periodograms of the residuals from both coplanar edge-on and inclined 2:1 MMR best fits do not show significant evidence for the existence of a third planet, contrary to previous results in Beaugé et al. (2008) where the periodogram of the residuals of their 2-planet best fit showed significant signal for the existence of a third planet. When we fit for a third planet, the best fit has of 0.59 and rms of . The fact that the is significantly lower than 1.0 and the rms is significantly lower than the estimated stellar jitter of hints that the 3-planet model over fits the current RV data. The best-fit model becomes unstable within hundreds of years due to the large . In the grid, two other local minima have been found, but they are unstable. Dynamical stability test in the grid shows that only fits with remain stable for 50,000 yr. Fits in grid have been explored, and we did not find any good fits associated with the Laplace resonance configuration.

For the 1:1 resonance configuration, only one minimum has been found in coplanar edge-on orbits, and several minima have been found in the mutually inclined fits. All coplanar fits and most of the mutually inclined fits are unstable, but we have found some stable fits with high mutual inclination in the grid. Only one minimum is stable with a small libration amplitude of , and dynamical behaviors of the stable fits are similar to the stable fits found by Goździewski & Konacki (2006). However, all fits we found have significantly higher () than the 2:1 MMR best fit, so they are ruled out.

In summary, based on the statistics and dynamical stability constraints, the 2:1 MMR coplanar inclined best fit is reported as the best fit for the HD 82943 planetary system. The best-fit parameters are listed in Table 6, and their uncertainties determined by bootstrap in Table 6 are suggested as the reported uncertainties. There is no evidence for either the 3-planet Laplace resonance fits or the 1:1 resonance fits. The HD 82943 planetary system contains two planets in 2:1 mean-motion resonance, and its dynamical state is well established. The resonant angles and of two nearly equal-mass planets are librating around with moderate libration amplitudes of about and , respectively.

It is interesting to show the differences between the Keplerian best fit, dynamical coplanar edge-on best fit and coplanar inclined best fit of 2:1 MMR in graphical form using radial velocity plots, rather than , so that one can have an intuitive evaluation of the improvements of the fitting. More importantly, if there are significant variations of the RV values from different fits after the last observed epoch of the Keck data sets, RV observations in the near future may provide more constraints on our best fit. A convenient method is to compare the residuals of two fits. For example, we plot the residuals of fit(a): , and then we plot a curve of the RV values of another fit (b) which is subtracted by the fit (a): . By evaluating how the curve fits the residuals compared to the zero line, we can know how fit(b) improves the fitting compared to fit(a).

First, we compare the Keplerian best fit to the coplanar inclined best fit. The upper panel of Fig. 21 shows the differences of the residuals from the Keplerian and the inclined best fit, in which the dots are the residuals of the Keplerian best fit and the curve represents the RV values of [fit()fit(Kep)]. The fluctuations of the curve oscillate around the zero line and do not show an obvious systematic trend in the observation time span (10 yr). This situation suggest that the improvement from the Keplerian best fit to the coplanar inclined best fit is primarily contributed by the short term mutual interactions of the planets. A large fraction of dots are apparently fitted better by the curve than the zero line. For example, the first dot and dots around BJD 2,454,500 and BJD 2,455,500 significantly deviate from zero but they are much closer to the curve. The orbital precession rate is on the order of only about (or in 10 yr) for the coplanar inclined best fit, so it is reasonable that we primarily see improvement from short term interactions. The fluctuations in the near future do not show significant peaks until BJD 2,458,000. Next we compare the coplanar edge-on best fit to the coplanar inclined best fit in the lower panel of Fig. 21, in the same format as the upper panel. In this comparison, the fluctuations are smaller than in the comparison associated with the Keplerian best fit in the observed time span, except for the large peaks at the beginning. Unlike the comparison with the Keplerian best fit, the improvements are not so obvious as a large fraction of the dots are not obviously fitted well by the curve. Interestingly, the curve shows large peaks at around BJD 2,456,800, which is about 3 yr after the last observed epoch in the Keck data. Thus future RV observations at around that time could provide more constraints on the inclinations and the true planetary masses of the HD 82943 planetary system.

During the course of the submission and review of this article, we learned of a complimentary investigation by Kennedy et al (in press) in which Herschel observations of HD 82943 detect a debris-disk with an inner edge (far beyond the planets studied in our analysis). The debris-disk appears to have a best-fit inclination of approximately degrees (to the plane of the sky), strongly supporting the inclination we deduced from our purely dynamical studies.

We thank NASA and NExScI for providing Keck time in the 2011A semester for the study of multiplanet systems (NExScI ID40/Keck ID# N141Hr). Data presented herein were obtained at the W. M. Keck Observatory from telescope time allocated to the National Aeronautics and Space Administration through the agency’s scientific partnership with the California Institute of Technology and the University of California. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. We thank Howard Isaacson, R.P. Butler, and S. Vogt for help with observing at the telescope, and the referee for helpful comments on the manuscript. X.T. and M.H.L. were supported in part by the Hong Kong RGC grant HKU 7034/09P. Contributions by M.J.P. and E.B.F. were supported by NASA Origins of Solar Systems grant NNX09AB35G only prior to August 8, 2011. J.A.J. was supported by generous grants from the David and Lucile Packard Foundation and the Alfred P. Sloan Foundation. J.T.W. was supported by NSF Astronomy and Astrophysics Grant AST-1211441. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. M.H.L. and E.B.F. also acknowledge the hospitality of the Isaac Newton Institute for Mathematical Sciences, where part of this work was completed.
JD Radial Velocity Uncertainty
(2450000) () ()
Data Set 1
2006.9130 43.90 1.56
2219.1210 22.08 1.36
2236.1262 28.57 1.34
2243.1295 36.80 1.35
2307.8391 45.39 1.51
2332.9834 9.84 1.72
2333.9558 10.83 1.63
2334.8726 0.77 1.49
2362.9717 25.84 1.66
2389.9438 47.74 1.57
2445.7387 55.28 1.55
2573.1473 49.82 1.42
2575.1400 47.46 1.36
2576.1437 51.51 1.55
2601.0664 24.35 1.57
2602.0728 17.51 1.47
2652.0009 19.62 1.55
2988.1092 91.26 1.40
3073.9287 2.98 1.56
3153.7544 0.00 1.41
3180.7448 62.77 1.36
3181.7416 53.85 1.51
Data Set 2
3397.9083 107.66 1.13
3480.7572 33.53 0.97
4084.0897 28.82 0.88
4139.0165 32.21 1.06
4398.1361 3.17 1.05
4428.1286 23.07 1.15
4454.1151 39.44 0.96
4455.0271 40.27 1.05
4456.0837 40.81 0.98
4492.9843 49.24 1.31
4544.9394 5.61 0.67
4545.9278 7.24 0.67
4602.8125 52.94 1.27
4603.7880 48.85 1.25
4807.1705 25.89 1.37
4811.1235 21.83 1.18
4847.0406 10.63 1.41
4929.8077 44.36 1.41
4963.8378 19.10 1.33
4985.7698 3.45 1.21
5134.0952 25.88 1.20
5172.1518 99.47 1.19
5174.0713 100.61 1.13
5189.1241 87.31 1.10
5197.9924 81.09 1.26
5229.0603 49.61 1.31
5252.0320 32.41 1.41
5255.8658 25.32 1.23
5285.8518 0.23 1.25
5289.8818 13.11 1.25
5312.7844 27.11 1.18
5320.7774 35.04 1.13
5342.7496 19.90 1.18
5372.7395 53.19 1.26
5522.0604 47.18 1.27
5556.0871 36.55 1.29
5558.0063 26.56 1.13
5585.0690 74.56 1.32
5605.9988 99.45 1.39
5633.7992 74.82 1.26
5672.8239 43.12 1.20
5700.7478 16.85 1.27
Table 1: Radial Velocities of HD 82943 from Keck
Parameter Inner Planet Outer Planet
Keplerian fit: ,
58.5 38.0
(days) 220.2 440.6
0.413 0.136
(deg) 114 52
(deg) 273 22
() 1.72 1.53
Coplanar edge-on dynamical fit: ,
57.8 37.8
(days) 221.0 438.5
0.407 0.096
(deg) 123 90
(deg) 276 342
() 1.71 1.53
Coplanar inclined dynamical fit: ,
54.3 39.8
(days) 219.3 442.4
0.425 0.203
(deg) 133 107
(deg) 256 333
() 4.78 4.80
(deg) 19.4
Table 2: Best fits of 2:1 MMR configuration
Parameter Planet 1 Planet 2 Planet 3
58.8 38.1 6.3
(days) 221.2 438.7 1077.4
0.412 0.068 0.402
(deg) 116 89 27
(deg) 279 343 148
() 1.73 1.55 0.314
Table 3: Best fit of edge-on 3-planet configuration
Parameter Planet 1 Planet 2
92.3 59.7
(days) 442.0 439.5
0.466 0.654
(deg) 117 125
(deg) 325 136
() 3.33 1.83
Table 4: Best fit of coplanar edge-on 1:1 resonance configuration
Parameter Planet 1 Planet 2
Fit (a): ,
95.0 63.0
(days) 450.5 436.7
0.516 0.678
(deg) 126 135
(deg) 322 129
(deg) 110.0 10.4
(deg) 0.0 327
() 3.58 10.47
Fit (b): ,
88.0 64.5
(days) 446.4 439.9
0.460 0.701
(deg) 111 130
(deg) 331 130
(deg) 88.1 8.1
(deg) 0.0 165.1
() 3.22 13.31
Fit (c): ,
89.0 65.9
(days) 455.3 442.0
0.480 0.641
(deg) 135 156
(deg) 317 122
(deg) 32.8 3.2
(deg) 0.0 161.1
() 6.07 37.07
Fit (d): ,
93.6 83.0
(days) 440.4 451.0
0.508 0.776
(deg) 125 141
(deg) 318 131
(deg) 10.0 162.3
(deg) 0.0 71
() 19.12 6.98
Table 5: Best fits of mutually inclined 1:1 resonance configuration
Parameter Value Uncertainties (1-)
Constant Bootstrap
() 54.4
(days) 219.3
(deg) 133
(deg) 256
() 39.8
(days) 442.4
(deg) 107
(deg) 333
(deg) 19.41
(AU) 0.7423
(AU) 1.1866
() 4.78
() 4.80

Note. – Semimajor axis and planetary mass are not direct fitting parameters in our model, so their uncertainties here are only determined by the bootstrap method. The parameters listed here are suggested as the best-fit parameters for the HD 82943 planetary system, and the uncertainties determined by bootstrap method are suggested as the reported uncertainties.

Table 6: Uncertainties of fitting parameters for the coplanar inclined 2:1 MMR best fit
Parameter Mean
Keplerian Fits Inclined Dynamical Fits Edge-on Dynamical Fits
Bootstrap MCMC Bootstrap DEMCMC Bootstrap DEMCMC