THE McDONALD OBSERVATORY PLANET SEARCH: NEW LONG-PERIOD GIANT PLANETS AND TWO INTERACTING JUPITERS IN THE HD 155358 SYSTEM

The velocities for HD 79498 and HD 197037 were obtained with the Smith Telescope's Tull Coudé Spectrograph (Tull et al. 1995) using a 18 slit, giving a resolving power of R = 60,000. Before starlight enters the spectrograph, it passes through an absorption cell containing iodine (I₂) vapor maintained at 50°C, resulting in a dense forest of molecular absorption lines over our stellar spectra from 5000 to 6400 Å. These absorption lines serve as a wavelength metric, allowing us to simultaneously model the instrumental profile (IP) and RV at the time of the observation. For each star, we have at least one high signal-to-noise ratio (S/N) iodine-free template spectrum, which we have deconvolved from the IP using the Maximum Entropy Method, and against which the shifts due to the star's velocity and the time-variant IP are modeled. Our reported RVs are measured relative to this stellar template and are further corrected to remove the velocity of the observatory around the solar system barycenter. All radial velocities have been extracted with our pipeline AUSTRAL (Endl et al. 2000), which handles the modeling of both the IP and stellar velocity shift.

Our RV data for HD 155358 and HD 220773 were taken with the High Resolution Spectrometer (HRS; Tull et al. 1998) on the queue-scheduled 9.2 m HET (Ramsey et al. 1998). As with the 2.7 m observations, the HET/HRS spectra are taken at R = 60,000 with an I₂ absorption cell. The fiber-fed HRS is located below the telescope in a temperature-controlled room. Separate stellar templates were obtained with HRS for these stars. Details of our HET observing procedure are given in Cochran et al. (2004).

While we use HRS spectra exclusively for obtaining RVs for HD 155358 and HD 220773, we have I₂-free spectra from the 2.7 m telescope for these stars as well. These spectra serve two purposes. First, they allow us to determine the stellar parameters for all five stars using the same instrumental setup. Also, the Tull 2.7 m coudé spectrograph provides Ca H and K indices, which contain information as to the activity levels of the stars.

Tables 1–4 list the relative velocities and their associated uncertainties for HD 79498, HD 155358, HD 197037, and HD 220773, respectively. Table 2 includes observations published in Cochran et al. (2007), but since all of our spectra have been re-reduced with our most up-to-date methods, the velocities presented here have a higher precision.

We determine the stellar parameters of our targets according to the procedure described in Brugamyer et al. (2011). The method relies on a grid of ATLAS9 model atmospheres (Kurucz 1993) in combination with the local thermodynamic (LTE) line analysis and spectral synthesis program MOOG⁵ (Sneden 1973). Using the measured equivalent widths of 53 neutral iron lines and 13 singly ionized iron lines, MOOG force-fits elemental abundances to match the measured equivalent widths according to built-in atomic line behavior. Stellar effective temperature is determined by removing any trends in equivalent widths versus excitation potential, assuming excitation equilibrium. Similarly, we compute the stellar microturbulent velocity ξ by eliminating trends with reduced equivalent width (≡W_λ/λ). Finally, by assuming ionization equilibrium, we constrain stellar surface gravity by forcing the abundances derived from Fe i and Fe ii lines to match.

We begin our stellar analysis by measuring a solar spectrum taken during daylight (through a solar port) with the same instrumental configuration used to observe our targets. The above procedure yields values of T_eff = 5755 ± 70 K, log g = 4.48 ± 0.09 dex, ξ = 1.07 ± 0.06 km s⁻¹, and log (Fe) = 7.53 ± 0.05 dex. We then repeat this analysis for our target stars, using the I₂-free template spectra. We note that our derived metallicities are, as is conventional, differential to solar. The high S/N and spectral resolution of our stellar templates allow us to make robust estimates of each star's effective temperature, log g, metallicity, and microturbulent velocity, which we include in Table 5. We also include photometry, parallax data, and spectral types from the ASCC-2.5 catalog (Version 3; Kharchenko & Roeser 2009), as well as age and mass estimates from Casagrande et al. (2011).

In addition to basic stellar parameters, we use the standard Ca H and K indices to evaluate the hypothesis that our observed RV signals are actually due to stellar activity. For HD 79498 and HD 197037 we have time-series measurements of the Mount Wilson S_HK index. For each of these stars, we obtain S_HK simultaneously with RV, so we include those values in the RV tables. From the average value of these measurements, we derive (via Noyes et al. 1984) log R'_HK, the ratio of Ca H and K emission to the integrated luminosity of the star, which we include in Table 5. The activity indicators are discussed in detail for each star below, but we note here that the sample overall appears to be very quiet, and we have no reason to suspect stellar variability as the cause of the observed signals.

As a second fail-safe against photospheric activity mimicking Keplerian motion, we calculate the bisector velocity spans (BVSs) for spectral lines outside the I₂ absorption region. As described in Brown et al. (2008), this calculation offers information regarding the shapes of the stellar absorption lines in our spectra. As photospheric activity that might influence RV measurements occurs, it alters the shapes of these lines, causing a corresponding shift in the BVSs. Our BVSs, then, offer a record of the activity of our target stars. Provided that our planetary signals are real, we expect the BVSs to be uncorrelated with our RV measurements. For each RV point, we have a corresponding BVS measurement computed from the average of the stellar lines outside the I₂ region.

To determine the orbital parameters of each planetary system, we first analyze each RV set using the fully generalized Lomb–Scargle periodogram (Zechmeister & Kürster 2009). We estimate the significance of each peak in the periodogram by assigning it a preliminary false-alarm probability (FAP) according to the method described in Sturrock & Scargle (2010). To conclude our periodogram analysis, we assign a final FAP for each planet using a bootstrap resampling technique. Our bootstrap method, which is analogous to the technique outlined in Kürster et al. (1997), retains the original time stamps from the RV data set and selects a velocity from the existing set (with replacement) for each observation. We then run the generalized Lomb–Scargle periodogram for the resampled data. After 10,000 such trials for each data set, the FAP is taken to be the percentage of resampled periodograms that produced higher power than the original RV set. Power in the Zechmeister & Kürster (2009) periodogram is given by Δχ² ≡ χ²₀ − χ_P², or the improvement of fit of an orbit at period P over a linear fit, so a "false positive" result in the bootstrap trial occurs when a random sampling of our measured RVs is better suited to a Keplerian orbit than the actual time series. In Table 6, we include both FAP estimates for comparison, but we adopt this bootstrap calculation as our formal confidence estimate.

Table 6. Derived Planet Properties

	HD 79498b	HD 155358		HD 197037ba	HD 220773b
		b	c
Period P (days)	1966.1 ± 41	194.3 ± 0.3	391.9 ± 1	1035.7 ± 13	3724.7 ± 463
Periastron passage T0	3210.9 ± 39	1224.8 ± 6	5345.4 ± 28	1353.1 ± 86	3866.4 ± 95
(BJD − 2,450,000)
RV amplitude K (m s−1)	26.0 ± 1	32.0 ± 2	24.9 ± 1	15.5 ± 1	20.0 ± 3
Mean anomaly M0b	329° ± 11°	129° ± 07	233° ± 09	186° ± 3°	226° ± 13°
Eccentricity e	0.59 ± 0.02	0.17 ± 0.03	0.16 ± 0.1	0.22 ± 0.07	0.51 ± 0.1
Longitude of periastron ω	221° ± 6°	143° ± 11°	180° ± 26°	298° ± 26°	259° ± 15°
Semimajor axis a (AU)	3.13 ± 0.08	0.64 ± 0.01	1.02 ± 0.02	2.07 ± 0.05	4.94 ± 0.2
Minimum mass Msin i (MJ)	1.34 ± 0.07	0.85 ± 0.05	0.82 ± 0.07	0.79 ± 0.05	1.45 ± 0.3
rms (m s−1)	5.13	6.14		8.00	6.57
Stellar "jitter" (m s−1)	2.76	2.49		2.02	5.10
FAP (from periodogram)	4.09 × 10−9	1.15 × 10−11	<10−14	8.53 × 10−8	⋅⋅⋅
FAP (from bootstrap)	<10−4	<10−4	<10−4	<10−4	⋅⋅⋅

Notes. ^aParameters include subtraction of linear RV trend in residuals. ^bEvaluated at the time of the first RV point reported for each system.

Download table as:  ASCIITypeset image

The periods identified in the periodograms are then used as initial estimates for Keplerian orbital fits, which we perform with the GaussFit program (Jefferys et al. 1988). We list the orbital parameters of each planet in Table 6 and describe the individual systems in detail below. As a consistency check, we also fit each orbit using the SYSTEMIC console (Meschiari et al. 2009), finding good agreement between our results in every case. An additional advantage of the SYSTEMIC console is that it computes a stellar "jitter" term—a measure of random fluctuations in the stellar photosphere—(e.g., Queloz et al. 2001) for each star, which we include in Table 6.

4.1.1. RV Data and Orbit Modeling

Our RV data for HD 79498 consist of 65 observations taken over 7 years from the 2.7 m telescope. The data have an rms scatter of 18.3 m s⁻¹ and a mean error of 4.15 m s⁻¹, indicating significant Doppler motion. Our RVs are listed in Table 1 and plotted in Figure 1.

Zoom In Zoom Out Reset image size

The periodogram for the HD 79498 RVs reveals a large peak around 1815 days, with a preliminary FAP of 4.09 × 10⁻⁹. Performing a Keplerian fit with P = 1815 days as an initial guess, we find a single-planet solution with the parameters P = 1966 days, e = 0.59, and K = 26.0 m s⁻¹, indicating a planetary companion with Msin i = 1.34M_J at a = 3.13 AU. Considering the high eccentricity and the width of the periodogram peak, these values are in good agreement with our initial guess. This fit produces a reduced χ² = 1.77 with an rms scatter of 5.13 m s⁻¹ around the fit. We note that we have not included the stellar "jitter" term in our error analyses, but we have verified that most of the χ² excess above unity can be attributed to our fitted value of 2.76 m s⁻¹. In a series of 10⁴ trials of our bootstrap FAP analysis, we did not find an improvement in Δχ² for any resampled data set, resulting in an upper limit of 10⁻⁴ for the FAP of planet b. Figure 1 shows our fit, plotted over the time-series RV data. We include the full parameter set for HD 79498b in Table 6.

Our analysis of HD 79498b admits a second, slightly different orbital solution of nearly equal significance. This solution converges at P = 2114 days, e = 0.61, and K = 26.2 m s⁻¹. The corresponding planetary parameters then become Msin i = 1.35M_J and a = 3.27 AU. For comparison, we include the plot of this fit in Figure 1. Although the shorter-period solution may appear to be driven mainly by the first data point, we note that our fitting routines converge to the 1966 day period regardless of whether that point is included and that the fitting statistics are still better for the shorter period with the point excluded. Clearly, the qualitative properties of the planet remain unchanged regardless of the choice of parameters. However, because of the slightly better fitting statistics (χ² = 1.77 versus χ² = 1.84, rms = 5.13 m s⁻¹ versus rms = 5.50 m s⁻¹) for the first model, we adopt it as our formal solution.

We computed the periodogram of the residuals around the fit of HD 79498b to search for additional signals. We see no evidence of additional planets in the system. We note that no additional signals appear in the residuals of the alternate fit for planet b either.

4.1.2. Stellar Activity and Line Bisector Analysis

4.1.3. Stellar Companions

4.2.1. RV Data and Orbit Modeling

We have obtained 113 RV points for HD 197037 over 10 years from the 2.7 m telescope. These data have an rms scatter of 13.1 m s⁻¹ and a mean error of 7.65 m s⁻¹. We report our measured RVs in Table 3 and plot them as a time series in Figure 2.

Zoom In Zoom Out Reset image size

Our periodogram analysis of HD 197037, shown in Figure 3(a), indicates a strong peak around P ∼ 1030 days. The power in this peak corresponds to a preliminary FAP of 8.53 × 10⁻⁸. Our one-planet Keplerian model yields parameters of P = 1036 days, e = 0.22, and K = 15.5 m s⁻¹, showing excellent agreement with the prediction of our periodogram. The inferred planet has a minimum mass of 0.79M_J and lies at an orbital separation of 2.07 AU. We overplot our model with the RVs in Figure 2.

Zoom In Zoom Out Reset image size

The addition of a single Keplerian orbit only reduces the residual rms of our RVs to 9.18 m s⁻¹, which is nearly a factor of two worse than our other one-planet fits. Furthermore, the periodogram of the residual RVs (Figure 3(a)) shows a significant increase in power at long periods. While any additional periods are too far outside our observational time baseline to properly evaluate, we can generate a preliminary two-planet fit with a second ∼0.7M_J planet with a period around 4400 days and an eccentricity 0.42. However, our current RV set can also be modeled as a single planet and a linear trend with slope −1.87 ± 0.3 m s⁻¹ yr⁻¹. Both fits give a χ² of 1.10 and an rms scatter of 8.00 m s⁻¹, so we adopt the more conservative planet-plus-slope model pending further observations. HD 197037 has no known common proper-motion companions within 30 arcsec, so it is very possible that we are seeing evidence for a distant giant planet or brown dwarf companion. The final orbital parameters reported in Table 6 are derived from the planet+slope model. We note that we see no additional signals in the residuals of either the two-planet or planet-slope fits.

4.2.2. Stellar Activity and Line Bisector Analysis

Even after accounting for the linear trend in the residuals around our fit to planet b, the scatter in our RVs is still higher than we typically expect from our 2.7 m data. Because HD 197037 is an earlier type (F7) than the other stars discussed in this paper, we expect a lower precision as a result of fewer spectral lines to determine velocities. Additionally, with log R'_HK = −4.53, it is the most active of the stars presented here. With an rms scatter in the BVS of 13 m s⁻¹, we must take particular care to ensure that we have not mistaken an activity cycle for a planetary signal. In Figure 4, we show our measured BVS and S_HK indices, which we use to evaluate the influence of stellar photospheric activity on the signal of HD 197037b. Using both a Pearson correlation test and least-squares fitting, we find no correlation to the RVs for the BVSs. While the S_HK indices appear somewhat correlated with the RVs to the unaided eye, the Pearson correlation coefficient of −0.26 indicates that the relation is not statistically significant. Periodogram analysis (Figure 3(b)) shows no power at the 1030 day peak in the BVS or S_HK time series. We do see a modestly (FAP = 0.01) significant peak at 19.1 days in the periodogram of our S_HK measurements, which we speculate may be the stellar rotation period. Although we are admittedly unable to completely rule out the possibility of stellar activity as the source of our observed RV signal, the lack of correlation between our activity indicators and the velocities and the absence of periodicity in S_HK and BVS around the fitted period lead us to the conclusion that a planetary orbit is the most likely cause.

Zoom In Zoom Out Reset image size

While the stellar activity measurements reinforce the planetary nature of the primary RV signal of planet b, the analysis of the long-term trend is not so clear. Our periodograms of the BVS and S_HK series both show an increase in power at very long periods, matching the behavior of the trend in the residual RVs, albeit at much lower power. Furthermore, we see a correlation between S_HK and the residual RVs around the one-planet fit with a Pearson correlation coefficient of −0.33, which is significant for the size of our data set. The residual RVs are uncorrelated with the BVSs, though, so the trends may be coincidental. HD 197037 will need continued monitoring to determine the true nature of this long-period signal.

4.3.1. RV Data and Orbit Modeling

Our RV data set (Table 4) for HD 220773 consists of 43 HET/HRS spectra taken over 9 years between 2002 July and 2011 August. The data have an rms of 11.7 m s⁻¹ and a mean error of 4.11 m s⁻¹.

Based on our data, we find evidence for a highly eccentric giant planet on a long-period orbit. Because of its high eccentricity, the period of planet b does not appear at significant power in our periodogram analysis. However, the high rms of our RV set, combined with the characteristic eccentric turnaround of the velocity time series (Figure 5), strongly indicates the presence of a planetary companion. In order to offset the lack of information from the periodogram, we have run GaussFit with a broad range of periods (2500–4500 days) and eccentricities (0.4–0.7). All fits converge unambiguously to a period of 3725 days with an eccentricity of 0.51, indicating a 1.45M_J planet at 4.94 AU. We include the plot of this model in Figure 5. The final fit gives a reduced χ² of 3.14 and an rms scatter of 6.57 m s⁻¹. Because the reduced χ² is so high, we have computed the model again after adding the 5.10 m s⁻¹ "jitter" term in quadrature to the errors listed in Table 4. While our fitted parameters and uncertainties do not change significantly from the values given in Table 6, the reduced χ² drops to 1.18, lending additional confidence to our solution. We do not currently see any evidence for additional signals in the residual RVs.

Zoom In Zoom Out Reset image size

To confirm the orbital parameters listed in Table 6, we use a genetic algorithm to explore the parameter space and evaluate the likelihood of the null hypothesis. The algorithm fits a grid of parameters (number of planets, masses, periods, eccentricities) to the RV set, allowing for a thorough exploration of how χ² behaves in response. Areas of parameter space that do not match the data are iteratively rejected, allowing the routine to converge on an optimal solution. We have performed 10,000 iterations of the algorithm with our RV data, considering periods between 3000 and 10,000 days. For this experiment, we have again added the stellar "jitter" term to our measurement errors. As shown in Figure 6, the genetic algorithm reaches χ² < 1.20 for a one-planet model with periods close to our fitted value in Table 6. This is a dramatic improvement over a zero-planet model, which yields χ² = 3.65.

Zoom In Zoom Out Reset image size

4.3.2. Stellar Activity and Line Bisector Analysis

Since the publication of Cochran et al. (2007), we have obtained 51 additional RV points for HD 155358 from HET. We have re-reduced all of our spectra with the latest version of AUSTRAL for consistency and include all 122 points in Table 2. The entire set of velocities, shown as a time series in Figure 7(a), covers approximately 10 years from 2001 June until 2011 August. These RVs have an rms scatter of 30.2 m s⁻¹, much higher than we expect from a mean error of just 5.38 m s⁻¹.

Zoom In Zoom Out Reset image size

We include the periodogram and window function for HD 155358 in Figure 8(a). As expected, we see a significant peak around 196 days, with a preliminary FAP of 1.15 × 10⁻¹¹. Fitting a one-planet orbit from this peak, we find an RV amplitude K = 32 m s⁻¹ with a period P = 195 days and an eccentricity of e = 0.23. Based on this fit, HD 155358b has a minimum mass of 0.85 Jupiter masses at a = 0.64 AU.

Zoom In Zoom Out Reset image size

Cochran et al. (2007) note that the periodogram for the residuals around the one-planet fit to planet b initially showed power around 530 days and 330 days. Comparing fits at both periods produced better results at the longer period for our data at that time. However, the periodogram of the residuals around planet b for our current data (Figure 8(a)) clearly indicates the shorter period as the true signal. The peak at 391 days has a preliminary FAP too low for the precision of our code (approximately 10⁻¹⁴), and we no longer see additional peaks at longer periods. While the window function for our sampling does show some power at the one-year alias, the period of this periodogram peak is sufficiently separated from the yearly alias (Figure 8(b)) that we are confident the signal is not due to our sampling.

Using 391 days as a preliminary guess for the period of planet c,  we performed a two-planet fit to our RVs. This updated fit changes the parameters of planet b to P = 194.3 days, Msin i = 0.85M_J, and e = 0.17, with the orbital separation remaining at 0.64 AU. Planet c then converges to a period of 391.9 days, with e = 0.16 and K = 25 m s⁻¹. The derived properties for planet c then become M = 0.82M_J and a = 1.02 AU. We plot this orbit over our RV data in Figure 7(a). The addition of planet c reduces our rms to 6.14 m s⁻¹ with a reduced χ² of 1.41. As a consistency check, we test a fit with parameters for planet c more closely matching those from Cochran et al. (2007) but find no satisfactory solution at this longer period. Holding the parameters for planet c fixed at the values derived in that study results in a reduced χ² of 11.2 and a residual rms of 15.4 m s⁻¹.

The BVSs for HD 155358 display an rms scatter of 25 m s⁻¹ and are uncorrelated with both the RVs and the residuals to the two-planet fit. Additionally, we see no significant peaks in the periodogram of the BVS time series. Also, while we present HET velocities here, we do have some 2.7 m spectra from a preliminary investigation of HD 155358, allowing us to examine the Ca H and K indices. S_HK shows very little activity, and we derive a log R'_HK of −4.54. The hypothesis that the two large-amplitude signals observed here are due to stellar activity can therefore conclusively be ruled out.

Even at the separation claimed in Cochran et al. (2007), planets b and c are close enough to interact gravitationally. Our updated orbital model now indicates they are actually much closer together. While their periods suggest that the planets are in a 2:1 mean-motion resonance (MMR), our preliminary orbital simulations showed the system to be unstable for a range of input values. We therefore decided to perform a highly detailed dynamical study of the system to investigate whether the orbits that best fit the data are indeed dynamically feasible. To do this, we performed over 100,000 unique simulations of the HD 155358 system using the Hybrid integrator within the n-body dynamics package MERCURY (Chambers 1999).

To systematically address the stability of the HD 155358 system as a function of the orbits of planets b and c, we followed Horner et al. (2011) and Marshall et al. (2010) and examined test systems in which the initial orbit of the planet with the most tightly constrained orbital parameters (in this case planet b) was held fixed at the nominal best-fit values. The initial orbit of the outermost planet 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 suggested to move on unusual orbits (e.g., Wittenmyer et al. 2012) and allow the construction of detailed dynamical maps for the system in orbital element phase space.

Keeping the initial orbit of the innermost planet fixed, we examined 31 unique values of semimajor axis for planet c, ranging from 0.96 AU to 1.08 AU, inclusive, in even steps. For each of these 31 initial semimajor axes, we studied 31 values of orbital eccentricity, ranging from the smallest value possible (0.0) to a maximum of 0.46 (corresponding to the best-fit value, 0.16, plus 3σ). For each of the 961 a–e pairs, we considered 11 values of initial mean anomaly and initial longitude of periastron (ω), resulting in a total suite of 116,281 (31 × 31 × 11 × 11) plausible architectures for the HD 155358 system.

In each simulation, the planet masses were set to their minimum (Msin i) values; the mass of the innermost planet was set to 0.85M_J and that of the outermost was set to 0.82M_J. As such, the dynamical stability maps obtained show the maximum stability possible (since increasing the masses of the planets would clearly increase the speed at which the system would destabilize for any nominally unstable architecture). 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.

The results of our simulations are shown in Figures 9 and 10. We present the multi-dimensional grid of orbital parameters over which we ran our simulations as two-dimensional cross sections, indicating the mean and median stability lifetimes over each of the runs at each grid point. So, for example, each grid point in our plots in a–e space reveals either the mean or the median of 121 unique simulations of an orbit with those particular a–e (i.e., 11 in ω times 11 in mean anomaly). Similarly, each point in the a–ω plots corresponds to the mean (or median) of 341 separate trials (31 in e times 11 in ω).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The first thing that is apparent from Figure 9(a), which shows the mean lifetime of the system as a function of semimajor axis and eccentricity, is that the stability of the system is a strong function of eccentricity. Solutions with low orbital eccentricity are typically far more dynamically stable than those with high eccentricity. In addition, the stabilizing effect of the 2:1 MMR between planets b and c can be clearly seen as offering a region of some stability to even high eccentricities between a ∼ 1.05 and 1.25 AU. However, it is apparent when one examines the median lifetime plot (Figure 9(b)) that a significant fraction of eccentric orbits in that region are dynamically unstable (the reason for the apparently low median lifetimes in that region). The reason for this is that the dynamical stability of orbits in this region, particularly for high eccentricities, is also a strong function of the initial longitude of periastron for planet c's orbit. Indeed, comparison of this figure to those showing the influence of the longitude of periastron for planet c reveals that the most stable regions therein lie beyond the 1σ error bars for the nominal orbit. For low eccentricities (e < 0.1), in the vicinity of the 2:1 resonance, the orbit is stable regardless of the initial longitude of periastron, but for high eccentricities (e > 0.15), the stable regions lie toward the edge of the allowed parameter space, making such a solution seem relatively improbable.

It is apparent from close examination of Figure 9(b) (which shows the median stability as a function of semimajor axis and eccentricity) that the key determinant of the stability of the system (particularly for non-resonant orbits) is actually the periastron distance of planet c. The sculpted shape of the stability plot, outward of a ∼ 0.98 AU, is very similar to that observed for the proposed planets around the cataclysmic variable HU Aquarii (Horner et al. 2011; Wittenmyer et al. 2012). As was found in that work, the dividing line between unstable and stable orbits again seems to fall approximately five Hill radii beyond the orbit of the innermost planet. Any orbit for planet c that approaches the orbit of b more closely than this distance will be unstable on astronomically short timescales (aside from those protected from close encounters by the effects of the 2:1 MMR between the planets). On the other hand, orbits that keep the two planets sufficiently far apart tend to be dynamically stable. We note, too, that the sharp inner cutoff to this broad region of stability, at around 0.99 AU, corresponds to the apastron distance of the innermost planet plus five Hill radii. Once again, orbital solutions that allow the two planets to approach closer than five times the Hill radius of the inner planet destabilize on relatively short timescales.

The results of our dynamical simulations show that a large range of dynamically stable solutions exist within the 1σ errors on the orbit of HD 155358c. Given the breadth of the resonant feature apparent in Figure 9(a) (the mean lifetime as a function of a and e), it seems clear that all orbits within 1σ of the nominal best fit will be strongly influenced by that broad resonance. As such, it seems fair to conclude that these two planets are most likely trapped within their mutual 2:1 MMR, although it is also apparent that non-resonant solutions also exist that satisfy both the dynamical and observational constraints. The stability restrictions are consistent with our orbital fit, but the parameters listed in Table 6 have not been modified to include the information obtained from the orbital simulations.

Using the [Fe/H] and T_eff derived herein, an [α/Fe] of 0.32 (Fuhrmann & Bernkopf 2008), and an age of 10.7 Gyr (Casagrande et al. 2011), we have fit Yonsei-Yale isochrones (Demarque et al. 2004) to HD 155358. With mass estimates ranging from 0.87 M_☉ (Lambert et al. 1991) to 0.92 M_☉ (Casagrande et al. 2011), HD 155358 has a luminosity somewhere between L_* = 1.14 L_☉ and L_* = 1.67 L_☉. If we assume that the location of the habitable zone scales as $\sqrt{L_{\ast }}$, then at an orbital separation of 1.06 AU, planet c lies within the habitable zone of its parent star (as defined by Kasting et al. 1993). Typically, Jovian planets are not considered to be potential habitats for Earth-like life (Lammer et al. 2009). However, it has been suggested that such planets could host potentially habitable satellites (e.g., Porter & Grundy 2011) or even that they could dynamically capture planet-sized objects during their migration as satellites or Trojan companions⁶ (e.g., Tinney et al. 2011). In fact, the presence of two giant planets on relatively close-in orbits may indirectly increase the water content on terrestrial satellites through radial mixing of planetesimals rich in ices (Mandell et al. 2007), a key requirement for such objects to be considered habitable (e.g., Horner & Jones 2010). Furthermore, any satellites sufficiently large to be considered habitable would also be subject to significant tidal heating from their host planet, which would likely act to increase their habitability when they lie toward the outer edge of Kasting's habitable zone (as they would have earlier in the life of the system). The induced tectonics would also potentially improve the habitability of any such moons (e.g., Horner & Jones 2010).

Our dynamical analysis suggests that the two planets in the HD 155358 system are most likely trapped in mutual resonance. It might initially seem that such orbits would be variable on the long term, such that they would experience sufficiently large excursions as to render any satellites or Trojan companions uninhabitable. However, we know from our own solar system that long-term resonant captures can be maintained on timescales comparable to the age of the solar system (e.g., the Neptune Trojans; Lykawka & Horner 2010). As such, it is reasonable to assume that, if the planets are truly trapped in mutual resonance, they could have been on their current orbits for at least enough time for any moons or Trojans of planet c to be considered habitable. However, the current best-fit eccentricity for planet c is sufficiently high that it lies at the outer limit of the range in e that would allow for a habitable exomoon (Tinney et al. 2011), which might limit the potential habitability of any exomoons in the system. That said, the 1σ range of allowed orbits extends to relatively low eccentricities (e = 0.06), which is certainly compatible with potential habitability. As such, it is possible that such moons, if present, could be potential habitats for life.