THE MAGELLAN PFS PLANET SEARCH PROGRAM: RADIAL VELOCITY AND STELLAR ABUNDANCE ANALYSES OF THE 360 au, METAL-POOR BINARY "TWINS" HD 133131A & B^*

In the study of exoplanets, how unique our solar system is in comparison to other planetary systems is one of the most compelling yet elusive questions. Planet-detecting transit surveys, especially the Kepler mission, as well as radial velocity (RV) surveys have revealed populations of planets unlike any in our solar system, of which the most numerous appear to be "super-Earths" that are larger than Earth but smaller than Uranus and Neptune (Borucki et al. 2011; Batalha et al. 2013), and that may or may not be gas-dominated (Ikoma & Hori 2012; Lopez & Fortney 2014; Marcy et al. 2014; Rogers 2015; Wolfgang & Lopez 2015). Both techniques for finding planets are biased toward short-period and large/massive planets, but the long baseline of RV observations (over 15 years in some cases) makes the RV technique more suitable for detecting long-period planets, like Jupiter and Saturn, than any current or planned transit surveys. Current theories of our solar system formation point to Jupiter and Saturn as important players in the dynamical shaping of the asteroid belt (Morbidelli et al. 2010; Walsh et al. 2011), and to their migratory dance toward and away from the Sun as the cause of a late delivery of volatile material to the inner rocky planets (Owen & Bar-Nun 1995; Morbidelli et al. 2000; Horner & Jones 2010). The architecture of our solar system is likely due in large part to Jupiter (and to a lesser extent, Saturn); thus, how unique our solar system is may depend on the occurrence of long-period giant planets.

Recently, the frequency of Jupiter analogs—defined as $5\lt { \mathcal P }\lt 20$ years, $0.3\lt { \mathcal M }\lt 3\,{{ \mathcal M }}_{{\rm{J}}}$, $e\lt 0.3$—was determined from a sample of over 1100 stars observed with Keck for RV variations due to planets (Rowan et al. 2016). These authors, correcting for relative observability around each of the stars, found that ∼3% of stars host Jupiter analogs, with their 10%–90% confidence intervals suggesting a frequency of 1%–4%. Wittenmyer et al. (2016) also recently published a Jupiter analog occurrence rate, corrected for incompleteness, based on 17 years of AAT observations of 202 solar-type stars, finding a Jupiter analog frequency of ${6.2}_{-1.6 \% }^{+2.8 \% }$. Their "Jupiter analog" definition included any planets with ${ \mathcal M }\gt 0.3\,{{ \mathcal M }}_{{\rm{J}}}$ and $a\gt 3\,\mathrm{au}$, and their confidence interval is only 65.7% around the peak of the posterior distribution function. Comparing "apples to apples," the two studies are consistent within errors.

These studies suggest that Jupiter analogs may indeed be rare. However, it is also interesting to consider long-period giant planets that fall outside this strict definition, particularly the e limit. Giant planets with higher eccentricities may have drastically influenced the formation/survival of any interior rocky planets (Thébault et al. 2002; Raymond 2006; Matsumura et al. 2013), or may be relics of a multi-planet system where one or more planets were ejected, potentially shedding light on the existence of free-floating planets (e.g., Rasio & Ford 1996; Weidenschilling & Marzari 1996). In fact, the working theory of giant planet evolution in the solar system proposes that an ice giant scattered off Jupiter to cause it to "jump" over a 2:1 mean motion resonance with Saturn, that this planet–planet scattering excited Jupiter's eccentricity (e.g., Tsiganis et al. 2005; Brasser et al. 2009; Morbidelli et al. 2010; Nesvorný 2011; Batygin et al. 2012), and that early-formed short-period planets would be scattered into the Sun (Batygin & Laughlin 2015). If long-period giant planets are more frequently eccentric, this would be an important clue to understanding our solar system formation and planet habitability in general (e.g., Kaib & Chambers 2016).

Here we report the detection of three giant planets on long, eccentric orbits from data collected primarily with the Planet Finder Spectrograph (PFS) on Clay/Magellan II at Las Campanas Observatory. In addition, we perform a precision stellar parameter and abundance analysis on the two host stars, which form a ∼360 au binary and are "twins" of identical spectral type. Interestingly, we find small but significant differences in the refractory element (with condensation temperatures $\gt 1000$ K) abundances between the two host stars. We discuss how these differences may be related to differences in planet formation/composition, and the nature of the host star binary.

HD133131 (HIP 73674), composed of a pair of bright (V = 8.40 and 8.42) G2V stars, was first recorded as a 7'' binary by Stock & Wroblewski (1972) on objective prism plates taken with the 24'' Curtis Schmidt Telescope at CTIO. The two stars are included in the Geneva–Copenhagen Survey (GCS), a comprehensive catalog providing kinematics and Galactic orbits, companion detections, and distances of nearly 17000 late-type stars in the solar neighborhood (Nordström et al. 2004; Holmberg et al. 2007, 2009). The latest large update to the catalog was by Casagrande et al. (2011), who used the infrared flux method (Casagrande et al. 2010) to derive new temperatures, metallicities (which are calibrated to high-resoultion spectroscopy), and ages. They report updated effective temperatures (T_eff) of 5791 and 5768 K and updated [Fe/H]⁶ values of −0.40 and −0.42 for HD 133131A and B, respectively, and report an age between 5.52 and 5.96 Gyr based on a Bayesian approach using BASTI (Pietrinferni et al. 2004, 2006, 2009) and Padova (Bertelli et al. 2008, 2009) isochrones and log(T_eff), absolute Johnson V magnitude, and metallicity. However, their age distributions are wide, and extend to ∼12.8 Gyr. The average distance reported in Nordström et al. (2004) and Holmberg et al. (2009) is 47 pc, adopted from the trigonometic distance from Hipparcos (original reduction from ESA 1997 and new reduction from van Leeuwen 2007). Tokovinin (2014a) included the HD 133131 system in his imaging survey of wide binaries, and found no evidence of a companion around HD 133131B between 0042 and 15.

Desidera et al. (2006a, 2006b) also included HD 133131A and B in their study of abundance differences between components of wide binaries, measured with the FEROS spectrograph at ESO La Silla. The separation between the two stars is 74 (from Hipparcos; ESA 1997), translating to a physical separation of 360 au. We calculate a period of ∼4240 years for HD 133131AB from the median separation (Washington Double Star Catalog⁷ ), parallax (van Leeuwen 2007), and assumed solar mass of both stars. Tokovinin (2014b), who compiled results from several of his own works, also reported a period of ∼4240 years for the HD 133131AB pair. This makes HD 133131A and B the smallest-separation "twin" binary system analyzed for high-precision abundance differences (the next closest is 16 Cyg AB, at ∼860 au; Eggenberger et al. 2003). The Desidera stellar properties of HD 133131A and B, listed in Table 6, were determined with a methodology very similar to that used here (described below), so represent the best comparison. Desidera et al. used their spectroscopic T_eff and abundances with isochrones of Girardi et al. (2002) and the bolometric corrections of Kurucz & Bell (1995) to iteratively derive masses of 0.95 ${{ \mathcal M }}_{\odot }$ and 0.93 ${{ \mathcal M }}_{\odot }$ for the A and B components, respectively, which are the stellar masses we assume in our RV search for planets below. They derive a much older age for the HD 133131 system, of ∼9.86 Gyr, from chromospheric activity measured in the Ca H & K lines (log $R^{\prime}$(HK)).

The RV observations of HD 133131A and B are part of the large Magellan Planet Search Program, which began in 2002 and is surveying a sample of ∼500 of the nearest stars (<100 pc). The survey was started with observations from the MIKE echelle spectrograph (Bernstein et al. 2003), mounted for a limited time on the Magellan I (Baade), but mostly on Magellan II (Clay), 6.5 m telescopes at Las Campanas Observatory. In 2010, the survey switched to using the Carnegie PFS (Crane et al. 2006, 2008, 2010), a temperature-controlled high-resolution echelle spectrograph built for precision RV observations, on Magellan II. Both instruments contain an iodine absorption cell (Marcy & Butler 1992) that imprints the reference iodine spectrum on the incoming starlight, providing a measurement of the instrument point-spread function (PSF) and a precise wavelength solution. Doppler shifts from the spectra are obtained using the technique described in Butler et al. (1996). In brief, the iodine region of the stellar spectrum (between ∼5000–6200 Å) is divided into 2 Å chunks, on which a forward modeling procedure is performed, providing an independent measurement of the instrument PSF, wavelength and Doppler shift. The measured velocity for each spectrum is calculated from the weighted mean of the independent chunk velocities; the reported internal uncertainty is the standard deviation of all of the chunk velocities measured from one spectrum. The weighted mean Doppler shift and the internal uncertainty for each observation are reported in Tables 1–3.

3.1. Magellan/MIKE RV Observations

3.2. Magellan/PFS Observations

Both HD 133131A and B were observed with PFS, the former observations ranging from 2010 February to 2015 September, and the latter from 2010 August to 2015 September. PFS has a more limited wavelength range than MIKE (3880–6680 Å), but still covers the entire iodine wavelength region, Ca ii H and K, and Hα. We use a 0.5 × 25 slit for target observations, providing R ∼ 80,000 in the iodine region. The total exposure times for the A component range from 285 to 720 s, and for the B component range from 282 to 800 s.

Similar to MIKE, PFS calibrations taken at the beginning of each night include 20–30 flat-field images, two iodine exposures, two rapidly rotating B star exposures, and one or two ThAr exposures. A modification of the MIKE pipeline is used for PFS raw frame reduction and spectral extraction.

Zoom In Zoom Out Reset image size

We use the SYSTEMIC console for Keplerian fitting of the RV data (Meschiari et al. 2009). First, we computed the error-weighted (w_j = 1/${\sigma }_{j}^{2}$), normalized Lomb–Scargle (L-S) periodogram (Gilliland & Baliunas 1987; Zechmeister & Kürster 2009) of each star's observations (MIKE+PFS for HD 133131A, PFS for HD 133131B), shown in Figure 2. The false alarm probabilities were computed by scrambling the data sets 10,000 times and sampling the periodogram at 80,000 frequencies, to calculate the probability that the power at each frequency could be exceeded by random chance. Overplotted in red on the periodograms is the spectral window, or periodogram of the sampling, which shows the expected peaks from the observing cadence, the sidereal and solar days, and the solar year (Dawson & Fabrycky 2010).

Zoom In Zoom Out Reset image size

To evaluate the stellar activity levels of HD 133131A and B, we measure Mt. Wilson S-index values in every RV observation spectrum; these values are listed in Tables 1–3. The S-index compares the flux in triangle-weighted bins with FWHMs of 1.09 Å centered on the Ca ii H&K lines (at 3968.47 and 3933.66 Å) to the flux in two rectangular 20 Å wide continuum regions centered on 3901 and 4001 Å (Duncan et al. 1991). This index is known to be correlated with spot activity on the stellar surface, and serves as a proxy for chromospheric activity that could cause RV shifts that mimic those induced by planets. Previous measurements of S-indices for stars in the Magellan Planet Search Program monitored with the MIKE spectrograph were detailed in Arriagada (2011). In Figure 1 we show how these values change with time in HD 133131A and B, specifically focusing on the PFS data, since we do not have MIKE data for HD 133131B. In this work we only measure the Ca ii H line and surrounding flux for PFS S-index values, effectively taking the ratio of the flux in the H line to the flux in the R continuum region, centered at 3996.5 Å as in Santos et al. (2000), referred to as S_COR in that work. The two large outliers in the S-index values for HD133131A, shown in the inset of Figure 1, are caused by small errors in the wavelength calibration, which is not well-constrained outside of the iodine region. The S-index is thus measured on the upward slope of the Ca ii H line instead of centered in the middle, causing an increase in its value. We include these outliers to avoid any "special" treatment of specific spectra, and the effect is small—the final median S-index value for HD 133131A changes by $\sim 3\times {10}^{-5}$ when we discount the two large outliers. We show in Figure 3 that the peaks in the L-S periodograms for the Mt. Wilson S-index measurements, which are a proxy for stellar magnetic activity modulation, are not as significant as, and do not correspond to, the peaks in the RV periodograms.

Zoom In Zoom Out Reset image size

The highest peaks in the L-S periodogram provide a first guess of the periods of planets around the two stars. For HD 133131A, the highest periodogram peak (that is not at an expected period due to time sampling) is at ∼660 days with a FAP of 1.28 × 10⁻⁸. For HD 133131B, the highest peak is at ∼6100 days with a false alarm probability (FAP) of $2.03\times {10}^{-2}$. The procedure that SYSTEMIC uses to model the RV signals is described in detail in Meschiari et al. (2009) and Vogt et al. (2015), but briefly, we begin by fitting a one-planet Keplerian model with six free parameters (period, mass, mean anomaly, eccentricity, longitude of pericenter, and a vertical offset to account for differences in the velocity zero point between data sets). Each RV measurement is represented by the predicted velocity, the formal (observing error) uncertainty, and an additional error term accounting for scatter about the fit (e.g., from underestimated measurement errors, stellar jitter, other astrophysical sources of variation) that is the same for each observation in the data set. The best-fit parameters are then derived by optimizing the log-likelihood of the model:

$$\begin{eqnarray}&&\mathrm{log}\,{ \mathcal L }=-\displaystyle \frac{1}{2}\left[{\chi }^{2}+\displaystyle \sum _{i=1}^{{N}_{o}}\mathrm{log}\,({e}_{{\rm{i}}}^{2}+{s}_{{\rm{i}}}^{2})+{N}_{o}\mathrm{log}\,(2\pi )\right]\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{{N}_{o}}{({V}_{i}-{v}_{i})}^{2}/({e}_{i}^{2}+{s}_{i}^{2}),\end{eqnarray} \tag{ 2 }$$

and V_i is the predicted velocity, v_i is the observed velocity, e_i is the formal error, and s_i is the additional error term.

Zoom In Zoom Out Reset image size

In a semi-automatic way, SYSTEMIC chooses the best parameters for a planet fit using a downhill simplex algorithm (AMOEBA; Nelder & Mead 1965; Press et al. 1992).⁸ The first planet fit is almost always the lowest FAP peak in the RV periodogram; this is the case for HD 133131A and B. After the first planet fit, the periodogram is re-calculated and the FAPs of the peaks reassessed. In the case of the B component, there are no remaining significant peaks (see Figure 5); stopping with a single planet fit, displayed in Figure 6, results in a final rms of 1.59 m s⁻¹. The detailed parameters of HD 133131Bb are given in Table 5, with the errors explained below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The next highest peak in the periodogram of HD 133131B after removal of planet b is at 5.88 days with a FAP of 4.63 × 10⁻¹, as shown in Figure 5. This is not high enough to merit a detection designation, and when including a planet at 5.88 days in our fit to the HD 133131B RV data, the reduced ${\chi }^{2}$ is significantly less than 1 (0.45), suggestive of over-fitting. However, with a best-fit ${ \mathcal M }\,\sin i$ of 0.018 [${{ \mathcal M }}_{{\rm{J}}}$], this potential planet would be sub-Neptune in mass, likely falling into the super-Earth planet regime, perhaps at the boundary between planets that are dominated by a volatile envelope and those that are not (e.g., Rogers et al. 2011; Lopez et al. 2012; Zeng & Sasselov 2013; Wolfgang & Lopez 2015). Only ∼15 other planets have both masses and periods less than or equal to the values for the potential HD 133131Bc planet. Adding a low-mass short-period planet to the system does not disrupt the stability—a 100,000 year BulirschStoer integration of the two planets' orbits shows no orbital or eccentricity overlap, and no significant changes in semimajor axis. Further implications of this potential planet are explored in Section 5; in any case HD 133131B warrants continued monitoring.

In the case of HD 133131A, after removing a one-planet fit, there is still significant power, with a FAP of $1.37\times {10}^{-4}$, at ∼3000 days as seen in Figure 7, left. We fit this as a second planet around the star, resulting in a periodogram with no remaining significant power (Figure 7, right). Thus we stop with a two-planet fit, displayed in Figure 8 (top), with a final rms of 9.38 m s⁻¹. The detailed parameters of HD 133131Ab and Ac are given in Table 4, with the errors explained below. We note that only including the 26 Magellan/PFS data points in the HD 133131A fits, and not the 17 Magellan/MIKE data points, decreases the rms to 1.82 m s⁻¹, and does not alter the fitted planet parameters outside their errors (see below), with the exception of the mean anomalies for both planets (which are basically reversed in the sans-MIKE data fits, and are not well constrained in either case).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Best-fit Solution for Planets b and c Around HD 133131A

Parameter	b			c, high e			c, low e
	Best-fit	Median	MAD	Best-fit	Median	MAD	Best-fit	Median	MAD
Period (days)	649	649	3	3407	3686	970	3925	3568	1084
${ \mathcal M }\,\sin i$ $({{ \mathcal M }}_{{\rm{J}}})$	1.43	1.42	0.04	0.48	0.44	0.14	0.63	0.42	0.15
Mean anomaly (deg)	268	267	11	56	79	82	106	101	87
Eccentricity	0.32	0.33	0.03	0.47	0.50	0.22	0.20[fixed]	0.49	0.22
ϖ (deg)	14	16	4.7	104	104	31	98	100	37
K (m s−1)	36.68	36.52	0.93	7.57	7.15	2.10	8.56	6.89	2.20
Semimajor axis (au)	1.44	1.44	0.005	4.36	4.59	0.82	4.79	4.49	0.92
Periastron passage time (JD)	2452326	2452327	23	2452277	2451701	1171	24523231	2452327	21
Stellar mass $({{ \mathcal M }}_{\odot })$	0.95
Reduced ${\chi }^{2}$	0.980						0.965
rms (m s−1)	9.38						9.39
Stellar jitter (m s−1)	13.7 (MIKE)
	1.92 (PFS)
Instrument rms (m s−1)	6.60 (MIKE)
Instrument rms (m s−1)	1.47 (PFS)
Data points	43
Span of observations (JD)	2452808.68 −
	2457268.49

Note. All elements are defined at epoch JD = 2452808.68. Uncertainties are reported in brackets.

Download table as:  ASCIITypeset image

Considering the low amplitude and incomplete phase coverage of this second planet around HD 133131A, we performed some additional tests as suggested by the referee to rule out stellar activity as the cause of the RV variation. First, we checked that excluding the two outliers in the S-index values measured from PFS (Figure 1) did not significantly alter the S-index periodogram. The heights of all of the peaks are slightly increased in the outlier-free periodogram (below, Figure 4), but none of the peaks is significant (no FAPs below 1). The peaks around 1000 days in both Figures 3 and 4 are at shorter periods (∼1280 days in Figure 3 and ∼940 days and ∼2350 days in Figure 4) than the ∼3500 day peak in the residuals of HD 133131A after one planet is removed (Figure 7).

Second, we checked the correlation between the S-index and the residuals in the RV data after fitting for a single planet (shown in Figure 7). We find no significant correlations—Spearman's correlation coefficient is −0.097, and the Pearson correlation coefficient is −0.15; this remains true after the removal of the same two S-index outliers as above (resulting in correlation coefficients of −0.045 and 0.11, respectively). These tests lend confidence to the planetary nature of the residual signal in the RV measurements of HD 133131A after fitting for one planet (b).

As with HD 133131B, we integrated the orbits of the b and c planets around HD 133131A over 100,000 years with the BulirschStoer extrapolation algorithm in SYSTEMIC. In this case, the fits decribed in columns 5–7 of Table 4 result in an unstable configuration, defined as a >10% change in semimajor axis in planet c (see Figure 9, top). However, if we manually reduce and fix the eccentricity of HD 133131Ac at 0.20, instead of the median value of ∼0.5, a similar orbit integration results in a more stable configuration that does not include eccentricity overlap between the two planets (Figure 9, bottom). We report the low-eccentricity fit parameters for planet c in the last three columns of Table 4; the parameters for planet b are unchanged. This fit is shown in the bottom panel of Figure 8. Though there is almost total agreement between the parameters for planet HD 133131Ac in the high- and low-eccentricity fits when errors are included, and there is no change within the errors to the parameters of planet HD133131Ab, we favor the low-eccentricity fit for planet c based on its greater stability.

4.1. Error Analysis

To characterize the marginal distribution of the parameters of the model above, we used the functionality of SYSTEMIC to implement the Markov chain Monte Carlo algorithm (MCMC; e.g., Ford 2005, 2006; Gregory 2011) fitting paired with flat priors on log P, log ${ \mathcal M }$, and the other parameters. In previous works using SYSTEMIC (e.g., Rowan et al. 2016: Vogt et al. 2015), the noise parameter s_j was fit with a modified Jeffrey function as a prior (see Vogt et al. 2015, Section 4); the MCMC routine then also returns a best-fit s_j "jitter" term and its marginal distribution.

Instead of relying solely on the MCMC resulting best-fit s_j values, here we outline an analytic approach that can be applied to any RV data set to estimate the stellar jitter term without including them as free parameters. In the case of HD 133131A, the Magellan/MIKE and Magellan/PFS data have significantly different errors (with MIKE errors being on average ∼4× as large as PFS errors), meaning that the s_j stellar jitter values for the two instruments will, in reality, be different. We estimated the s_j jitter term separately for each data set with the following procedure, and added it in quadrature to the formal error of each RV data point to get new errors, which we used as the initial input into SYSTEMIC.

The purpose of the s_j is to compensate for any noise, beyond the formal errors, that separates the data from the final best fit, as defined by a reduced ${\chi }^{2}=1$. The contribution to the ${\chi }^{2}$ of each data point is ${\sigma }_{\mathrm{fit}}^{2}{/\sigma }_{\mathrm{inst}}^{2}$, where ${\sigma }_{\mathrm{fit}}$ is the residual on the fit, and ${\sigma }_{\mathrm{inst}}$ is the instrumental (+jitter) error. To find the appropriate s_j, we wanted the sum of these values, divided by the number of data points minus the number fitted parameters, to equal 1:

$$\begin{eqnarray}&&\left[\displaystyle \sum _{i}^{{N}_{o}}\displaystyle \frac{{\sigma }_{\mathrm{fit},i}^{2}}{({\sigma }_{\mathrm{inst},i}^{2}+{s}_{j}^{2})}\right]\times \displaystyle \frac{1}{{N}_{o}-{N}_{\mathrm{param}}}=1.\end{eqnarray} \tag{ 3 }$$

We separately fit the 26 PFS and 17 MIKE RV data points to find the residuals, for the ${\sigma }_{\mathrm{fit},i}$ values. To decide how ${N}_{\mathrm{param}}=12$ should be divided between the PFS and MIKE data points, we estimated the effective weight for each data set as the square of the ratio of the rms values of each data set, (1.47 m s⁻¹/6.60 m s⁻¹)${}^{2}\,=\,$0.05, so that that number of free parameters allocated to the MIKE data are 0.05 × 12 ∼ 1, and PFS is allocated 11 free parameters. The final s_j value that solved the equation above for the PFS data was 1.92 m s⁻¹, and for the MIKE data it was 13.70 m s⁻¹. These values are reported in Table 4.

Even though the HD 133131B data are only from PFS, we can apply the same procedure to keep the analysis consistent, this time allocating all six parameters to the PFS data. We derive a s_j value of 1.03 m s⁻¹ for HD 133131B, as reported in Table 5.

Zoom In Zoom Out Reset image size

Table 5.  Best-fit One-planet Solution for HD 133131B

Parameter	b
	Best-fit	Median	MAD
Period (days)	6119	5769	831
${ \mathcal M }\,\sin \,i$ $({{ \mathcal M }}_{{\rm{J}}})$	2.50	2.50	0.05
Mean anomaly (degree)	302	299	6
Eccentricity	0.62	0.61	0.04
ϖ (degree)	103	103	3
K (m s−1)	37.29	37.41	0.65
Semimajor axis (au)	6.40	6.15	0.59
Periastron passage time (JD)	2450298	2450644	828
Stellar mass $({{ \mathcal M }}_{\odot })$	0.93
Reduced ${\chi }^{2}$	0.998
rms (m s−1)	1.59
Stellar jitter (m s−1)	1.03
Instirument rms (m s−1)	1.49
Data points	25
Span of observations (JD)	2455428.51− 2457268.50

Note. All elements are defined at epoch JD = 2455428.51. Uncertainties are reported in brackets.

Download table as:  ASCIITypeset image

Armed with appropriate total errors for each RV data point, we ran the SYSTEMIC MCMC algorithm with two chains, skipping the first 1000 samples (burn in), until the potential scale reduction factor R was equal to 1.1 (Gelman & Rubin 1992; Brooks & Gelman 1997), which we took to indicate convergence. The resulting marginal distributions of the orbital elements are shown in Figures 10, 11 (high-eccentricity fit for planet c), 12 (low-eccentricity fit for planet c), and 13, with the best-fit values shown with a red dot. We summarize these distributions by reporting the best fit (highest likelihood), median, and median absolute deviation (MAD) for each parameter in Tables 4 and 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We confirmed our simple, analytic estimation of s_j values for each data set by running the same MCMC algorithm in SYSTEMIC on input data without the s_j term added in quadrature to the formal errors. The marginal posterior distributions for the planet parameters were all consistent with those derived above, and the distributions of the s_j values were peaked around the values we estimated. The median±MAD s_j value for HD 133131A PFS is 1.66 ± 0.63, and for MIKE it is 13.8 ± 3.15; for HD 133131B PFS it is 0.63 ± 0.58.

4.2. Planet Detectability

Given the interesting nature of the potential planet HD 133131Bc at 5.88 days, along with the stellar abundance differences described in Section 5.2, we explored what limits our current data could place on additional planet detection around HD 133131A and B. The RV semi-amplitude, $K=({\upsilon }_{r,\max }-{\upsilon }_{r,\min })/2$, can be expressed in terms of the interesting quantities ${{ \mathcal M }}_{P}$ and P,

$$\begin{eqnarray}&&{K}_{* }=\displaystyle \frac{8.95\,\mathrm{cm}\,{{\rm{s}}}^{-1}}{\sqrt{1-{e}^{2}}}\displaystyle \frac{{{ \mathcal M }}_{P}\,\sin \,i}{{M}_{\oplus }}{\left(\displaystyle \frac{{{ \mathcal M }}_{* }+{{ \mathcal M }}_{P}}{{{ \mathcal M }}_{\odot }}\right)}^{-2/3}{\left(\displaystyle \frac{P}{\mathrm{year}}\right)}^{-1/3}\end{eqnarray} \tag{ 4 }$$

which we use with the following assumptions to approximate our observational detection limit in this case: zero eccentricity, 90° or 45° inclination, 0.95 ${{ \mathcal M }}_{\odot }$ stellar mass, and a K* limit of 2 m s⁻¹, based on our stellar jitter and instrument rms values. In Figure 14, we show this estimated detection limit in planet mass–period space; the solid curve is assuming a 90° inclination and the dashed curve a 45° inclination for the planetary orbit. As expected, a 0.018 ${{ \mathcal M }}_{{\rm{J}}}\sim 5.7\,{{ \mathcal M }}_{\oplus }$ planet at 5.88 days ∼0.016 years is at the limit of detectability with our current data, even assuming the maximum potential RV signal. However, we could expect to detect a ∼20 ${{ \mathcal M }}_{\oplus }$ planet in an orbit of a few hundred days, for instance, given our current data and RV precision and the assumptions made above.

Zoom In Zoom Out Reset image size

5.1. Stellar Parameters

Here we derive stellar parameters—effective temperature (T_eff), surface gravity (log g), metallicity ([Fe/H]), and microturbulent velocity (ξ)—with iodine-free spectra from PFS; the instrument and data reduction are described below. A traditional ionization–equilibrium balance method was implemented, in which the correlations between [Fe i/H] and the excitation potential (χ) and between [Fe i/H] and the reduced equivalent width [log($\mathrm{EW}/\lambda$)] are minimized. The difference between the mean iron abundances measured from Fe i and Fe ii lines is also minimized. This procedure is done iteratively to derive the "best" spectroscopic stellar paramters. We used the publicly available Qoyllur-quipu (q²) Python package⁹ , described in detail in Ramírez et al. (2014), a wrapper for the spectral analysis code MOOG (Sneden 1973; Sobeck et al. 2011), which takes equivalent widths (EWs) as input to derive stellar parameters and abundances.

EW measurements were made by manually fitting Gaussian functions to observed line profiles using the splot task in IRAF¹⁰ , with the iron line list from Ramírez et al. (2014) that includes Fe i lines from a wide range of excitation potentials and 16 Fe ii lines. The EWs are measured differentially with respect to the Sun on a line-by-line basis; our solar spectrum comes from PFS observations on 2016 January 7 UT of reflected sunlight from the asteroid Vesta, and we assume ${T}_{\mathrm{eff},\odot }=5777$ K, log g_⊙ = 4.44, [Fe/H]${}_{\odot }$ = 0, and ${\xi }_{\odot }$ = 1 km s⁻¹. This differential approach reduces the impact of uncertainties in stellar models and in atomic data, as they are canceled out in each star−Sun measurement. The continuum regions for each line were chosen to be "clean" and the same for the two stars and the Sun; this is possible because HD 133131A and B are both inactive G2 stars and the data are taken with the same spectrograph set-up. EWs were translated into abundances using the abfind driver in MOOG with 1D-LTE model atmospheres linearly interpolated from the marcs grid.

The derived [Fe i/H] and [Fe ii/H] values for HD 133131A and B have line-to-line scatters of 0.037/0.032 dex (71 Fe i lines) and 0.041/0.041 dex (16 Fe ii lines), respectively. q² automatically calculates errors on T_eff, log g, and ξ as in Epstein et al. (2010) and Bensby et al. (2014), and [Fe/H] errors by adding the other stellar parameter errors in quadrature (under the assumption that they are independent) with the standard error of the mean line-to-line scatter ($\sigma /\sqrt{N-1}$) of the [Fe/H] values. The errors reported here are a reflection only of how well the minimization criteria above are met, and are still subject to systematic errors, e.g., 3D and NLTE effects (Asplund 2005) not captured in our 1D-LTE analysis. Comparisons of abundances measured in 1D versus 3D stellar atmosphere models can differ at the 0.01–0.02 dex level, and have been shown to increase with lower [Fe/H] (e.g., Ramírez et al. 2008; Bergemann et al. 2012; Magic et al. 2014). However, these differences are measured in an absolute sense, so that they are subject to systematic uncertainties in a way that our strictly differential analysis is not, especially because HD 133131A and B are very similar to the Sun except for their metallicities (Table 6).

Table 6.  Stellar Parameters

Star	${T}_{\mathrm{eff}}$	log g	$[\mathrm{Fe}/{\rm{H}}]$	ξ
	(K)	(cgs)	(dex)	(km s−1)
Desidera et al. (2006b)
HD 133131A	5745	4.46	−0.33	1.05
HD 133131B	5739	4.46	−0.36	1.05
Δ(A-B)	+6	0	+0.03	0
Δ(A-B)a	+5 ± 15	⋯	+0.032 ± 0.015	⋯
This work, Solar Reference
HD 133131A	5799 ± 19	4.39 ± 0.050	−0.306 ± 0.016	1.10 ± 0.040
HD 133131B	5805 ± 15	4.41 ± 0.045	−0.281 ± 0.013	1.12 ± 0.030
Δ(A-B)	−6 ± 24	−0.02 ± 0.064	−0.025 ± 0.021	−0.02 ± 0.050
This Work, HD 133131A Reference
HD 133131A (same as solar ref.)	5799 ± 19	4.39 ± 0.050	−0.306 ± 0.016	1.10 ± 0.040
HD 133131B	5811 ± 12	4.42 ± 0.032	−0.275 ± 0.010	1.11 ± 0.020
Δ(A-B)	−12 ± 22	−0.03 ± 0.059	−0.031 ± 0.019	−0.01 ± 0.045
This Work, HD 133131B Reference
HD 133131A	5796 ± 13	4.40 ± 0.032	−0.311 ± 0.011	1.12 ± 0.030
HD 133131B (same as solar ref.)	5805 ± 15	4.41 ± 0.045	−0.281 ± 0.013	1.12 ± 0.030
Δ(A-B)	−9 ± 20	−0.01 ± 0.055	−0.030 ± 0.017	0.00 ± 0.042

Note.

^aAs reported in their Table 11, final analysis with propogated errors.

Download table as:  ASCIITypeset image

The stellar parameters derived here from a differential analysis with the Sun are listed in the second section of Table 6, and agree moderately well with the parameters found by Desidera et al. 2006a), who use a similar ionization–equilibrium balance to derive T_eff and [Fe/H] but derive their log g values using photometric information. Interestingly, we find the stars to be slighly closer to solar in T_eff and [Fe/H], and also that the A component is cooler and more metal-rich than the B component, the opposite of Desidera. However, given our errors, we find the A and B stellar parameters are insignificantly different in everything except [Fe/H] (with Δ(A-B) [Fe/H] = −0.025 ± 0.021).

The metal-poor nature of this system puts both stars outside the "solar twin" realm, and although it is nearby (∼50 pc), the system may have experienced slight differences in Galactic chemical evolution (especially if is much older; see above). As in previous stellar "twin" papers (e.g., Ramírez et al. 2015; Saffe et al. 2015; Teske et al. 2016), here we also derive more precise stellar parameters by comparing the two stars strictly against each other. In previous works where such a strict differential method is used to compute stellar parameters, only the "cool"–"hot" (here, A-B) case is derived. In this paper we compare both (A-B) and (B-A) stellar parameters, since both stars are so similar to the Sun. In each case (A-B, B-A), we assume the (B, A) solar reference parameters and measure only the (A, B) T_eff, log g, Δ[Fe/H], and ξ. These strictly A versus B differential parameters are listed in the third and fourth sections of Table 6. In both cases, the alternative parameters are very similar to the solar reference parameters, but the errors are reduced, particularly in the (A-B) case (last section of the table).

The median S-index value for HD 133131A is 0.155, while the median for HD 133131B is 0.154 (see Tables 1–3). Taking these median S-index values and converting them to ${R}_{\mathrm{HK}}^{\prime }$ using the Noyes et al. (1984) calibration results in ${R}_{\mathrm{HK}}^{\prime }=-4.913$ for A and −4.919 for B. Several less recent chromospheric age relations (Soderblom et al. 1991, Equations (1) and (3); Donahue 1993) indicate that stars with ${R}_{\mathrm{HK}}^{\prime }\sim -4.9$ should be $\gtrsim 9.5\,\mathrm{Gyr}$ old, as does Figure 1 of Pace (2013, for the T_eff of HD 133131A & B derived here). However, the age–metallicity relation of Rocha-Pinto et al. (2000), for [Fe/H] $\sim \,-0.3$, suggests a slightly younger age of ∼7 Gyr. The recent work of Bergemann et al. (2014) measuring the age–metallicity relation of Milky Way disk stars as part of the Gaia-ESO survey also suggests an age closer to ∼9.5 Gyr for a star with the [Fe/H], [Mg/Fe] (∼0.03) and T_eff of HD 133131A & B (see their Figures 6 and 7). The kinematics of the pair suggest with ∼93% probability thin disk membership, based on UVW velocities (calculated from the coordinates, proper motions, and parallax of van Leeuwen 2007, and the absolute RV of Gontcharov 2006) and the relations of Reddy et al. (2006, Equations (1) and (2)). We conclude that the pair are coeval, and likely older than the Sun, perhaps over twice as old.

5.2. Stellar Abundances

In the spectra of Vesta, HD 133131A, and HD 133131B we measured absorption lines of 19 elements in addition to Fe, and combined these EW measurements with our derived stellar parameters in a curve-of-growth analysis within MOOG to derive stellar abundances. The procedure was the same as that used in Teske et al. (2016), involving an examination of the normalized, Doppler-corrected spectra of every line with the splot task in IRAF and choosing regions for the continuum that were clean and the same in all three spectra. Carbon abundances were derived from both C i and CH features, and lines from both neutral and singly ionized species were measured for Sc, Ti, and Cr. We applied hyperfine structure corrections to the V, Mn, Co, Cu, Rb, Y, and Ba abundances. The measured EW values for each element, including Fe, are listed in Table 7. The abundances and total errors for each abundance, including line-to-line scatter and the errors propogated from each parameter uncertainty, are listed in Table 8. In the case of CH, Al i, and Zr ii, where only one line is measured, we adopted the largest line-to-line scatter value among the other elements with >3 lines, and added this value in quadrature with the errors propagated from the stellar parameter errors.

Table 7.  Measured Lines and Equivalent Widths

Ion	λ	χ	log gf	EW${}_{\odot }$	EW HD 133131A	EW HD 133131B
	(Å)	(eV)	(dex)	(mÅ)	(mÅ)	(mÅ)
Fe i	4389.245	0.052	−4.583	71.5	65.2	61.6
Fe i	4602.001	1.608	−3.154	72.4	59.6	63.1
Fe i	4690.140	3.69	−1.61	57.7	44.1	45.7
Fe i	4788.760	3.24	−1.73	70.1	58.2	57.8
Fe i	4799.410	3.64	−2.13	35.6	22.5	24.5
Fe i	4808.150	3.25	−2.69	28.1	17.5	18.0
Fe i	4950.100	3.42	−1.56	76.8	59.2	60.5
Fe i	4994.129	0.915	−3.08	103.6	89.7	91.7
Fe i	5141.740	2.42	−2.23	90.0	72.1	73.4
Fe i	5198.710	2.22	−2.14	97.3	85.0	85.2
Fe i	5225.525	0.11	−4.789	73.2	60.0	61.3
Fe i	5242.490	3.63	−0.99	88.0	72.7	73.8
Fe i	5247.050	0.087	−4.961	65.9	53.2	54.6
Fe i	5250.208	0.121	−4.938	68.0	52.8	54.5
Fe i	5295.310	4.42	−1.59	29.6	17.2	18.5
Fe i	5322.040	2.28	−2.89	61.5	46.5	47.9
Fe i	5373.710	4.47	−0.74	62.4	48.1	49.8
Fe i	5379.570	3.69	−1.51	60.9	46.6	48.3
Fe i	5386.330	4.15	−1.67	32.3	20.0	21.5
Fe i	5441.340	4.31	−1.63	31.2	19.5	20.1

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Table 8.  Derived Stellar Abundances

Species	Tc	A-Vesta Params		B-Vesta Params		B-A Params		A-B Params
		Δ[X/H]	error	Δ[X/H]	error	Δ[X/H]	error	Δ[X/H]	error
	(K)	(dex)	(dex)	(dex)	(dex)	(dex)	(dex)	(dex)	(dex)
C i	40	−0.250	0.042	−0.253	0.021	−0.005	0.025	0.007	0.025
CHa	40	−0.298	0.091	−0.293	0.086	0.014	0.057	−0.013	0.058
O i b	180	−0.242	0.037	−0.247	0.032	−0.002	0.021	0.007	0.022
Na i	958	−0.290	0.025	−0.292	0.014	0.002	0.016	−0.001	0.017
Mg i	1336	−0.279	0.024	−0.250	0.018	0.033	0.013	−0.033	0.013
Al i a	1653	−0.283	0.089	−0.280	0.084	0.006	0.056	−0.005	0.056
Si i	1310	−0.281	0.009	−0.260	0.008	0.024	0.006	−0.022	0.007
S i	664	−0.266	0.018	−0.251	0.021	0.014	0.025	−0.011	0.025
Ca i	1517	−0.256	0.017	−0.241	0.015	0.021	0.012	−0.022	0.013
Sc i	1659	−0.255	0.044	−0.255	0.024	0.007	0.028	−0.003	0.028
Sc ii	1659	−0.271	0.025	−0.267	0.019	0.009	0.015	−0.004	0.016
Ti i	1582	−0.272	0.021	−0.236	0.017	0.043	0.012	−0.040	0.013
Ti ii	1582	−0.266	0.025	−0.245	0.025	0.026	0.017	−0.021	0.018
V i	1429	−0.341	0.029	−0.299	0.026	0.049	0.016	−0.046	0.017
Cr i	1296	−0.313	0.017	−0.288	0.014	0.030	0.011	−0.028	0.012
Cr ii	1296	−0.317	0.025	−0.305	0.022	0.015	0.012	−0.011	0.012
Mn i	1158	−0.476	0.034	−0.449	0.028	0.034	0.012	−0.033	0.013
Fe i	1334	−0.305	0.037	−0.281	0.032	0.030	0.024	−0.030	0.024
Fe ii	1334	−0.309	0.041	−0.280	0.041	0.031	0.038	−0.028	0.038
Co i	1352	−0.302	0.025	−0.276	0.02	0.032	0.012	−0.028	0.013
Ni i	1353	−0.316	0.014	−0.294	0.012	0.028	0.009	−0.026	0.009
Cu i	1037	−0.344	0.064	−0.331	0.051	0.020	0.018	−0.020	0.020
Zn i	726	−0.340	0.024	−0.349	0.023	−0.005	0.007	0.005	0.008
Y ii	1659	−0.414	-0.034	−0.375	0.030	0.045	0.015	−0.040	0.016
Zr ii a	1741	−0.360	0.091	−0.308	0.087	0.057	0.057	−0.050	0.058
Ba ii	1455	−0.334	0.059	−0.323	0.060	0.020	0.011	−0.021	0.015

Notes.

^aThese elements have only one line measured, so the line-to-line to scatter value adopted in their error calculation is the greatest line-to-line scatter from the rest of the elements with >3 lines measured. ^bThe errors on the oxygen abundances were conservatively estimated as the average of all the other elemental abundance errors, within each parameter set.

Download table as:  ASCIITypeset image

The only elemental abundance not derived directly from EW measurements was [O/H]. The oxygen triplet at 7775 Å is not included in the wavelength range of PFS spectra, so we performed a synthesis analysis on the [O i] line at 6300 Å using the MOOG synth driver, as described in Teske et al. (2014). We adopt a conservative error of the average of the other elemental abundance errors, for each set of parameters (A-Vesta, B-Vesta, B-A, A-B; see Table 8).

In Figures 15 and 16, we show the abundance differences derived with respect to Vesta, and derived for each star with respect to the other, plotted against the 50% condensation temperatures (T_c) from Lodders (2003) for solar compostion gas. In each case, the corresponding stellar parameters were used to derive the abundances—e.g., in the A-B abundances, the A-B stellar parameters from Table 6 were used. The dotted lines in Figure 15 show the weighted means of the A-Vesta (orange, lower value) and B-Vesta (blue, higher value) abundances; HD 133131A is more metal-poor overall. There are slight differences between the weighted means of elements above and below T_c = 1000K, with a ∼0.01 dex decrease in more refractory (higher T_c) element abundances, but no obvious pattern like that seen in Meléndez et al.'s (2009) analysis of the Sun versus other solar twins. In that work, Meléndez et al. suggested that the deficit of refractory elements in the Sun with respect to other solar twins was a signature of small planet formation. More recent work (e.g., Adibekyan et al. 2014; Nissen 2015; Schuler et al. 2015; Spina et al. 2016) has shown that trends with T_c may not be (entirely) due to planet formation, but instead to the birthplace and age of a star.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can instead examine the Δ[X/H] values, derived from a strict comparison of one star in the HD 133131 system versus the other, to avoid potential causes of abundance differences like age and birthplace; we omit the uncertainty introduced by the Sun and comparing HD 133131 A and B to a dissimilar star (at least in [Fe/H]). In Figure 16, these Δ[X/H] values are plotted against T_c, with a green dashed line showing zero difference and dotted lines indicating the weighted means of elements with ${T}_{c}\lt 1000$ K and ${T}_{c}\gt 1000$ K, the approximate demarcation between volatile and refractory elements. We show both A-B (orange circles) and B-A (blue stars) cases to demonstrate the good agreement between both stellar parameter derivations; the EWs are the same in both abundance determinations. In previous stellar "twin" abundance studies (Ramírez et al. 2011; Liu et al. 2014; Mack et al. 2014; Tucci Maia et al. 2014; Ramírez et al. 2015; Saffe et al. 2015; Teske et al. 2015, 2016), the Δ[X/H] abundances are always quoted from the "cold–hot" star comparison. There is now an obvious decrease (A-B)/increase (B-A) in the Δ[X/H] values moving from volatile to refractory elements. The weighted means of the ${T}_{c}\lt 1000$ K abundances are 0.0032/–0.0027 dex, whereas the weighted means of the ${T}_{c}\gt 1000$ K abundances are −0.0255/0.0274 dex.

As demonstrated by the RV observations and planet orbital fitting in Sections 3 and 4, HD 133131A hosts two planets of moderate eccentricity, a 1.43 ± 0.04 ${{ \mathcal M }}_{{\rm{J}}}$ minimum mass planet at 1.44 ± 0.005 au and a 0.63 ± 0.15 ${{ \mathcal M }}_{{\rm{J}}}$ minimum mass planet at 4.79 ± 0.92 au (assuming the low-eccentricity configuration for planet c). HD 133131B hosts one planet of relatively high eccentricity (0.62 ± 0.04), with a mimimum mass of 2.50 ± 0.05 ${{ \mathcal M }}_{{\rm{J}}}$ orbiting at 6.40 ± 0.59 au. The two stars are separated by only ∼360 au (Desidera et al. 2006a), making them the most closely separated "twin" pair (${\rm{\Delta }}{T}_{\mathrm{eff}}\lesssim 100\,{\rm{K}}$) with detected planets (the next closest are 16 Cyg A/B and HD 80606/7 at ∼1000 au). This system is even more rare in that both "twins" host planets. Interestingly, although previous studies have found that close-in giant planet host stars are more likely to have stellar companions with separations between 50 and 2000 au versus field stars (Ngo et al. 2015), broader planet occurrence rates seem to decrease for stellar companions as a function of separation; Wang et al. (2014b) found that stars with a companion at 100 au are 2.6 ± 1 times less likely to host a planet, and stars with a companion at 1000 au are 1.7 ± 0.5 times less likely to host a planet. Furthermore, stellar multiplicity seems to be negatively correlated with the presence of multiple transiting planets, although it is unclear for separations similar to that of HD 133131A and B whether this is due to the suppression of planet formation (Roell et al. 2012; Wang et al. 2014a, 2014b; Touma & Sridhar 2015; Kraus et al. 2016) or the inclination of planets being perturbed by stellar companions (Wang et al. 2014b, see Figure 7; Wang et al. 2015).

It is also noteworthy that both stars are relatively metal-poor (with [Fe/H] ∼ −0.30 dex); their metallicities fall below ∼90% of measured planet host metallicites, and they join the sample of only six other such metal-poor exoplanet host stars in binary systems (based on exoplanets.org). Thus, based only on what has been measured, the probability of HD 133131A/B hosting planets is <1%. The most recent functional forms of the giant planet metallicity correlation predict stars with [Fe/H] ∼ −0.30 dex are ∼3.5 times less likely to host planets than solar metallicity stars, and ∼8.5 times less likely to host planets than stars with [Fe/H] = 0.30 dex (Mortier et al. 2013). Thus, these giant planets around metal-poor stars seem to be rare, at least based on current statistics of 0.1–25 ${{ \mathcal M }}_{{\rm{J}}}$ planets at 5–5000 day periods. The metal-poor nature of HD 133131A and B also places them in a sparsely populated part of planet e versus host star [Fe/H] space—only 5 (11) other exoplanets have both $e\gt 0.3$ ($e\gt 0.2$) and host star [Fe/H] ≤ −0.30 (out of 80 [Fe/H] ≤ −0.30 host stars, based on exoplanets.org).

As discussed in Section 2, HD 133131A and B have different abundances relative to each other, for elements with ${T}_{c}\gt 1000$ K, by about 0.025 dex, with A more metal-poor than B (see Figure 16). Past studies of "twin" planet hosting stars have used such measured abundance differences to try to deduce information about planet composition and/or formation (Ramírez et al. 2011, 2015; Tucci Maia et al. 2014; Teske et al. 2015, 2016). We can explore similar arguments here. However, these explorations depend significantly on the assumed convection zone masses of stars, which change over the lifetime of the star (e.g., Bahcall et al. 2001).

For instance, we can estimate whether the difference in mass between the two planets around HD 133131A versus around B can be explained by the observed abundance differences using the corrected formula from Ramírez et al. (2011):

$$\begin{eqnarray}&&{\rm{\Delta }}[M/H]=\mathrm{log}\left[\displaystyle \frac{{(Z/X)}_{{cz}}{{ \mathcal M }}_{{cz}}+{(Z/X)}_{P}{{ \mathcal M }}_{p}}{{(Z/X)}_{{cz}}({{ \mathcal M }}_{{cz}}+{{ \mathcal M }}_{p})}\right].\end{eqnarray} \tag{ 5 }$$

The current convection zone mass of star A or B ${{ \mathcal M }}_{{cz}}$ is estimated to be 0.026 ${{ \mathcal M }}_{\odot }$ from the relation in Pinsonneault et al. (2001).¹¹ The planet metallicity, ${(Z/X)}_{p}$, is estimated to be 0.1, although this value can vary from ∼0.03–0.5 for planets with masses ∼0.4–2.5 ${{ \mathcal M }}_{{\rm{J}}}$ (see Thorngren et al. 2015, Figure 8.). The mass difference between Ab+Ac and Bb, ${{ \mathcal M }}_{p}$ is 0.44 ${{ \mathcal M }}_{{\rm{J}}}$ (assuming a low e c planet, and ignoring the effects of inclination in ${ \mathcal M }\,\sin \,i$, which on average should be the same for both stars), and we estimate the stellar convection zone metallicity ${(Z/X)}_{{cz}}$ by scaling Asplund et al.'s (2009) ${(Z/X)}_{\odot }=0.134$ to the system [Fe/H] (−0.30 dex). Then the required depletion/enhancement expected in HD 133131A/B is 0.09 dex, three times larger than the observed abundance difference (although this value can range from 0.04 to 0.13 if we include ${{ \mathcal M }}_{{ \mathcal P }}$ errors). If instead we assume ${{ \mathcal M }}_{{cz}}$ is 0.15 ${{ \mathcal M }}_{\odot }$, estimating a solar-type star's convection zone size at 15–20 Myr according to the standard model of Serenelli et al. (2011), then the required depletion/enhancement in A/B to explain the planet mass difference is only 0.02 dex, similar to the observed ∼0.025 dex difference. We can instead perform the reverse calculation, and ask what mass difference explains a 0.025 dex abundance difference, which results in 0.11 ${{ \mathcal M }}_{{\rm{J}}}$ for ${{ \mathcal M }}_{{cz}}=0.026{{ \mathcal M }}_{\odot }$ and 0.65 ${{ \mathcal M }}_{{\rm{J}}}$ for ${{ \mathcal M }}_{{cz}}=0.15{{ \mathcal M }}_{\odot }$; the latter value is closer to the observed mass difference between the planets around HD 133131A and HD 133131B. Referring back to Figure 14, a 0.11 ${{ \mathcal M }}_{{\rm{J}}}\sim 35\,{{ \mathcal M }}_{\oplus }$ planet with a period of up to ∼3 years should be detectable in our data, and a 0.65 ${{ \mathcal M }}_{{\rm{J}}}\sim 200\,{{ \mathcal M }}_{\oplus }$ planet with a period of up to ∼tens of years should be detectable in our data, including the assumptions of zero e, edge-on orientation, and a ${K}_{* }\sim 2$ m s⁻¹ limit to our RV precision. The non-detection of these planets, given our detection limit approximation and current data, suggests that planets of these masses are on longer orbits, if they exist.

In the estimates above, the motivation is to account for the abundance differences between the two stars with refractory-rich mass "missing" from HD 133131A or "added to" HD 133131B at some point during the planet formation process. The results change by an order of magnitude depending on the assumed convection zone sizes of the stars, which is related to the time at which refractory material is accreted/depleted. Perhaps the good match with observations with an assumed larger ${{ \mathcal M }}_{{cz}}$, and thus earlier accretion/depletion, is a sign that this is the more likely scenario. We note that none of the simple "toy model" calculations above account for the effect of gravitational settling of heavy elements (e.g., Bahcall et al. 1995), enhanced mixing (e.g., Pinsonneault et al. 1989), or radiative levitation (e.g., Pinsonneault 1997) that are known to occur in the Sun.

Instead of assuming that the HD133131 planets have the same heavy element mass (${(Z/X)}_{p}=0.1$), we can improve this estimate by considering the relation between planet mass, and total heavy element mass derived by Thorngren et al. (2015):

$$\begin{eqnarray}&&{{ \mathcal M }}_{z}=(46\pm 5.2){{ \mathcal M }}_{{ \mathcal P }}{}^{0.59\pm 0.073}\end{eqnarray} \tag{ 6 }$$

where ${{ \mathcal M }}_{z}$ is the heavy element mass in ${{ \mathcal M }}_{\oplus }$ in a planet of mass ${{ \mathcal M }}_{{ \mathcal P }}$ in ${{ \mathcal M }}_{{\rm{J}}}$. Including the errors in the constants and the uncertainties in our derived ${ \mathcal M }\,\sin \,i$ values, planet Ab, Ac, and Bb could have ${{ \mathcal M }}_{z}$ ranging from 48–69, 30–43, and 65–95 ${{ \mathcal M }}_{\oplus }$, respectively. Taking the middle of each of these ranges, the heavy element mass in Ab+Ac ∼ 95 ${{ \mathcal M }}_{\oplus }$, and in Bb is ∼80 ${{ \mathcal M }}_{\oplus }$, or ${\rm{\Delta }}{ \mathcal M }=15\,{{ \mathcal M }}_{\oplus }\sim 0.05$ ${{ \mathcal M }}_{{\rm{J}}}$ of heavy element material that should be in the HD 133131A planets (versus the star) and not in the HD 133131B planet (but in the star). If we use Equation (6) to solve for ${{ \mathcal M }}_{{ \mathcal P }}$, assuming ${{ \mathcal M }}_{z}=15$ ${{ \mathcal M }}_{\oplus }$, this $=0.15{{ \mathcal M }}_{{\rm{J}}}={{ \mathcal M }}_{p}$ in the denominator of Equation (5). Also assuming 15 ${{ \mathcal M }}_{\oplus }$ for ${(Z/X)}_{P}{{ \mathcal M }}_{p}$ in Equation (5) , the resulting Δ[M/H] is then 0.10 dex for ${{ \mathcal M }}_{{cz}}=0.026$ ${{ \mathcal M }}_{\odot }$, and 0.02 dex for ${{ \mathcal M }}_{{cz}}=0.15$ ${{ \mathcal M }}_{\odot }$. The observed Δ[M/H] (0.025) falls in between these values, again perhaps giving a glimpse of how the timing of planet formation and stellar convection zone size are related. Overall, the Δ[M/H] between HD 133131A and B may be an artifact of the different interior compositions of their planets. Again referring back to Figure 14, our detection limit estimates suggest that a 0.04 ${{ \mathcal M }}_{{\rm{J}}}\sim 13\,{{ \mathcal M }}_{\oplus }$ planet or a 0.11 ${{ \mathcal M }}_{{\rm{J}}}\sim 35\,{{ \mathcal M }}_{\oplus }$ planet must be in orbits with periods longer than ∼a few dozen days or ∼3 years, respectively, since we do not detect them.

Alternatively, if values of ${{ \mathcal M }}_{z}$ for Ab, Ac, and Bb are on the low end of our estimates (48–69, 30–43 and 65–95 ${{ \mathcal M }}_{\oplus }$), the two systems could have the same amount of heavy element mass sequestered in their respective planets, leaving the Δ[M/H] unaccounted for. There may have been small planets or rocky material that was once in either system, but due to planet–planet scattering was ejected (in the case of A) or accreted on to the star (in the case of B). Ford & Rasio (2008) discuss in detail planet–planet scattering in systems with two unequal-mass planets, and find that collisions are less frequent versus scattering between equal-mass planets. Ejection also seems more likely, versus planet collision, given the relatively large periods of the detected planets (Petrovich 2015). The "jumping Jupiter" scenario, to account for the observed angular momentum deficit in the terrestrial planets despite Jupiter crossing its 2:1 mean motion resonance with Saturn (Brasser et al. 2009; Minton & Malhotra 2009; Morbidelli et al. 2009, 2010; Walsh & Morbidelli 2011; Agnor & Lin 2012), suggests that (1) a planet like Uranus or Neptune scattered off Jupiter (Thommes et al. 1999; Nesvorný 2011; Batygin et al. 2012), and (2) that the angular momentum diffusion of the terrestrial planets could have been significant (Agnor & Lin 2012; Brasser et al. 2013). In fact, a large suite of N-body simulations show a $\geqslant$ 85% probability that at least one terrestrial planet is lost (accreted on to the star or ejected) in the process of the solar system giant planet instability (Kaib & Chambers 2016). Ford & Rasio (2008) also point out that an important parameter for determining whether giant planets have circular or eccentric orbits may be the timing of the final planet–planet scattering event, whether there is still enough material in the disk to damp eccentricities (see also Ida & Lin 2008; Raymond et al. 2009). It is important to keep in mind the relatively close separation of HD 133131A and B (∼360 au), and that this could have affected the disk sizes, orientations, stabilities/lifetimes (see Zhou et al. 2012 for a brief overview).

There is always the possibility that the Δ[X/H] differences observed between the two stars are not related to planet formation. While we avoid some complicating factors like different degrees of dust cleansing by luminous stars (Önehag et al. 2011, 2014) and different ages (Adibekyan et al. 2014; Nissen 2015; Spina et al. 2016) by analyzing "twin" stars in a binary, there is a chance that gas–dust segregation (Gaidos 2015) may have happened differently in the two stellar disks, and/or that the stellar composition difference is an artifact of formation of the stars. Larger surveys of binary stars suggest that [Fe/H] differences $\gtrsim 0.03$ are atypical (Desidera et al. 2004, 2006a), and here Δ[Fe/H]${}_{{\rm{A}} \mbox{-} {\rm{B}}}$ is −0.030 ± 0.017. Gratton et al.'s (2001) study of a wider suite of abundances in six wide binaries found four with abundances equal at the ∼0.012 dex level, but suggested the other two stars could be planet hosts based on their larger abundance differences. Furthermore, the two stars in the Gratton sample with the lowest ${\rm{\Delta }}{T}_{\mathrm{eff}}$ (at 11 ± 32 and 90 ± 13 K), have Δ[Fe/H] values of 0.000 ± 0.019 and −0.005 ± 0.009, respectively. Our even smaller ${\rm{\Delta }}{T}_{\mathrm{eff}}=-9\pm 20$ K suggests the Δ[Fe/H]${}_{{\rm{A}} \mbox{-} {\rm{B}}}$ difference we measure is meaningful. We also note that while Nordström et al. (2004) list the $v\,\sin \,i$ values for HD 133131A and B as 4 and 7 km s⁻¹, respectively, we see no difference in the width of absorption lines in the two components across many unblended lines, indicating differences in rotation are not resposible for differences in abundance measurements.

Speculation about the origin of the Δ[M/H] abundance differences aside, the context of the Thorngren et al. (2015) study again emphasizes how unusual the HD 133131A and B planets are, based on their observed masses and host star metallicities. Though not well defined, there is a positive relationship between planet heavy element mass and host star [Fe/H]; a bootstrap regression gives ${{ \mathcal M }}_{z}=(31.4\pm 3.4)\,\times {10}^{0.48\pm 0.047[\mathrm{Fe}/{\rm{H}}]}$. This relation predicts ∼23 ${{ \mathcal M }}_{\oplus }$ of heavy element mass in the HD 133131A and B planets, much lower than most of the estimates from Equation (6) above. Understanding how giant planet heavy element mass depends on other elements may help shed light on the disagreement between these two predictions. Indeed, HD 133131A and B both have higher abundances of C, O, Mg, Si, and Ti than Fe; further studies of variation in these host stars abundances with giant planet heavy enrichment are needed.

In this work we present the detection of three giant planets, two orbiting the A component of HD 133131 and one orbiting the B component. This is the first planet detection paper based primarily on Magellan/PFS RV observations, and demonstrates the ∼1 m s⁻¹ precision of this instrument over its six-year baseline. The outer planets around HD 133131A and B also push the boundary of RV planet detection—only around a dozen other RV planets have been detected at periods $\geqslant 3600$ days. A full analysis of the frequency of eccentric giant planets at long periods is beyond the scope of this paper, but we note that of the 52 planets detected at ${ \mathcal P }\gt 5$ years and with ${ \mathcal M }\gt 0.3{{ \mathcal M }}_{{\rm{J}}}$ (∼Saturn mass), 18 have $e\gt 0.3$ and 30 have $e\gt 0.2$ (exoplanets.org); the detections presented here add two (three in the $e\gt 0.2$ case) planets to that group, implying an uncorrected frequency of ∼40% (or 60% for $e\gt 0.2$).

Our careful differential analyses, both comparing the stars to the Sun (A-Vesta, B-Vesta) and comparing the two stars to each other (B-A, A-B), provide precise stellar parameters and abundances that reveal that HD 133131A and B are solar "twins" in T_eff, log g, and ξ, but are more metal-poor ([Fe/H] for both ∼−0.30) and likely older than the Sun (∼9.7 Gyr). Additionally, we find a small but significant depletion of high-T_c ($\gt 1000$ K) elements in HD 133131A versus B, and explore how this could be related to differences in planet formation, interactions between the planets and the two stars, and/or the composition of the planets orbiting the two stars. This system is the smallest separation "twin" binary system for which such a detailed abundance analysis has been conducted. Overall, HD 133131A and B are especially noteworthy because they are metal-poor but are orbited by multiple giant planets, contradictory to the predicted giant planet–metallicity correlation. The planets detected here will be important benchmarks in studies of host star metallicity, binarity of host stars, and long-period giant planet formation.

RPB gratefully acknowledges support from NASA OSS Grant NNX07AR40G, the NASA Keck PI program, and from the Carnegie Institution of Washington. MD gratefully acknowledges the support of CONICYT-PFCHA/Doctorado Nacional, Chile. The work herein is based on observations obtained at Las Campanas Observatory of the Carnegie Institution of Science. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France, and the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. We thank the referee for a detailed reading of our paper and helpful comments, which improved the content and quality of the paper.

Facility: Magellan:Clay (MIKE, PFS).

Software: Qoyllur-quipu (Ramírez et al. 2014), IRAF, MOOG (Sneden 1973 Sobeck et al. 2011), SYSTEMIC (Meschiari et al. 2009).