A publishing partnership

Articles

A SECOND GIANT PLANET IN 3:2 MEAN-MOTION RESONANCE IN THE HD 204313 SYSTEM

, , , , , , , , , and

Published 2012 July 5 © 2012. The American Astronomical Society. All rights reserved.
, , Citation Paul Robertson et al 2012 ApJ 754 50 DOI 10.1088/0004-637X/754/1/50

0004-637X/754/1/50

ABSTRACT

We present eight years of high-precision radial velocity (RV) data for HD 204313 from the 2.7 m Harlan J. Smith Telescope at McDonald Observatory. The star is known to have a giant planet (Msin i = 3.5 MJ) on a ∼1900 day orbit, and a Neptune-mass planet at 0.2 AU. Using our own data in combination with the published CORALIE RVs of Ségransan et al., we discover an outer Jovian (Msin i = 1.6 MJ) planet with P ∼ 2800 days. Our orbital fit suggests that the planets are in a 3:2 mean motion resonance, which would potentially affect their stability. We perform a detailed stability analysis and verify that the planets must be in resonance.

Export citation and abstract BibTeX RIS

1. INTRODUCTION

HD 204313, a Sun-like V = 8 star observable in both hemispheres, has been a target for multiple radial velocity (RV) surveys. Ségransan et al. (2010) announced the detection of the first member of the star's planetary system with the discovery of HD 204313b, a Jovian-class (Msin i ∼ 4 MJ) planet on a long-period (P ∼ 5 yr) orbit. More recently, the HARPS survey revealed an interior Neptune-mass planet with P = 35 days (Mayor et al. 2011).

At [Fe/H] = 0.18 (as measured by Ségransan et al. 2010), HD 204313 follows the observed trend of gas giant hosts being generally metal-rich (e.g., Fischer & Valenti 2005). Furthermore, planet c adds to the mounting evidence that Neptune- and lower-mass planets are extremely common around main-sequence stars (Howard et al. 2010; Mayor et al. 2011; Wittenmyer et al. 2011a). In many ways, HD 204313 represents a "typical" planetary system, according to current observations.

Since 1987, we have used the 2.7 m Harlan J. Smith Telescope at McDonald Observatory for a long-baseline RV planet survey (Cochran & Hatzes 1993). An upgrade to our two-dimensional coudé spectrograph in 1998 gave us access to the full optical wavelength range of our I2 absorption cell, enabling us to monitor hundreds of FGK stars with ∼6 m s−1 precision over 7–13 years. One of the primary scientific objectives of the survey is to obtain a census of Jupiter analogs—giant planets in long-period orbits (see Wittenmyer et al. 2006, 2011b, for a complete discussion of Jupiter analogs and early detection limits from the McDonald Observatory RV survey). We have recently announced three giant planets in long-period orbits (Robertson et al. 2012), demonstrating that we have the time baseline and sensitivity to detect long-period giants. In the core-accretion theory of giant planet formation (Pollack et al. 1996; Lissauer 1995), surface-density enhancement by ices facilitates the formation of ∼10–15 M cores. The ice line, beyond which ices are present in the protoplanetary disk, has been estimated to lie at 1.6–1.8 AU in a minimum-mass solar nebula (Lecar et al. 2006). For the case of HD 204313, the inclusion of published CORALIE velocities from Ségransan et al. (2010) gives us a total time baseline of 12 years, extending our sensitivity comfortably beyond the ice line into the formation locations of gas giant planets. In this paper, we present HD 204313d, another Jupiter analog exterior to planet b, and describe its orbital parameters and evolution.

2. OBSERVATIONS AND DATA REDUCTION

Our RV data for HD 204313 are all taken from the 2.7 m Smith telescope between 2003 July and 2011 June, resulting in an eight-year time baseline. We use the Tull coudé spectrograph (Tull et al. 1995) with a 1farcs8 slit, yielding a resolving power R = 60,000. Our RV measurement procedure and reduction code AUSTRAL is discussed in detail in Endl et al. (2000). In short, immediately before starlight enters the slit, it passes through an I2 absorption cell regulated at 50°C, which superimposes thousands of molecular absorption lines over the object spectra in the spectral region between 5000 and 6400 Å. Using these lines as a wavelength standard, we simultaneously model the time-variant instrumental profile and Doppler shift relative to an I2-free template spectrum. The resulting RVs are corrected for the motion of the observatory around the solar system barycenter. We report our RV data for HD 204313 in Table 1.

Table 1. Radial Velocities for HD 204313

BJD−2,450,000 Radial Velocity Uncertainty SMW
  (m s−1) (m s−1)  
2840.892291 −8.48 7.26 0.1799 ± 0.0215
2933.666733 −11.81 4.58 0.1757 ± 0.0231
3930.951397 25.41 4.08 0.1725 ± 0.0262
4347.740522 76.44 6.37 0.2100 ± 0.0327
4376.744143 88.50 4.07 0.1849 ± 0.0217
4401.617979 86.26 4.55 0.1769 ± 0.0250
4606.945812 62.97 5.54 0.1499 ± 0.0233
4663.924822 46.99 5.37 0.1728 ± 0.0251
4663.938000 39.25 2.76 0.2807 ± 0.0466
4703.749028 44.40 5.24 0.1592 ± 0.0206
4732.798800 40.00 3.96 0.1590 ± 0.0214
4752.702889 40.71 6.14 0.1677 ± 0.0217
4782.617353 32.02 6.04 0.1706 ± 0.0228
4986.953099 −3.09 5.05 0.1604 ± 0.0265
5023.896398 −10.58 6.18 0.1640 ± 0.0233
5049.798314 −14.30 6.87 0.1725 ± 0.0246
5101.700768 −23.34 4.92 0.1682 ± 0.0223
5172.538410 −20.96 7.60 0.1688 ± 0.0250
5470.683390 −51.54 4.11 0.1479 ± 0.0222
5496.568795 −54.50 3.56 0.1583 ± 0.0275
5526.556625 −54.25 6.77 0.2355 ± 0.0439
5529.556590 −51.40 5.55 0.1704 ± 0.0290
5548.544742 −46.72 3.67 0.1508 ± 0.0232
5722.947515 −30.38 7.86 0.1689 ± 0.0294
5761.927802 −18.40 5.41 0.2000 ± 0.0296
5790.814330 −16.89 4.45 0.2248 ± 0.0355
5791.852748 −16.59 4.50 0.1901 ± 0.0282
5811.740128 −20.36 6.16 0.1496 ± 0.0251
5812.750735 −25.75 3.56 0.1498 ± 0.0229
5817.780745 −33.40 6.91 0.2204 ± 0.0354
5838.671924 −13.10 5.43 0.1492 ± 0.0231
5840.686936 −7.06 3.68 0.1768 ± 0.0259
5841.664865 −10.98 5.72 0.1626 ± 0.0346
5842.697923 −13.76 5.16 0.1673 ± 0.0367
5845.589786 −7.23 4.97 0.1803 ± 0.0262
5846.724380 −18.09 3.57 0.2182 ± 0.0372

Download table as:  ASCIITypeset image

3. STELLAR CHARACTERIZATION

We seek to independently verify the stellar atmosphere parameters for HD 204313 derived by Ségransan et al. (2010). Using our I2-free stellar template, we measure the equivalent widths of 61 Fe i lines and 17 Fe ii lines. We feed these equivalent widths into the MOOG5 local thermodynamic equilibrium (LTE) line analysis and spectral synthesis program (Sneden 1973). By utilizing a grid of ATLAS9 model atmospheres (Kurucz 1993), MOOG derives heavy-element abundances to match the measured equivalent widths. We then determine effective temperature Teff by removing any trends in abundances versus excitation potential (assuming excitation equilibrium), and compute microturbulent velocity ξ by eliminating trends with reduced equivalent width (≡ Wλ/λ). Stellar surface gravity is obtained by forcing the abundances measured with Fe i and Fe ii lines to match (assuming ionization equilibrium). Our measured abundances are differential with respect to the Sun. Using a solar port, we have taken a solar spectrum using the same instrumental setup used for our RV observations, and run the above analysis for the Sun. For reference, we obtain values of Teff = 5780 ± 70 K, log g = 4.50 ± 0.09 dex, ξ = 1.16 ± 0.06 km s−1, and log epsilon (Fe) = 7.52 ± 0.05 dex for the Sun. Full details of our stellar analysis can be found in Brugamyer et al. (2011).

Our stellar parameters for HD 204313 are given in Table 2. For values our routine does not calculate, we include catalog values from Kharchenko & Roeser (2009, Version 3) and Casagrande et al. (2011). Our computed values agree extremely well with those presented in Ségransan et al. (2010). Of particular interest is our measured [Fe/H] of 0.24 ± 0.06, which confirms the metal-rich nature of HD 204313.

Table 2. Stellar Properties for HD 204313

Spectral Type G5 V
Va 8.006 ± 0.014
BVa 0.695 ± 0.02
MV 4.63 ± 0.03
Parallaxa 21.06 ± 1.04 mas
Distance 47 ± 0.3 pc
Teff 5760 ± 100 K
log g 4.45 ± 0.12
[Fe/H] 0.24 ± 0.06
ξ 1.20 ± 0.15 km s−1
Massb 1.02 M
Ageb 7.20 Gyr
log R'HK −4.65 ± 0.03

Notes. aKharchenko & Roeser (2009). bFrom Casagrande et al. (2011), maximum likelihood estimate using Padova isochrones.

Download table as:  ASCIITypeset image

4. ORBIT MODELING

Over the eight-year period from 2003 July to 2011 June, we collected 36 RV points for HD 204313. Our data are plotted as a time series in Figure 1(a). The rms scatter about the mean of these velocities is 40 m s−1, with an average internal error of 5.21 m s−1. When analyzed alone, our RVs show the high-amplitude (K ∼ 70 m s−1) signal expected as a result of HD 204313b, with a period around 5.5 years. We compute Keplerian orbital solutions using the GaussFit modeling program (Jefferys et al. 1988) and the SYSTEMIC console (Meschiari et al. 2009), finding excellent agreement between the one-planet solutions from both routines. However, a one-planet fit to our data gives a period more than 100 days longer than the period reported in Ségransan et al. (2010), a discrepancy more than three times the combined 1σ uncertainties in the orbital period for the two models. Additionally, our fit includes a long-period linear trend. We note that Mayor et al. (2011) also find a period for planet b considerably longer than the originally published value.

Figure 1.

Figure 1. (a) Top: radial velocity data for HD 204313. Points in black are our 2.7m observations, while points in red are CORALIE observations from Ségransan et al. (2010). The best-fit orbit model is shown as a blue line. Bottom: residuals to a two-planet fit. (b) RVs after subtracting our fit to planet d from the velocities in (a). The blue line shows our Keplerian model for planet b. (c) RVs after subtracting our fit to planet b from the velocities in (a). The blue line shows our Keplerian model for planet d.

Standard image High-resolution image

We then compute a one-planet model using the CORALIE RVs as well as our own. The combined RV set includes 132 RVs taken over 10 years. The resulting parameters are closer to the previously published solution; we find a period of 2000 days, with eccentricity 0.16 and a minimum mass of 4.36 Jupiter masses. However, this fit is still considerably discrepant from the Ségransan et al. (2010) solution, and leaves a residual rms scatter of 11.0 m s−1, and a reduced χ2 = 5.66.

We have computed the fully generalized Lomb–Scargle periodogram (Zechmeister & Kürster 2009) for the combined CORALIE-McDonald data set, and the residual RVs around the one-planet fit. The resulting power spectrum is shown in Figure 2. The periodogram of the residual RVs shows significant peaks around 340 days, 395 days, and a broad peak between 2700 and 6700 days. We calculate a false-alarm probability (FAP) for these peaks using the method described in Sturrock & Scargle (2010), and find an FAP of approximately 5 × 10−5 for the long-period peak, while the 395 day and 340 day peaks have FAPs of 4 × 10−4 and 1.5 × 10−3, respectively.

Figure 2.

Figure 2. From top: (a) Generalized Lomb–Scargle periodogram for the combined CORALIE/McDonald RVs of HD 204313. (b) The same periodogram for the residual RVs after subtracting a one-planet fit. (c) Periodogram of the residual RVs after subtracting a two-planet fit. (d) Periodogram of our time sampling (the window function). The dashed lines indicate the approximate power level for an FAP of 0.01, computed from Equation (24) of Zechmeister & Kürster (2009).

Standard image High-resolution image

We attempt to fit an additional planet to the residuals at each of the periods identified in the periodogram. Our fitting routine produces unsatisfactory solutions at the two shorter periods, but converges to a fit with an outer giant (Msin i = 1.68 MJ) planet at 2831 days. The period of planet b in this solution is consistent with the original published result, although the eccentricity is higher than in a one-planet fit. The resultant two-planet fit is included in Figure 1(a), and we give residual plots of the individual planets in Figures 1(b) and (c). As with the one-planet solution, GaussFit and SYSTEMIC agree nicely on the orbital parameters and their uncertainties. The parameters of our final orbital model are listed in Table 3. The addition of planet d removes the need to include a linear slope. Although we attempted to fit an outer planet to each of the CORALIE and McDonald RV sets individually, our routines failed to converge for either set. Evidently, both data sets are required to achieve the time baseline needed to detect planet d, a fact reinforced by our periodogram analysis. When examining the residual RVs to our one-planet fit for each data set individually, we see only a monotonic increase in power at long periods for the CORALIE data and insignificant power in the McDonald data. Only when the data are combined, and the total time baseline exceeds a full orbit of planet d, does the power spectrum show a clearly defined peak around the period of that planet.

Table 3. Two-planet Orbital Solution for the HD 204313 System

Orbital Parameter Planet b Planet d
Period P (days) 1920.1 ± 25 2831.6 ± 150
Periastron passage T0 (BJD−2,450,000) 2111.6 ± 28 6376.9 ± 176
RV amplitude K (m s−1) 57.0 ± 3 23.7 ± 4
Mean anomaly M0a 300° ± 0fdg4 137° ± 2°
Eccentricity e 0.23 ± 0.04 0.28 ± 0.09
Longitude of periastron ω 298° ± 6° 247° ± 16°
Semimajor axis a (AU) 3.04 ± 0.06 3.93 ± 0.14
Minimum mass Msin i (MJ) 3.55 ± 0.2 1.68 ± 0.3
CORALIE RV offset (m s−1) −19.3
2.7 m RV offset (m s−1) 29.8
rms (m s−1) 7.80
Stellar "jitter" (m s−1) 5.46

Note. aEvaluated at the time of the first RV point reported in Ségransan et al. (2010).

Download table as:  ASCIITypeset image

We note that we have not included planet c (Mayor et al. 2011) in our analysis. With a reported RV amplitude of just 3.28 m s−1, the signal of this short-period planet is below the sensitivity limit of our data and that of CORALIE. Indeed, our periodogram of the residuals to our two-planet solution (Figure 2) shows no additional signals. Furthermore, the inclusion of a third planet with the orbital elements published for planet c does not significantly change our orbital solution. However, the rms (7.79 m s−1) and reduced χ2 (2.98) of our two-planet model are still higher than we expect given the precision of the Tull spectrograph and the CORALIE data. While this is reflected in our model as a relatively high level of stellar "jitter" (5.46 m s−1), our stellar activity analysis (see below) suggests that HD 204313 should not be so active. Although it would be ideal to include the orbit of planet c in our model to verify this hypothesis, we are unable to do so because Mayor et al. (2011) include neither their measured RVs nor their complete orbital fit for the HD 204313 system. We nevertheless conclude that the additional scatter around our fit is most likely due to the unresolved planet c, and potentially additional low-mass companions. We refer to the outer planet as HD 204313d in acknowledgement of the inner Neptune-mass planet.

Eggenberger et al. (2007) report a companion star 6farcs2 to the south of HD 204313, although they admit a significant probability of a chance alignment. At a distance of 47 pc, the angular separation indicates a minimum physical distance of 583 AU between the two objects. The companion is approximately 9 mag fainter in the near infrared (Eggenberger et al. 2007), and is therefore much less massive than HD 204313 if they are in fact bound. If we overestimate the mass of this object at 0.5 M and assume it is associated with HD 204313, the resulting RV slope due to the companion is 0.28 m s−1 yr−1, which is roughly equal to our 1σ uncertainty level of 0.2 m s−1 yr−1 for a slope in the combined data set. It is therefore safe to conclude that the second star is not influencing our modeling of the planetary system.

5. STELLAR ACTIVITY AND LINE BISECTOR ANALYSIS

While we do not anticipate that stellar activity should produce RV signals of the amplitudes of planets b and d, it is nevertheless important to understand how changes in the atmosphere of HD 204313 may influence our velocity measurements, particularly with the amount of scatter seen around our fit. We examine stellar activity simultaneously with RV through line bisector analysis of stellar lines outside the I2 region and changes in the Ca H and K indices.

Changes in the stellar photosphere (starspots, etc.) may produce changes in the measured RVs. However, these processes will also alter the shapes of the individual stellar absorption lines. Following the method of Brown et al. (2008), we calculate the bisector velocity span (BVS) for each of our spectra. The BVS is sensitive to these subtle changes in line shapes, and therefore a reliable indicator of activity the stellar photosphere.

Similarly, if stellar activity is producing RV signals, those signals should also appear in the Ca H and K indices. For each RV point in Table 1, we have computed the Mount Wilson SHK index, which we list alongside the velocities. From Noyes et al. (1984) we use SHK to derive log R'HK, the ratio of Ca H and K emission to the bolometric luminosity of the star. From log R'HK we obtain a more general idea of the overall activity level of HD 204313.

All examinations show HD 204313 to be an extremely quiet star. The results of our activity analyses are shown in Figure 3. In Figure 3(a), we plot BVS and SHK versus our measured RVs and their residuals around the one-planet fit. In both cases, there is no significant correlation, suggesting that photospheric activity is not influencing our velocities. Periodograms of SHK and BVS (Figure 3(b)) show no periodicity for either index. Furthermore, we measure an rms of only 17 m s−1 for the BVS, and log R'HK = −4.65. It is not surprising, then, that we see no signals or correlations in any of our activity indicators. With three planets now known, HD 204313 is rapidly becoming a rich planetary system. Its low activity level makes it an ideal candidate for follow-up observations to search for additional low-mass companions.

Figure 3.

Figure 3. (a) Left: bisector velocity spans plotted against our measured RVs (top) and residual RVs to a one-planet fit (bottom) for HD 204313. Right: SHK indices plotted against our measured RVs (top) and residual RVs to a one-planet fit (bottom) for HD 204313. (b) Generalized Lomb–Scargle periodograms for the BVS (top) and SHK indices (bottom) of our spectra for HD 204313. The dashed lines indicate the approximate power level for an FAP of 0.01.

Standard image High-resolution image

6. DYNAMICAL STABILITY ANALYSIS

A number of recent studies have highlighted the need for observational detections of multiple exoplanet systems to be supported by dynamical simulations that test whether the orbits of the proposed planets are dynamically feasible (e.g., Horner et al. 2011, 2012b; Wittenmyer et al. 2012; Hinse et al. 2012). Such studies are particularly important when the planets in question appear to move on orbits close to mutual mean-motion resonance (e.g., Robertson et al. 2012), an architecture that can yield either extreme stability or instability, depending on the precise orbits of the planets involved. In the case of HD 204313, the best-fit orbits for planets b and d suggest that they may well be trapped in mutual 3:2 mean-motion resonance. Furthermore, as shown in Figure 4, the planets' orbital paths nearly cross, which could potentially lead to collisions. As such, we chose to perform a highly detailed dynamical study of the orbits of planets b and d to investigate whether the orbits that best fit the data are dynamically feasible.

Figure 4.

Figure 4. Face-on orbital diagram of the giant planets in the HD 204313 system. The ellipses shown are derived from the model in Table 3 (the open square in Figure 5), with the lines from the star pointing toward the periastron of each planet. The locations of planets b and d are adopted from the mean anomalies in Table 3.

Standard image High-resolution image

Following earlier work (Horner et al. 2011; Marshall et al. 2010; Robertson et al. 2012), we used the Hybrid integrator within the n-body dynamics package MERCURY (Chambers 1999) to examine test systems in which the initial orbit of planet b was held fixed at the nominal best-fit values (in this case, a = 3.04 AU, e = 0.23). The initial orbit of planet d was then systematically changed from one simulation to the next, such that scenarios were tested for orbits spanning the full ±3σ error ranges in semimajor axis, eccentricity, longitude of periastron and mean anomaly. Such tests have already proven critical in confirming or rejecting planets thought to follow unusual orbits (e.g., Horner et al. 2011; Wittenmyer et al. 2012), and allow the construction of detailed dynamical maps for the planetary system studied, in orbital element phase space.

We examined 31 unique values of semimajor axis for planet d, ranging from 3.51 AU to 4.35 AU, inclusive, in even steps. For each of these 31 initial semimajor axes, we studied 31 values of orbital eccentricity, ranging across the full ±3σ range (e = 0.01–0.55). For each of the resulting 961 ae pairs, we considered 11 values of initial longitude of periastron (ω), and 5 values of initial mean anomaly (M0), resulting in a total suite of 52,855 (31 × 31 × 11 × 5) plausible architectures for the HD 204313 system.

In each of these simulations, the masses of the two planets studied were set to their minimum (Msin i) values. The mass of planet b was therefore set to 3.55 MJ, while that of planet d was set to 1.68 MJ. To first order, the more massive the planets, the more strongly they will perturb one another, and so setting their masses to the minimum allows us to maximize the potential stability of the planetary system. In other words, we expect our resulting dynamical maps to show the maximal stability of the orbits tested. The dynamical evolution of the two planets was then followed for a period of 100 million years, or until one of the planets either collided with the central star, was transferred to an orbit that took it to a distance of at least 10 AU from the central star, or collided with the other planet. Collisions were modeled by assuming a density of 1.33 g cm−3—equal to the average density of Jupiter—for each planet and computing a radius accordingly, so that our code registered a collision if an actual physical encounter occurred. The time of such events was recorded, allowing us to construct a dynamical map of the planetary system, shown in Figure 5(a). The figure shows the mean lifetime of the HD 204313b-d system as a function of the initial semimajor axis, a, and eccentricity, e, of planet d. Each individual initial ae pair was tested a total of 55 times, each of which featured a different initial combination of ω–M0. The lifetimes shown are the mean value of the 55 individual lifetimes obtained from those runs.

Figure 5.

Figure 5. (a) The mean dynamical lifetime of the HD 204313b–d planetary system, as a function of the initial semimajor axis (a) and eccentricity (e) of planet d. The lifetimes are shown on a logarithmic scale, ranging from 102 years (blue) to 108 years (red). The location of the nominal best-fit orbit for planet d is denoted by the open square, with the 1σ uncertainties shown by the solid lines radiating from that point. It is immediately obvious that the vast majority of the ae space tested is highly unstable, with only a narrow region of stability centered on the mutual 3:2 mean-motion resonance between planets b and d. (b) The mean lifetime of the HD 204313b–d system, as a function of the semimajor axis, a, and longitude of periastron, ω, of planet d's orbit. The lifetime shown at each location in a–ω space is the mean of 155 individual runs, which tested 31 different orbital eccentricities and 5 different mean anomalies for that particular a–ω combination. The color scheme is the same as in (a), and the nominal best-fit orbit for planet d is again marked by the unfilled square, with the 1σ errors on that fit shown by the solid lines that radiate from that point.

Standard image High-resolution image

Aside from resonant solutions, the entire ae phase space of allowed orbits for planet d is extremely unstable, with collisions between planets b and d often occurring within the first few hundred years of the simulations. It is also clear that, even within the resonance, some subset of the solutions are dynamically unstable (hence, the reason the mean lifetime in the stable region is somewhat less than 108 years). At the highest eccentricities permitted for the orbit of planet d, no stable solutions exist, but there is a broad region of stability within the 1σ errors on the best-fit orbit. As can be seen in Figure 5(b), the stability of orbits in the vicinity of the 3:2 mutual mean-motion resonance between the planets is a strong function of the longitude of periastron ω for planet d (we note here that planet b's initial longitude of periastron was 298 deg). Qualitatively, the strong ω dependence is reflective of the fact that the 3:2 resonance provides stability by ensuring planets b and d never simultaneously approach a true anomaly ν ∼ 300°, where their orbital paths allow very small separations. For configurations outside the stable a–ω space, the resonance becomes destructive. Once again, the stable region extends throughout the 1σ uncertainties on the best-fit orbit of planet d, with the location of the best-fit orbit lying close to the region of greatest stability.

Here again, our solution suffers from a lack of information regarding planet c. Fortunately, at P = 35 days and Msin i = 17 M, plausibility arguments suffice to rule out destabilizing interactions due to this inner planet. As demonstrated in Horner et al. (2011), planets tend to be stable when separated by ∼5 Hill radii. With planets b and c separated by nearly 12 Hill radii (measured from planet b), we expect little mutual influence. A similar argument is presented for KOI 961 (Muirhead et al. 2012). Nevertheless, we have performed a small number of simulations in which we include a planet with M = 17 M, a = 0.21 AU, and e = 0.17 (Mayor et al. 2011) to the stable configurations nearest to our best-fit orbital solution. In all cases, the stability of the system is unaffected; while the exact values of a and e for the giant planets are slightly different at each time step when planet c is included, the periods over which the orbital parameters vary remain unchanged, and the long-term evolution is the same regardless of whether c is included. We therefore conclude that excluding planet c from our larger analysis does not significantly affect our results.

Taken in concert, the results shown in Figure 5 reveal that dynamically stable coplanar solutions for the orbit of planet d require that it be trapped in mutual mean-motion resonance with planet b. Given that the nominal best-fit orbit lies perfectly within the region spanned by that resonance, and that a significant fraction of the 1σ error ellipse for planet d is dynamically stable in both ae and a–ω space, we find that our dynamical results are broadly in support of the existence of planet d, and may even be used to more tightly constrain its orbit.

7. DISCUSSION

With the addition of planet d, HD 204313 joins the growing list of stars hosting multiple gas giant planets. At G5 V spectral type and at virtually equal mass to the Sun, having two planets with masses ∼ 2–4 times that of Jupiter makes HD 204313 somewhat of an outlier on the correlation between stellar mass and giant planet fraction/mass (Johnson et al. 2011a). However, we confirm the star's super-solar metallicity measured by Ségransan et al. (2010), thereby offsetting the slight discrepancy with stellar mass.

As of 2012 May, there are 12 exoplanet systems in the exoplanets.org (Wright et al. 2011) database which contain giant planets believed to be in low-order resonances. However, only HD 45364 (3:2; Correia et al. 2009; Rein et al. 2010) and HD 200964 (4:3; Johnson et al. 2011b) host multiple gas giants in mean-motion resonances closer than 2:1. HD 204313 therefore joins a very small subset of the known planet systems. Furthermore, the minimum masses of planets b and d are considerably higher than either the HD 45364 or HD 200964 planets, making their continued stability even more remarkable. Interestingly, in addition to having the 3:2 resonance in common, HD 45364b/c and HD 204313b/d both have mass ratios close to the ∼3:1 Jupiter/Saturn mass ratio. These systems are thus valuable as a comparison to the Nice model (Tsiganis et al. 2005) for the formation of the outer solar system. In particular, Batygin et al. (2012) invoke a 3:2 resonance between Jupiter and Saturn in simulations which successfully reproduce the orbital configurations of the four outer solar system planets and the Kuiper Belt.

The importance of 3:2 mean-motion resonances within the solar system extends beyond the possible interactions between Jupiter and Saturn during the system's early evolution. Beyond the orbit of Neptune lie the Plutinos (named after the dwarf planet (134340) Pluto, the first known member). These objects, of which several hundred are currently known, are trapped within the 3:2 Neptunian mean-motion resonance. In Figure 6, we plot the Plutino distribution in ae space. The Plutino population contains objects with a wide range of eccentricities and inclinations, with the most eccentric objects crossing the orbit of Neptune, and some moving on orbits that can range as close as halfway between the orbits of Uranus and Neptune. The inclinations of the Plutinos range from 0 deg to over 30 deg. This wide distribution of orbital elements has been used to decipher the migration history of Neptune—the idea being that, as that giant planet migrated outward, objects were captured into the 3:2 mean-motion resonance and swept along with the planet, their orbits becoming ever more excited as they were carried along (e.g., Malhotra 1995). As a result, the Plutinos nicely map the extent of the stable region of the 3:2 mean-motion resonance with Neptune. Though the stable Plutinos do not range to quite as extreme eccentricities as are supported for the orbit of HD 204313d (as a result of the influence of Uranus on the evolution of the most eccentric members), it is striking that the region of stability occupied by the Plutinos is very similar to that obtained by our dynamical integrations, as can be seen when comparing Figures 5(a) and 6. We note in passing that the Hilda family of main belt asteroids are trapped in 2:3 mean motion resonance with Jupiter, orbiting with periods of ∼8 years. Despite their sometimes high orbital eccentricities (again up to, and in excess of, 0.3), these objects are protected from close encounters with the massive planet by the mean-motion resonance they occupy.

Figure 6.

Figure 6. Distribution in ae space of the population of trans-Neptunian objects between 35 and 43 AU. Note the concentration of objects just beyond 39 AU—the Plutinos, trapped in 3:2 mean-motion resonance with Neptune.

Standard image High-resolution image

Unlike our recent results for the 2:1 mean-motion resonance in the HD 155358 system (Robertson et al. 2012), our dynamical simulations for HD 204313 do not permit coplanar orbits outside the 3:2 resonance. Evidently, such a small period ratio is only stable when protected by the resonance. Furthermore, the existence of unstable orbits within the permitted parameter space shows that location within a resonance is not a guarantee of stability. Finally, it is interesting that while for both HD 155358 and HU Aquarii (Horner et al. 2011; Wittenmyer et al. 2012), the inner cutoff for stability in ae space is equal to the inner planet's apastron distance plus five Hill radii, our stable solutions for HD 204313 allow planet d to have values of a much smaller than this value. Once again, such a conclusion is supported by our knowledge of our solar system, in which resonant configurations ensure the stability of vast populations of objects on orbits that would otherwise be highly unstable. Prominent examples include the aforementioned Hildas (e.g., Franklin et al. 1993; Grav et al. 2012) and Plutinos (Malhotra 1995; Friedland 2001), along with the Jovian and Neptunian Trojans (Morbidelli et al. 2005; Sheppard & Trujillo 2006; Lykawka & Horner 2010; Horner et al. 2012a). In addition to the resonant exoplanets mentioned earlier, these populations reinforce the idea that such resonant scenarios are a common outcome from the planet formation process.

We note that the results of our stability analysis not only confirm the validity of our orbital model, but in fact place tighter constraints on the system's configuration than our fitting uncertainties alone. As long-term RV planet surveys such as ours become increasingly sensitive to systems with multiple long period companions, it is likely that additional systems with gas giants in close resonances will be discovered. Such dynamical simulations are therefore extremely valuable for understanding the true architecture of these systems, for which there may otherwise be considerable uncertainty as to the orbital parameters.

It is important to note that our claim that planets b and d are trapped in mutual 3:2 mean-motion resonance is based on simulations that assumed that their orbits are coplanar. However, as can be seen from the examples of the Hildas and Plutinos within our own solar system, resonant orbits can be dynamically stable for a wide range of mutual inclinations. It might instinctively seem that the coplanar case would actually be the least stable configuration, and therefore that mutually inclined orbits might allow a broader range of stable solutions. However, we note that in Horner et al. (2011), the authors considered a wide range of orbital inclinations in an attempt to address the apparent instability of the proposed HU Aquarii planetary system, and found that increasing the mutual inclination of the planets in question did little to remedy their instability. That said, it would certainly be interesting, in future work, to examine the influence of the mutual inclination of the orbits of the planets of the HD 204313 system. Fortunately, with predicted astrometric displacements of 0.432 mas and 0.264 mas for planets b and d, respectively, both planets should be accessible to astrometric measurements with the Hubble Space Telescope Fine Guidance Sensor (Nelan et al. 2010). Plus, the inclusion of inclination constraints would make HD 204313 a unique opportunity for comparison to the compact, multi-resonant planet systems discovered by the Kepler spacecraft (Holman et al. 2010; Lissauer et al. 2011; Cochran et al. 2011). If astrometry shows the HD 204313 planets to be coplanar, it would be strongly suggestive that similar migration mechanisms can result in systems as different as HD 204313 and the aforementioned Kepler planets. HD 204313 should therefore be considered a high-priority target for current and future astrometric surveys.

J.H. gratefully acknowledges the financial support of the Australian government through ARC Grant DP0774000. R.W. is supported by a UNSW Vice-Chancellor's Fellowship. M.E. and W.D.C. acknowledge support by the National Aeronautics and Space Administration under Grants NNX07AL70G and NNX09AB30G issued through the Origins of Solar Systems Program. A. E. Simon has been supported by the Hungarian OTKA Grants K76816, K83790, and MB08C 81013, the "Lendulet" Program of the Hungarian Academy of Science. This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org.

Footnotes

Please wait… references are loading.
10.1088/0004-637X/754/1/50