© ESO, 2015

1. Introduction

As the search for new exoplanets with the radial velocity (RV) technique continues, our detection capabilities have increased dramatically in two dimensions: smaller signals, such as those induced by short-orbit terrestrial-mass planets, can now be found thanks to the increased precision of modern RV spectrographs, and longer-period signals can also be identified thanks to an ever-increasing time baseline of data collection. Both of these detection capabilities are critical if we wish to find planetary systems comparable to our own solar system.

At the cutting edge of long-period planet detection is the search for Jupiter analogs. In terms of signal size, Jupiter induces by far the largest RV shift on the Sun of any planet and therefore should be the first planet detected in a twin of the solar system. However, the need for more than a decade of observational coverage with a consistently precise RV zero point complicates the process of discovering such a planet, as does the possibility of stellar magnetic activity cycles inducing RV variations on a similar timescale (Saar & Donahue 1997).

In spite of these challenges, a number of Jupiter analog candidates have been discovered from long-duration planet surveys (see, e.g., Vogt et al. 2000; Naef et al. 2003; Wright et al. 2007, 2008; Howard et al. 2010; Boisse et al. 2012; Marmier et al. 2013; Wittenmyer et al. 2014; Feng et al. 2015). These surveys also demonstrate the prevalence of solar-like magnetic activity cycles, which can manifest in the stellar RVs (Lovis et al. 2011; Zechmeister et al. 2013). These cycles can generally be traced by their influence on activity diagnostics such as the Ca II H and K lines (Baliunas et al. 1995) and asymmetries of the spectral cross-correlation function (Dravins 1985). With careful attention to these diagnostics, we can minimize the influence of activity cycles and recover planet signatures around cycling stars (Dumusque et al. 2011).

Over the last four years, we have been conducting an RV planet search with the HARPS spectrograph (Mayor et al. 2003) targeting a specific subset of 63 stars selected by the proximity of their fundamental parameters to those of the Sun (T_eff within ~100 K, log (g) within ~0.1 dex, and [Fe/H] within ~0.1 dex). These solar twin stars offer a unique opportunity to determine their compositions at a precision as good as 2% (see, e.g., Meléndez et al. 2009; Ramírez et al. 2009; Bedell et al. 2014). We aim to combine this high-precision chemical analysis with the m s^-1-level detection capabilities of HARPS to investigate the star-planet connection at an unprecedented level of detail.

In this paper we present our first detection from the planet search program: a Jupiter twin around the solar twin HIP11915. HIP11915 was previously established as a solar twin in Ramírez et al. (2014) based on its fundamental properties (T_eff = 5760 ± 4 K, log g = 4.46 ± 0.01, and [Fe/H] = −0.059 ± 0.004). Its age, estimated from isochrone fitting to be 4.0 ± 0.6 Gyr, is also in keeping with the solar value (Ramírez et al. 2014). HIP11915 is a bright (m_V = 8.6) star with relatively low activity (), making it an excellent target for high-precision RV monitoring. Based on a Keplerian signal in our HARPS RVs, we determine the planet candidate’s orbital parameters. We also investigate the influence of stellar activity on our RV measurements.

2. Data

All of the RV data used in this paper were obtained with the HARPS spectrograph on the ESO 3.6 m telescope at La Silla Observatory and were processed with the dedicated HARPS pipeline. Forty-three data points spanning a time range from October 2009 to January 2015 came from our HARPS large program dedicated to solar twin stars (program ID 188.C-0265). An additional 13 measurements were taken from the publicly available archive data and span a time range from October 2003 to November 2013 (program IDs 072.C-0488, 089.C-0732, 091.C-0034, 092.C-0721, and 183.C-0972).

For the 43 measurements taken through our program, the star was observed in high-accuracy mode with the ThAr lamp as a simultaneous reference source. Each observation used a 900 s exposure time to minimize potential RV noise from the five-minute p-mode oscillations of solar-type stars. To obtain a RV measurement, we selected a binary mask based on the G2 spectral type for cross-correlation by the pipeline. Owing to their generally high signal-to-noise ratios (median S/N over the spectral range on the order of 100), the pipeline-produced RV error estimates were generally below 1 m s^-1. We expect the actual uncertainties to be somewhat higher owing to stellar noise and potential undiagnosed instrumental effects. To account for these factors, we added the pipeline-estimated errors in quadrature with a baseline error level of 1 m s^-1, following Dumusque et al. (2011) and in accordance with the scatter in our own HARPS survey’s quietest stars.

The 13 archive measurements have exposure times ranging from 120 to 900 s. All were taken in high-accuracy mode without a source on the simultaneous reference fiber, so that even the low S/N measurements are not in danger of spectral contamination by the simultaneous reference source. The data were processed by the pipeline and we modified the error estimates as described above. An additional 1.1 m s^-1 was added in quadrature to the errors on all exposures under 15 min, based on the expected p-mode signal for a solar twin (Bazot et al. 2011).

We analyzed all of the reduced HARPS spectra to measure the Ca II H and K activity indices S_HK and ) at each observation epoch. This was done according to the methods of Lovis et al. (2011). Errors on S_HK were estimated from photon noise.

Complete RV and S_HK data are available in Table 3, along with the HARPS pipeline-estimated activity indicators bisector inverse span (BIS) and full width at half maximum of the cross-correlation function (FWHM).

3. Analysis and results

The RMS of the RVs is 6.5 m s^-1, indicating the presence of a signal stronger than the predicted stellar jitter of 2.1 m s^-1 (Wright 2005). Visual inspection and a periodogram analysis of the data clearly show a long-period trend in the stellar RVs. We ran a grid search across a range of periods from 0 to 10 000 days and found that the minimum chi-squared value lies near 3600 days.

We inspected the data for signs of stellar activity manifesting in the RVs. A strong correlation between RV and the activity index S_HK is seen in the measurements from recent years, although there is no correlation in the archival data. Similar RV-correlated behavior is seen for other activity indicators including the bisector inverse slope (BIS), FWHM, V_span, and the biGaussian fit differential (as defined by Queloz et al. 2001; Boisse et al. 2011; Figueira et al. 2013). Doing a weighted linear least-squares fit and removing the trend with any of these activity indicators from the RVs results in reducing the correlations while preserving the long-period planet signal with high significance in a periodogram analysis. We additionally attempted the prescribed method of Dumusque et al. (2011) of fitting a long-period Keplerian first to an activity diagnostic and then (with period and phase fixed) the RVs. This method failed to remove the long-period power from the RVs regardless of whether S_HK, BIS, or FWHM was used as the diagnostic. The failure of these methods to successfully “clean” the data of its long-period signal suggests that the signal originates from a true Doppler shift rather than stellar activity, although activity is also present in the data.

3.1. Best model: Keplerian signal + S_HK correlation

In the results presented here, we account for stellar activity by fitting a linear trend between RV and S_HK simultaneously with the Keplerian signal. This way the strength of the correlation between activity indicator and RV is free to change as the scatter induced by the Keplerian signal is removed. We use S_HK as the primary tracer of activity owing to the unreliability of the BIS diagnostic for slow-rotating stars and the demonstrated efficiency of S_HK in removing activity cycle signals (Boisse et al. 2011; Meunier & Lagrange 2013). Other potential diagnostics and alternative model choices are explored in Sect. 3.2.

We ran a Markov chain Monte Carlo (MCMC) analysis to fit the RV data to a model consisting of three components: a Keplerian signal, a linear correlation term with the S_HK values, and an offset term. This gives seven free parameters: the Keplerian parameters (P, K, e, ω + M₀, ω - M₀), where P is the orbital period, K is the RV semi-amplitude, e is the eccentricity, ω is the argument of periastron, and M₀ is the mean anomaly at reference date t₀ = 2 457 000.0; a slope parameter α which relates the activity-induced change in RV to the corresponding change in S_HK; and an offset parameter, C. We include an eighth free parameter, σ_J, representing the jitter that must be added in quadrature with the estimated RV error to make the measurement uncertainties consistent with the scatter in the residuals. This jitter physically corresponds to unknown Gaussian error that may arise from stellar noise or undiagnosed instrumental variations. When inputting the data, we subtracted a sigma-clipped median value from RV, S_HK, BIS, and FWHM for simplicity. This means that the offset parameter C is partially set by the difference between the median RV of the data set and the systemic RV of the star, but also includes the RV offset in the linear fit with the median-subtracted S_HK values (or any other correlations used in the model; see Sect. 3.2).

All of the fit parameters were given uniform priors except for P and K, which were sampled in log-space with a log-uniform prior. We ran five parallel chains and used the requirement of a Gelman-Rubin test statistic below 1.01 for all parameters to avoid non-convergence (Gelman & Rubin 1992). The best-fit orbital solution, plotted in Fig. 1, corresponds to an approximately Jupiter-mass planet on a 10-year orbit with eccentricity consistent with zero (see Table 1). As seen in the MCMC posterior distributions (Fig. 5), the chief source of uncertainty in the period measurement is in its degeneracy with eccentricity and the associated argument of periapse, parameters that are inevitably poorly constrained when only one orbit has been observed. The best-fit jitter term is consistent with the predicted stellar jitter of 2.1 m s^-1 from the formula of Wright (2005).

After subtracting the best-fit planet signal, we searched for additional signals in the residuals with a generalized Lomb-Scargle periodogram analysis (Zechmeister & Kürster 2009). No significant periodicities were identified in the residuals to the planet fit or in the residuals to the planet and the fitted S_HK trend. We place constraints on the presence of additional planets in the system using a test based on the bootstrapping and injection of a planet signal into the residuals. For this test, we use the residuals to the planet fit only, and inflate the error bars (with jitter included) such that a flat-line fit to the residuals gives a reduced chi-squared of one. We then run a set of trials in which the RV data are resampled with replacement at the times of observation, a planet signal with period and semi-amplitude fixed and the other orbital parameters randomized is injected, and the reduced chi-squared of a flat-line fit to these simulated data is measured. If the chi-squared statistic has increased with three-sigma significance in 99 out of 100 trials at a given period and velocity semi-amplitude, we consider such a planet signal to be excluded based on its inconsistency with the scatter in the actual RVs. We run this trial until upper limits on planet presence are found for all of the periods in a log-sampled grid. The results rule out the presence of gas giants within a 1000-day orbital period, leaving open the possibility of inner terrestrial planets (Fig. 2).

Fig. 1 Best-fit orbital solution for HIP11915b. Data points are the observed RVs with offset and SHK correlation terms subtracted. One residual point is located outside of the plot range, as denoted by the arrow.

Fig. 2 Exclusion limits on potential planets in the RV data after subtraction of the Keplerian signal. Any signal in the parameter space located above the lines plotted would have induced a larger scatter in the RV residuals with three-sigma significance if present. Potential trends with activity indicators were not removed for this test to avoid underestimating the RV scatter.

3.2. Alternative model choices

Fig. 3 Correlations for the activity tracers SHK, BIS, and FWHM with RV. The upper panels shows the observed RVs. The lower panels show the residuals to the 3800-day signal fit (as presented in Sect. 3.1). Dashed lines are linear least-squares best fits to the data.

It is not obvious from theory that the Keplerian + S_HK correlation model used above should be the best model choice. Other activity indicators exist and have been used with success to remove activity trends from RVs in the past. Chief among these are BIS and FWHM, both characteristics derived from the RV cross-correlation function (CCF) which measure the average spectral line asymmetry and/or broadening at the epoch of the observation (see, e.g., Queloz et al. 2001, 2009). Either BIS or FWHM could potentially be used in place of S_HK. Also, because these indicators arise from distortions to the rotationally broadened stellar absorption lines while S_HK probes chromospheric emission, we could combine S_HK with a CCF-derived measurement to simultaneously track two separate physical indicators of activity.

We considered several alternative models before settling on the Keplerian + S_HK model as the best choice. A summary of the models tried is given in Table 2. In all cases, the fitting was done with an MCMC as described above. The inclusion of a Keplerian signal is strongly justified, with the Keplerian + S_HK fit having a false alarm probability on the order of 10^-15 compared to the best fitting activity model (S_HK + FWHM). While the lowest χ² value is obtained for the fit to a Keplerian, S_HK, and FWHM, the addition of the FWHM term is not fully justified by the reduction in scatter around the fit. An F-test gives a formal false alarm probability of 13% for the Keplerian + S_HK + FWHM model compared to the simpler Keplerian + S_HK model. The Keplerian + S_HK fit has a false alarm probability of 0.3% when compared to the Keplerian only model. Using BIS in place of S_HK technically gives a slightly better fit, but we chose to use S_HK in the main analysis because of BIS’s previously noted unreliability for slow rotators. Regardless of the choice of activity indicator(s), all Keplerian fit parameters remain within the uncertainty intervals quoted in Table 1.

4. Planet or activity cycle?

The question of the role of stellar activity role in the observed RVs is of special importance given the proximity of the signal period to the 11-year solar magnetic cycle. The signal size is also consistent with the theoretically predicted 10–11 m s^-1 signal for a solar activity cycle (Meunier & Lagrange 2013). However, an activity cycle should manifest not only in the observed RVs, but also in other diagnostics. We expect a positive correlation in RV with the activity index S_HK as well as with other diagnostics such as BIS and the FWHM of the cross-correlation function (Lovis et al. 2011; Zechmeister et al. 2013). While such correlations are seen, they appear insufficient to fully explain the observed RVs.

For at least two other solar twins in our sample, similar Jupiter-like RV variations are seen; however, in these cases subtracting a linear fit with S_HK entirely removes the long-period power from a periodogram and reduces the RMS of RVs to a level consistent with the expected stellar jitter. In contrast, the RV signal inferred from the variation of S_HK or any other activity tracer for HIP11915 is much smaller in size than the total RV signal, so that significant power remains at long periods in the periodogram.

A large sample of stars with activity cycles have been studied with HARPS by Lovis et al. (2011), who derived empirical predictions for the expected correlation between RV and other activity indictors with as a function of stellar temperature and metallicity. The correlations seen for HIP11915 are generally lower than the predicted values: C_RV = 30 ± 1 m s^-1 (unit )^-1, C_FWHM = 58 ± 2 m s^-1 (unit )^-1, and C_BIS = 22.4 ± 0.7 m s^-1 (unit )^-1 from their Eqs. (9)–(12), while the C values recovered from our data are 11, 45, and 16 respectively. This disagreement suggests that the strength of activity seen in the data is insufficient to explain the RV signal.

Fig. 4 Evolution of the activity indices SHK, BIS, and FWHM with time. Errors on BIS and FWHM are approximated as the errors on the RVs.

An alternative interpretation is that we are observing activity overlaid with a planet signal in the RVs. In this view, the correlation with activity indicator(s) in the RV residuals is independent of the Keplerian signal. One piece of evidence in favor of this interpretation comes from inspecting the correlation plots before and after removal of the 3800-day Keplerian signal. If this signal is truly unrelated to the stellar activity, it should add random scatter about the correlation. Indeed, subtracting the Keplerian signal from the data reduces the scatter about the best-fit trend between RV and S_HK from 7.7 to 3.0 m s^-1 and has a comparable effect on the correlations with other activity indices (Fig. 3).

Inspecting the behavior of S_HK with time sheds some light on the evolution of HIP11915’s activity level. As depicted in Fig. 4, the stellar activity level exhibited an incoherent scatter until mid-2013, at which point it increased from the equivalent of to ~–4.82 by mid-2014. A concurrent rise in BIS and FWHM is also seen at this epoch. This level of variability is consistent with the behavior of the Sun, which has a mean ) of −4.94 and an amplitude of variation around 0.1 dex over a full activity cycle (Hall et al. 2007). Given the recent increase in Ca emission, we could then expect a corresponding increase in short-duration activity phenomena such as spots and plages to manifest in the more recent RVs. This would add power to the correlations with line asymmetry indicators.

Ideally, we would be able to remove activity phenomena such as spots and plages by averaging our measurements over the timescale of the stellar rotation period, on the order of a month. Unfortunately the sampling of the data is too uneven to average effectively over the rotation period. We average instead over an entire observing season, which is a regime in which the time coverage should be distributed evenly enough to remove activity signals at periods shorter than a cycle. The resulting data from this seasonal binning trace the 3800-day signal well, but show no statistically significant correlations with activity indicators.

In short, we have considered a variety of activity diagnostics and methods of cycle removal and the 3800-day signal remains robustly present in the data. The coincidence in phase between the rise in S_HK and the 3800-day signal maximum, however, calls for further monitoring of the star before we can definitively rule out the possibility of an activity cycle inducing the RV signal.

5. Summary

We have detected a Keplerian signal corresponding to a gas giant planet with a 10-year orbital period, consistent with a Jupiter twin, around the solar twin HIP11915. While we cannot conclusively exclude an activity cycle interpretation, commonly employed activity tracers fail to account for the signal. If the 3800-day signal is entirely induced by an activity cycle, the failure of various techniques to remove it bodes poorly for the prospect of correcting for such cycles in future long-period planet searches.

If this signal is truly planetary in origin, the HIP11915 system is a close analog to the solar system. Our analysis shows the planet HIP11915b to be a close match to Jupiter both in mass and orbital period, and its host star is extraordinarily similar to the Sun. Beyond having fundamental physical properties close to the Sun, initial spectroscopic analysis suggests that HIP11915 is also a solar twin in the sense that its detailed abundance pattern matches the solar pattern (Melendez et al., in prep.). The presence of a Jupiter twin and a solar-like composition both make HIP11915 an excellent prospect for future terrestrial planet searches.

Acknowledgments

We thank Jason Wright and Gaspare Lo Curto for helpful discussions. We are also indebted to the many scientists and engineers whose past work on and observations with HARPS made such a long-period planet discovery possible. M.B. is supported by the National Science Foundation (NSF) Graduate Research Fellowships Program (grant No. DGE-1144082). J.M. would like to acknowledge support from FAPESP (2012/24392- 2) and CNPq (Bolsa de Produtividade). J.B. and M.B. acknowledge support for this work from the NSF (grant No. AST-1313119). J.B. is additionally supported by the Alfred P. Sloan Foundation and the David and Lucile Packard Foundation.

References

Online material

Fig. 5 Constraints on the fit to the RVs from the MCMC analysis. Parameters are defined in Sect. 3. Plot made using code from Foreman-Mackey et al. (2014).

All Tables

All Figures

	Fig. 1 Best-fit orbital solution for HIP11915b. Data points are the observed RVs with offset and SHK correlation terms subtracted. One residual point is located outside of the plot range, as denoted by the arrow.
In the text

	Fig. 2 Exclusion limits on potential planets in the RV data after subtraction of the Keplerian signal. Any signal in the parameter space located above the lines plotted would have induced a larger scatter in the RV residuals with three-sigma significance if present. Potential trends with activity indicators were not removed for this test to avoid underestimating the RV scatter.
In the text

	Fig. 3 Correlations for the activity tracers SHK, BIS, and FWHM with RV. The upper panels shows the observed RVs. The lower panels show the residuals to the 3800-day signal fit (as presented in Sect. 3.1). Dashed lines are linear least-squares best fits to the data.
In the text

	Fig. 4 Evolution of the activity indices SHK, BIS, and FWHM with time. Errors on BIS and FWHM are approximated as the errors on the RVs.
In the text

	Fig. 5 Constraints on the fit to the RVs from the MCMC analysis. Parameters are defined in Sect. 3. Plot made using code from Foreman-Mackey et al. (2014).
In the text