A PLANETARY SYSTEM AROUND THE NEARBY M DWARF GJ 667C WITH AT LEAST ONE SUPER-EARTH IN ITS HABITABLE ZONE

The Doppler detection of extrasolar planets is achieved by measuring the periodic radial velocity (RV) variations induced in a star by the presence of orbiting low-mass companions. The Doppler signature of a gas giant planet with orbital parameters similar to Jupiter is about 10 m s⁻¹ over a period of 11 years. By comparison, Earth's reflex barycentric pull on the Sun corresponds to 8 cm s⁻¹. The most precise spectrographs can deliver long-term stability at the level of 1–3 m s⁻¹(Mayor et al. 2011; Vogt et al. 2010). This precision is sufficient to detect candidates of a few Earth masses in tight orbits around low-mass stars (M dwarfs). A key requirement in reaching such accuracy is the extraordinary calibration of both the wavelength scales of the spectrum and the instrumental point-spread function. The two most successful methods used to date are (1) the iodine cell technique (Butler et al. 1996) and (2) the construction of stabilized spectrographs fed by optical fibers (Baranne et al. 1996).

In the first of these two approaches, the star's light is passed through a transparent cell containing iodine gas at low pressure. The absorption spectrum of iodine is imprinted on the star's light, tracing precisely the same optical path as encountered by the star light in traversing the spectrometer. The combined spectrum of the star and iodine is then modeled to obtain a simultaneous point-spread function, wavelength, and Doppler shift solution for the stellar spectrum (Butler et al. 1996). The stabilized spectrograph approach relies on the construction of a vacuum-sealed spectrometer fed with optical fibers that produces a fairly constant instrumental profile and long-term wavelength stability. Each night, stabilized spectrographs are calibrated in wavelength by feeding the light from a wavelength standard source (e.g., Th/Ar lamp) through the same fiber as the science targets. This approach is exquisitely implemented by HARPS installed on the European Southern Observatory (ESO) 3.6 m telescope at La Silla Observatory (Pepe et al. 2003). HARPS is a cross-dispersed echelle spectrograph with a resolving power of λ/δλ = 110000, and over the years has demonstrated 1 m s⁻¹ long-term stability (Pepe et al. 2011; Mayor et al. 2011).

Even though HARPS is probably the most precise astronomical spectrometer ever built, the Cross-Correlation Function (CCF) data analysis method that has been commonly used to analyze this data is suboptimal in the sense that it does not exploit the full Doppler information in the stellar spectrum (Pepe et al. 2002). For this reason, instead of using the CCF RVs provided by the ESO archive, we use a least-square template-matching method to derive new RV measurements. Thanks to the instrumental stability and the excellent wavelength calibration provided by the HARPS-ESO data reduction software, the model required to match each observation to a high signal-to-noise-ratio (S/N) template only needs to include a Doppler offset and a multiplicative polynomial to correct for the flux variability across each echelle order. The template is obtained by co-adding all the spectra after a preliminary RV measurement is obtained using the highest S/N observation. The least-squares matching technique has been used on HARPS data before. An example is the RV measurements on GJ 1214 (V = 14.57) used to derive the mass of the transiting super-Earth reported in Charbonneau et al. (2009). The performance and description of our software tool, HARPS-TERRA (Template Enhanced Radial velocity Re-analysis Application) on a representative sample of stars can be found in Anglada-Escudé & Butler (2012).

One of the M dwarfs with many public HARPS spectra but no published detections is GJ 667C. The ESO archive contains 143 observations of this star obtained between 2004 June and 2008 October. Typical exposure time is between 900 and 1500 s and the average S/N is 64 at 6100 Å. At a 2009 conference, a planet candidate orbiting this star was announced with a period P of ∼7 days. Also recently, Bonfils et al. (2011) reported the detection of a plausible signal with P ∼ 28 days similar to one of the candidates we report here. However, neither detailed analysis nor any data were provided therein. We use our new HARPS-TERRA software to derive new RV measurements on GJ 667C (Anglada-Escudé & Butler 2012). The root mean square (rms) of these RVs is 3.89 m s⁻¹, which is significantly larger than the median internal precision (∼1.1 m s⁻¹) and the typical rms found on other stable M dwarfs (Bonfils et al. 2011). We found that the CCF RV measurements provided by the ESO archive were noisier (rms ∼ 4.3 m s⁻¹), causing the signals reported in Section 4 to appear less significant. In order to obtain more secure detections, we obtained 21 new measurements with the Planet Finder Spectrograph (PFS) between 2011 August and 2011 October. PFS is a cross-dispersed echelle spectrograph recently installed at the 6.5 m Clay Magellan Telescope at Las Campanas Observatory, and uses the iodine cell technique to obtain RV measurements at 1–2 m s⁻¹ precision (Crane et al. 2010).

GJ 667C has also been observed with High Resolution Echelle Spectrometer (HIRES)/Keck (Vogt et al. 1994) using the iodine cell method for just over a decade. Such a long time of observation should deliver tighter constraints on long-period signals. However, in 2004 August, the HIRES CCD array was replaced and the data obtained prior to this upgrade are of markedly inferior quality. Post-fix HIRES measurements show similar scatter to the HARPS and PFS RVs, so only 20 post-fix HIRES observations were used in our analysis. We emphasize that the signals discussed here were first detected using HARPS-TERRA measurements only, and that the contribution of the PFS and HIRES data was to improve the sampling cadence and increase the significance of the detections.

According to Skiff (2010), GJ 667C (HR 6426C) is classified as an M1.5 dwarf. The star is a common proper motion companion to the K3V+K5V binary GJ 667AB (HR 6426AB; Kron et al. 1957). At the distance to the system (∼6.8 pc), the minimum physical separation between GJ 667C and GJ 667AB is ∼230 AU. The metallicity of GJ 667AB has been measured before (e.g., Perrin et al. 1988) and amounts to [Fe/H] = −0.59  ± 0.10, meaning that the system is metal-poor compared to the Sun. The same studies show that the GJ 667AB pair is well within the main sequence, indicating an age between 2 and 10 Gyr (e.g., Cayrel de Strobel 1981). The membership of GJ 667 to a Galactic population is unclear. Although its low metallicity points to a thick disk membership, its total velocity in the Galactic local standard of rest is rather low (44.6 ± 1.5 km s⁻¹), which is the characteristic of thin disk kinematics (see Figure 3 in Bensby et al. 2003). The most recent parameters of GJ 667AB can be found in Cvetkovic & Ninkovic (2011) and Tokovinin (2008). Using the empirical relations given by Delfosse et al. (2000), the Hipparcos parallax of the GJ 667AB pair (van Leeuwen 2007) and K-band photometry from the Two Micron All Sky Survey (Skrutskie et al. 2006), we derived a mass of 0.310 ± 0.019 M_☉ for GJ 667C. The luminosity and effective temperature of GJ 667C are derived from the models in Baraffe et al. (1998) by assuming the aforementioned mass, metallicity, and an age of 5 Gyr.

Keplerian orbital fits to the combined RV data were obtained using the SYSTEMIC interface (Meschiari et al. 2009), which allows the interactive least-squares adjustment of complex multiplanetary systems to several data sets. To determine whether there was a significant periodicity remaining in the data, we used a custom-made version of a least-squares periodogram (Cumming 2004) that adjusts a separate zero-point offset to each instrument (HARPS, PFS, and HIRES).

To quantify the significance of a new signal, we estimated its false alarm probability (FAP) empirically. We created 10⁵ synthetic sets by randomly permutating the RV measurements over the same observing epochs (while retaining membership within each instrument). We then computed the periodogram of each synthetic set. A false alarm was identified when a synthetic data set produced a periodogram peak higher than the power of the signal under inspection. The number of false alarms was then divided by the number of simulations to derive the FAP, which was used as a measure of the probability that a spurious detection arose due to an unfortunate arrangement of the noise. This method is widely used to assess the likelihood of periodic signals and a detailed description can be found elsewhere (e.g., Cumming 2004). Since a few tens of M dwarfs have been intensively followed up at 1 m s⁻¹ precision (Bonfils et al. 2011) and to minimize the chances of detecting a false positive, only signals with an FAP <1% were added to the solution (dotted lines in Figure 1).

Zoom In Zoom Out Reset image size

The first detected signal was an extremely significant periodicity at 7.2 days (see Figure 1). No false alarms were found in any of the 10⁵ synthetic sets indicating an FAP < 0.001%. The signal corresponds to a planet with a minimum mass (Msin i) of 5.2 M_⊕ (GJ 667Cb) and a slightly eccentric orbit.

After subtracting a Keplerian solution for GJ 667Cb, a secular trend was the next most significant signal. The magnitude of the trend (∼1.8 m s⁻¹ yr⁻¹) is compatible with the gravitational pull from the GJ 667AB pair (maximum value is ∼3.6 m s⁻¹ yr⁻¹) but could also be caused by an additional unseen long-period companion. The corresponding FAP of this signal is 0.055%, a very significant detection. A tentative solution with a period of 7000 days provides a slight improvement to the fit due to some curvature detected when combining HIRES+PFS measurements (see Figure 2). We estimate that one more year of observations is required to determine if the signal is due to an additional low-mass companion or due to the gravitational pull of the GJ 667AB binary.

Zoom In Zoom Out Reset image size

The next signal also has a very low FAP (0.034%), and implies a planet with a period of 28.15 days and Msin i ∼ 4.5 M_⊕ (GJ 677Cc). Although the period is close to the lunar aliasing frequency (∼ 27.3 days), the orbital phase coverage is complete thanks to the multi-year time span of the observations. Because of its small amplitude, eccentric solutions cannot be ruled out (Shen & Turner 2008; O'Toole et al. 2009). A Monte Carlo Markov Chain (MCMC) analysis (Ford 2005) indicates that this eccentricity (with 98% confidence) must be less than 0.27. At a semimajor axis of 0.123 AU, the stellar flux S reaching the top of its atmosphere is 90% of the solar flux received by Earth (S₀). Using L and T_eff for the host star, the relations given in Kane & Gelino (2011) provide boundaries of the zone at which liquid water could exist on an Earth-like planet (also called liquid water habitable zone, or HZ). In the canonical model presented in Kasting et al. (1993) and updated in Selsis et al. (2007), the inner and outer boundaries depend on the fractional cloud coverage for the putative planet and are displayed as thick gray rectangles in Figure 3. Even though these limits are uncertain, GJ 667Cc comfortably falls in this HZ and also satisfies the empirical limits set by an unhabitable Venus and a possibly habitable early Mars (Selsis et al. 2007). Let us remark that the ultimate capability of GJ 667Cc to support liquid water depends on properties that are not yet known (e.g., albedo, atmospheric composition, and interior dynamics). Detailed studies using realistic climatic and geodynamical models (e.g., Heng & Vogt 2011; Wordsworth et al. 2011) are needed to better assess its chances of supporting life.

Zoom In Zoom Out Reset image size

After fitting for GJ 667Cc, we found a group of candidate periodicities between 75 and 105 days. This time domain is strongly affected by aliases. For example, HARPS data alone favor a period of 91 days, which is uncomfortably close to a very clear signal detected in two activity indices (Section 5). When combining all the data, the 75 day periodicity had the lowest FAP (0.021%), indicating that such signal could not be ignored in the analysis. A tentative orbital solution assuming a circular orbit is given in Table 1. A putative super-Earth with P ∼ 75 days could also support liquid water if its atmosphere contained high concentrations of CO₂ (e.g., GJ 581d in Wordsworth et al. 2011). We reiterate that, until more observations become available, this signal should be considered with due caution.

In a more thorough analysis (not presented here for brevity), we examined other orbital solutions with up to four signals at alternative periods (including, among others, signals at periods of 105, 120, and 33 days). All these attempts delivered significantly poorer fits, extreme eccentricities, and planetary systems that were unstable on timescales shorter than 1 Myr.

We also performed long-term N-body simulations based on some of our fits using the Hybrid simplectic integrator included in the Mercury integration package (Chambers 1999). We included the first-order partial post-Newtonian correction in the central star's gravitational potential as in Lissauer & Rivera (2001) and used a time step of 0.2 days. We assumed a coplanar system throughout. Since the ratio of the periods of the inner two planets is near 4:1, we checked the four critical angles (involving only the periastron longitudes, ω) for the 4:1 mean motion resonance (MMR). We also examined the difference in the periastron longitudes. We found that the critical angles for the MMR circulate, whereas the difference in the longitudes of periastron librates about 180 deg with amplitude ∼90 deg. Results of our simulations also show that the eccentricities of the planets librate between 0.0 and 0.235 for companion b and between 0.04 and 0.265 for companion c in opposite senses as a result of angular momentum conservation. Thus, the system appears to be protected by a secular resonance between the two inner planets in which (1) the orbits can become nearly anti-aligned when the eccentricity of b is small, and (2) the periastron longitudes are nearly perpendicular when the eccentricity of b is large. Most of the time, the system is in the second configuration (in good agreement with the fitted parameters in Table 1). This stabilizing mechanism appears to function for at least the first 25 Myr of our simulations. Further research into the dynamical evolution of this system is warranted once the nature of the other two signals (long-period trend and the 75/91 day candidate) is better understood.

Here, we show the analysis of the time series of three activity indicators: the Bisector span (BIS), the full width at half-maximum (FWHM) of the CCF, and the Ca ii H+K S-index in the Mount Wilson system (S-index). Measurements of the BIS and the FWHM were provided by the HARPS-ESO pipeline. The S-index (Baliunas et al. 1995) was directly measured on the blaze-corrected spectra using the definitions given by Lovis et al. (2011). Since the BIS and FWHM could not be obtained in the iodine cell approach, we limited our analysis to the HARPS observations only. Briefly, the BIS is a measure of the asymmetry of the average spectral line and should correlate with the RV if the observed offsets are caused by spots or plages rotating with the star (Queloz et al. 2001). The FWHM is a measure of the width of the mean spectral line and its variability is usually associated with changes in the convective patterns on the surface of a star. Lovis et al. (2011) found that the FWHM is the most informative index when investigating the correlation of stellar activity with RV variability on cool stars (T_eff < 4600 K). The S-index is an indirect measurement of the chromospheric emission which depends upon the intensity of the stellar magnetic field. Because the strength of the magnetic field affects the efficiency of convection, the S-index could also correlate with observed RV variability. Since the connection between activity and RV jitter on M dwarfs is not well understood at the few m s⁻¹ level (Lovis et al. 2011), we limit our analysis to evaluate if any of the indices has a periodicity that could be related to any of the detected candidates.

While the BIS did not show any significant periodicity, the S-index and especially the FWHM did show significant power around 105 days (see periodograms in Figure 4). To obtain meaningful periodograms, one outlier point had to be removed in the FWHM (a 9.4σ outlier at JD = 2,454,677.66) and another one from the S-index (a 13σ outlier at JD = 2,454,234.79). The FAPs of such signals were obtained by applying same empirical method used for the RVs. We found that the 105 day peak in the S-index had an FAP of 0.2% and no false alarms were found in any of the synthetic sets generated for the FWHM (FAP<0.001%). Given the low metallicity and the age estimates of the GJ 667AB binary, GJ 667C should have a slow rotation rate (Irwin et al. 2011) compatible with the 105 day signal observed in these two indices. Because of the significant aliases affecting the ∼100 day time domain, it is likely (but not conclusive) that GJ 667Cd is a spurious signal resulting from localized magnetic activity rotating with the stellar surface. In addition to more RV observations to improve the phase sampling, a photometric follow-up could help to reveal the true nature of GJ 667Cd.

Zoom In Zoom Out Reset image size

We have derived precision RV measurements from public HARPS spectra using a least-squares matching approach on the M dwarf GJ 667C and thereby detected the Doppler signature of (at least) two planets. Additional observations with PFS and HIRES confirm the detection of these signals and further constrain the orbital parameters. Even though the public CCF Doppler measurements are not as precise, the CCF method still provides useful information on stellar activity that can be used to investigate the origin of candidate signals.

GJ 667Cc is the super-Earth candidate most securely detected within the liquid water HZ of another star. As for other proposed candidates (e.g., GJ 581d and HD 85512b announced in Mayor et al. 2009; Pepe et al. 2011, respectively), its actual capability of supporting liquid water depends on many physical properties that are yet unknown. Using the relations given by Charbonneau et al. (2007), the reported candidates have non-negligible probabilities of transiting in front of the star (∼2.7%, 1.1%, and 0.6% for planets b, c, and d, respectively). Even though these probabilities are low, the estimated transit depth assuming a density similar to Earth, is about 0.3%, which can be measured using small aperture telescopes. Statistical extrapolations based on Doppler, transit, and microlensing surveys indicate that such planets should be abundant around main-sequence stars (Mayor et al. 2009; Borucki et al. 2011; Cassan et al. 2012). With the new generation of optical and infrared spectrographs, many nearby M dwarfs will be efficiently surveyed for low-mass planets. If the detection rate holds, very soon we may have a real chance of searching for spectroscopic signatures of water and life on one of these worlds.

We thank the constructive comments given by the anonymous referee. This research has been partially supported by the Carnegie Postdoctoral Program to G.A.E.; ICM P07-021-F, FONDAP 15010003, and BASAL-CATA PFB-06 grants to D.M. and P.A.; NSF grant AST-0307493 to S.S.V.; NASA NNX07AR40G and NASA Keck PI program grants to R.P.B.; NASA grants NNA04CC08A and NNX09AN05G to N.H.; ARC grant DP774000 to C.G.T.; and Fondecyt grant 3110004 to J.S.J. This work is based on data obtained from the ESO Science Archive Facility. Observations were obtained from Las Campanas Observatory and W. M. Keck Observatory. W. M. Keck Observatory is operated jointly by Univ. of California and California Institute of Technology. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.