Habitablezone superEarth candidate in a sixplanet system around the K2.5V star HD 40307^{⋆}
^{1} University of Hertfordshire, Centre for Astrophysics Research, Science and Technology Research Institute, College Lane, AL10 9AB, Hatfield, UK
email: m.tuomi@herts.ac.uk
^{2} University of Turku, Tuorla Observatory, Department of Physics and Astronomy, Väisäläntie 20, 21500 Piikkiö, Finland
email: mikko.tuomi@utu.fi
^{3} Universität Göttingen, Institut für Astrophysik, FriedrichHundPlatz 1, 37077 Göttingen, Germany
email: guillem.anglada@gmail.com
^{4} Lohrmann Observatory, Technical University Dresden, 01062 Dresden, Germany
^{5} UCO/Lick Observatory, Department of Astronomy and Astrophysics, University of California at Santa Cruz, Santa Cruz, CA 95064, USA
^{6} Department of Terrestrial Magnetism, Carnegie Institute of Washington, Washington, DC 20015, USA
Received: 21 August 2012
Accepted: 15 October 2012
Context. The K2.5 dwarf HD 40307 has been reported to host three superEarths. The system lacks massive planets and is therefore a potential candidate for having additional lowmass planetary companions.
Aims. We rederive Doppler measurements from public HARPS spectra of HD 40307 to confirm the significance of the reported signals using independent data analysis methods. We also investigate these measurements for additional lowamplitude signals.
Methods. We used Bayesian analysis of our radial velocities to estimate the probability densities of different model parameters. We also estimated the relative probabilities of models with differing numbers of Keplerian signals and verified their significance using periodogram analyses. We investigated the relation of the detected signals with the chromospheric emission of the star. As previously reported for other objects, we found that radial velocity signals correlated with the Sindex are strongly wavelength dependent.
Results. We identify two additional clear signals with periods of 34 and 51 days, both corresponding to planet candidates with minimum masses a few times that of the Earth. An additional sixth candidate is initially found at a period of 320 days. However, this signal correlates strongly with the chromospheric emission from the star and is also strongly wavelength dependent. When analysing the red half of the spectra only, the five putative planetary signals are recovered together with a very significant periodicity at about 200 days. This signal has a similar amplitude as the other new signals reported in the current work and corresponds to a planet candidate with Msini ~ 7 M_{⊕} (HD 40307 g).
Conclusions. We show that Doppler measurements can be filtered for activityinduced signals if enough photons and a sufficient wavelength interval are available. If the signal corresponding to HD 40307 g is a genuine Doppler signal of planetary origin, this candidate planet might be capable of supporting liquid water on its surface according to the current definition of the liquid water habitable zone around a star and is not likely to suffer from tidal locking. Also, at an angular separation of ~46 mas, HD 40307 g would be a primary target for a future spacebased directimaging mission.
Key words: methods: statistical / methods: numerical / techniques: radial velocities / stars: individual: HD 40307
Appendix A is available in electronic form at http://www.aanda.org
© ESO, 2012
1. Introduction
Current highprecision spectrographs, such as the High Accuracy Radial Velocity Planet Searcher (HARPS; Mayor et al. 2003) and the High Resolution Echelle Spectrograph (HIRES; Vogt et al. 1994), enable detections of lowmass planets orbiting nearby stars. During recent years, radial velocity (RV) planet searches have revealed several systems of superEarths and/or Neptunemass planets around nearby stars (e.g. Mayor et al. 2009a,b; Lovis et al. 2011; Pepe et al. 2011; Tuomi 2012).
The system of three superEarths orbiting HD 40307 has received much attention because the planets appear in dynamically packed orbits close to mean motion resonances (Mayor et al. 2009a). This has been used as an argument to suggest that lowmass planets may be found in highly compact multiple systems that are still stable in longterm, e.g. a possibility of having ten planets with masses of 17 M_{⊕} within a distance of 0.26 AU on stable orbits (Funk et al. 2010). However, the physical nature of these companions as scaledup versions of the Earths is not entirely clear (Barnes et al. 2009a). Their masses, between those of Earth and Neptune, suggest they are Neptunelike protogas giants that could not accumulate enough gas before it was blown away by the newly born star. On the other hand, recent transit observations of hot superEarths around bright nearby stars (Léger et al. 2009; Batalha et al. 2011; Winn et al. 2011) indicate that a good fraction of these hot superEarth mass objects can have rocky compositions.
In this article we reanalyse the 345 HARPS spectra publicly available through the ESO archive using a newly developed software tool called HARPSTERRA (templateenhanced radial velocity reanalysis application; AngladaEscudé & Butler 2012). Instead of the classic crosscorrelation function method (CCF) implemented by the standard HARPSESO data reduction software (HARPSDRS), we derive Doppler measurements by leastsquares matching of each observed spectrum to a high signaltonoise ratio template built from the same observations. A description of the method and the implementation details are given in AngladaEscudé & Butler (2012). In addition to an increase in precision (especially for K and M dwarfs), this method allows us to perform additional analyses and tests beyond those enabled by the CCF data products provided by the HARPSDRS. As an example, it allows us to reobtain the RV measurements using only a restricted wavelength range. As we show in Sect. 5, this capability can be instrumental in ruling out the planetary nature of prominent signals correlated with stellar activity.
We rely on the Bayesian framework when estimating the orbital parameters supported by the data, determining the significances of the signals, and the modelling of the noise in the measurements. In previous studies, radial velocities received within an interval of an hour or so have been commonly binned together in an attempt to reduce the noise caused by stellar surfacerelated effects, i.e. stellar oscillations and granulation, and other factors within this timescale (Dumusque et al. 2010). In principle, this would enable the detections of planets smaller than roughly 5 M_{⊕} with HARPS over a variety of orbital distances, even at or near the stellar habitable zone (Dumusque et al. 2010, 2011a). In our approach, and instead of binning, we apply a selfconsistent scheme to account for and quantify correlated noise in the Bayesian framework and use Bayesian model probabilities to show that a solution containing up to six planets is clearly favoured by the data, especially when the redmost part of the stellar spectrum is used in the RV analysis. Only the confluence of refinements in these data analysis methods (reanalysis of the spectra and Bayesian inference) allows the detection and verification of these lowamplitude signals.
We start with a brief description of the stellar properties of HD 40307 (Sect. 2) and describe the statistical modelling of the observations and the data analysis techniques we used (Sect. 3). In Sect. 4 we describe the properties of the RV measurements and perform a detailed Bayesian analysis that identifies up to three new candidate signals. We discuss the stellar activity indicators and their possible correlations with the RV signals in Sect. 5. In this same section, we find that one of the candidates is spuriously induced by stellar activity by showing that the corresponding periodic signal (P ~ 320 days) is strongly wavelength dependent. When the RVs obtained on the redmost part of the spectrum are analysed (Sect. 6), the 320day signal is replaced by a signal of a superEarthmass candidate with a 200day period with a minimum mass of about ~7 M_{⊕} orbiting within the liquid water habitable zone of HD 40307. The analysis of the dynamical stability of the system (Sect. 7) shows that stable solutions compatible with the data are feasible and the potential habitability of the candidate at 200 days (HD 40307 g) is discussed in Sect. 8. We give some concluding remarks and discuss the prospects of future work in Sect. 9.
2. Stellar properties of HD 40307
We list the basic stellar properties of HD 40307 in Table 1. This K2.5 V star is a nearby dwarf with a Hipparcos parallax of 77.95 ± 0.53 mas, which implies a distance of 12.83 ± 0.09 pc. It is somewhat smaller (M_{ ⋆ } = 0.77 ± 0.05 M_{⊕}; Sousa et al. 2008) and less luminous (log L_{star} / L_{⊙} = −0.639 ± 0.060; Ghezzi et al. 2010) than the Sun. The star is quiescent (; Mayor et al. 2009a) and relatively metalpoor with [Fe/H] = −0.31 ± 0.03 (Sousa et al. 2008). It also lacks massive planetary companions, which makes it an ideal target for highprecision RV surveys aiming at finding lowmass planets. According to the calibration of Barnes (2007), HD 40307 likely has an age similar to that of the Sun (~4.5 Gyr).
Stellar properties of HD 40307.
3. Statistical analyses
3.1. Statistical models
We modelled the HARPS RVs using a statistical model with a moving average (MA) term and two additional Gaussian white noise components consisting of two independent random variables.
The choice of an MA approach instead of binning is based on the results of Tuomi et al. (2012b) and accounts for the fact that uncertainties of subsequent measurements likely correlate with one another at timescales of an hour in an unknown manner. We limit our analysis to MA models of third order (MA(3) models) because higher order choices did not improve the noise model significantly. Effectively, the MA(3) component in our noise model corresponds to binning. However, unlike when binning measurements and artificially decreasing the size of the data set, this approach better preserves information on possible signals in the data. The two Gaussian components of the noise model are the estimated instrument noise with zeromean and known variance (nominal uncertainties in the RVs) and another with zeromean but unknown variance corresponding to all excess noise in the data. The latter contains the white noise component of the stellar surface, usually referred to as stellar “jitter”, and any additional instrumental systematic effects not accounted for in the nominal uncertainties. Keplerian signals and white noise component were modelled as in Tuomi et al. (2011).
In mathematical terms, an MA(p) model is implemented on measurement m_{i} as (1)where r_{k}(t_{i}) is the superposition of k Keplerian signals at epoch t_{i} and γ is the reference velocity. The random variable ϵ_{i} is the Gaussian white noise component of the noise model with zeromean and variance , where is the (fixed) nominal uncertainty of the ith measurements and is a free parameter describing the magnitude of the jitter component. Finally, the free parameters of the MA(p) model are denoted as φ_{j},j = 1,...,p – they describe the amount of correlation between the noise of the ith and i − jth measurements. The exponential term in Eq. (1) ensures that correlations in the noise are modelled on the correct timescale. Specifically, using hours as units of time, the exponential term (that is always < 1 because t_{i} > t_{i − j}) vanishes in few hours.
To demonstrate the impact of binning on this data set, Sect. 4.1 shows the analyses of binned RVs with the common assumption that the excess noise is purely Gaussian. This corresponds to using the nightly average as the individual RV measurements and setting φ_{j} = 0 for all j in Eq. (1).
3.2. Bayesian analyses and detection thresholds
To estimate the model parameters and, especially, their uncertainties as reliably as possible, we drew random samples from the parameter posterior densities using posterior sampling algorithms. We used the adaptive Metropolis algorithm of Haario et al. (2001) because it can be used to receive robust samples from the posterior density of the parameter vector when applied to models with multiple Keplerian signals (e.g. Tuomi et al. 2011; Tuomi 2012). This algorithm is simply a modified version of the famous MetropolisHastings Markov chain Monte Carlo (MCMC) algorithm (Metropolis et al. 1953; Hastings 1970), which adapts the proposal density to the information gathered from the posterior density. We performed samplings of models with k = 0,1,...,7 Keplerian signals.
The samples from the posterior densities were then used to perform comparisons of the different models. We used the oneblock MetropolisHastings (OBMH) method (Chib & Jeliazkov 2001; Clyde et al. 2007) to calculate the relative posterior probabilities of models with differing numbers of Keplerian signals (e.g. Tuomi & Kotiranta 2009; Tuomi 2011; Tuomi et al. 2011; Tuomi 2012). We performed several samplings using different initial values of the parameter vector and calculated the means and the corresponding deviations as measures of uncertainties of our Bayesian evidence numbers P(m  ℳ_{k}), where m is the measurement vector and ℳ_{k} denotes the model with kKeplerian signals.
The prior probability densities in our analyses were essentially uniform densities. As in Tuomi (2012), we adopted the priors of the RV amplitude π(K_{i}) = U(0,a_{RV}), reference velocity π(γ) = U( − a_{RV},a_{RV}), and jitter π(σ_{J}) = U(0,a_{RV}), where U(a,b) denotes a uniform density in the interval [a,b] . Since the observed peaktopeak difference in the raw RVs is lower than 10 m s^{1}, the hyperparameter a_{RV} was conservatively selected to have a value of 20 m s^{1}. The priors of the longitude of pericentre (ω) and the mean anomaly (M_{0}) were set to U(0,2π), in accordance with the choice of Ford & Gregory (2007) and Tuomi (2012). We used the logarithm of the orbital period as a parameter of our model because, unlike the period as such, it is a scaleinvariant parameter. The prior of this parameter was set uniform such that the two cutoff periodicities were T_{min} and T_{max}. These hyperparameters were selected as T_{min} = 1.0 days and T_{max} = 10T_{obs} because we did not expect to find signals with periods less than 1 day. Also, we did not limit the period space to the length of the baseline of the HARPS time series (T_{obs}), because signals in excess of that can be detected in RV data (Tuomi et al. 2009) and because there might be longperiod signals apparent as a trend with or without curvature in the data set.
Unlike in traditional Bayesian analyses of RV data, we did not use uniform prior densities for the orbital eccentricities. Instead, we used a semiGaussian as with the corresponding normalisation, where the hyperparameter σ_{e} was chosen to have a value of 0.3. This value decreases the posterior probabilities of very high eccentricities in practice, but still enables them if they explain the data better than lower ones (Tuomi 2012; Tuomi et al. 2012a).
For the MA components φ_{j}, we selected uniform priors as π(φ_{j}) = U( − 1,1), for all j = 1,2,3. This choice was made to ensure that the MA model was stationary, i.e. timeshift invariant – a condition that is satisfied exactly when the values are in the interval [–1, 1].
Finally, we did not use equal prior probabilities for the models with differing numbers of Keplerian signals. Instead, following Tuomi (2012), we set them as P(ℳ_{k}) = 2P(ℳ_{k + 1}), which means that the model with k Keplerian signals was always twice as probable prior to the analyses than the model with k + 1 Keplerian signals. While this choice makes our results more robust in the sense that a posterior probability that exceeds our detection threshold is actually already underestimated with respect to equal prior probabilities, there is a physical motivation as well. We expect that the dynamical interactions of planets in any given system make the existence of an additional planet less probable because there are fewer dynamically stable orbits. This also justifies the qualitative form of our prior probabilities for the eccentricities.
Our criterion for a positive detection of kKeplerian signals is as follows. First, we require that the posterior probability of a kKeplerian model is at least 150 times greater than that of the k − 1Keplerian model (Kass & Raftery 1995; Tuomi 2011, 2012; Tuomi et al. 2011; Feroz et al. 2011). Second, we require that the radial velocity amplitudes of all signals are statistically significantly different from zero. Third, we also require that the periods of all signals are wellconstrained from above and below. These criteria were also applied in Tuomi (2012).
We describe the parameter posterior densities using three numbers, namely, the maximum a posteriori (MAP) estimate and the limits of the corresponding 99% Bayesian credibility sets (BCSs) or intervals in one dimension (e.g. Tuomi & Kotiranta 2009).
3.3. Periodogram analysis
As is traditionally the case when searching for periodic signals in time series, we used leastsquares periodograms (Lomb 1976; Scargle 1982) to probe the next most significant periods left in the data. In particular, we used the leastsquares periodograms described in Cumming (2004), which adjust for a sine wave and an offset at each test period and plot each test period against the Fratio statistic (or power) of the fit. While strong powers likely indicate the existence of a periodic signal (though strong powers may be caused by samplingrelated features in the data as well), the lack of them does not necessarily mean that there are no significant periodicities left (e.g. Tuomi 2012). This is especially so in multiKeplerian fits due to strong correlations and aliases between clearly detected signals and yetundetected lower amplitude companions (e.g. AngladaEscudé et al. 2010). The reason is that residuals must necessarily be calculated with respect to a model that is assumed to be correct, which is clearly not the case when adding additional degrees of freedom, i.e. additional planetary signals, to the model. Therefore, determining the reliability of a new detection based on goodnessoffit comparisons is prone to biases, which effectively reduces the sensitivity and reliability of these detections (Tuomi 2012).
While periodograms of the residuals are very useful, they do not properly quantify the significance of the possibly remaining periodicities and, therefore, we used them as a secondary rather than a primary tool to assess the significance of new signals. The analytic falsealarm probability (FAP) thresholds as derived by Cumming (2004) are provided in the figures as a reference and for illustrative purposes only. We used the same periodogram tools to assess the presence of periodicities in the time series of a few activity indicators.
4. Analysis of the RV data
The 345 measurements taken on 135 separate nights were obtained over a baseline of ~ 1900 days. In contrast to the discussion in Mayor et al. (2009a), we could not confirm a longperiod trend using our new RVs and we did not detect evidence for a trend in the new CCF RVs obtained using the HARPSDRS. As shown in AngladaEscudé & Butler (2012, Fig. 3), changes in the continuum flux accross each echelle order (also called blaze function) induce RV shifts of several m s^{1} if not properly accounted for. This effect was reported to affect HARPS measurements in Pepe (2010) and appears to have been fixed by HARPSDRS v3.5 based on the release notes^{1} (issued on 29 October, 2010). With respect to HD 40307 in particular, we found that when RVs were derived without blaze function correction, a strong positive drift (~2 ma^{1} yr^{1}) was left in the Doppler time series. We speculate that the trend reported earlier may be caused by this blaze function variability. A notable feature of the HD 40307 data set is that the epochs have a long gap of 638 days between 3055–3693 [JD2 450 000]. This feature can effectively decrease the phasecoverage of the data on longer periods and complicates the interpretation of periodograms due to severe aliases.
4.1. Analysis of binned data
In a first quicklook analysis, we worked with the nightly averages of the radial velocities as obtained from HARPSTERRA using the standard setup for K dwarfs. In this setup, all HARPS echelle apertures are used and a cubic polynomial is fitted to correct for the blaze function of each aperture. After this correction, the weighted means of all RV measurements within each night e are calculated. As a result, we obtain internal uncertainties of the order of 0.3–0.4 m s^{1} for the HARPSTERRA RVs. Because of stellar and/or instrumental systematic errors, we observed that these individual uncertainties are not representative of the real scatter within most nights with five or more measurements. Using three of those nights we estimate that at least 0.6 m s^{1} must be added in quadrature to each individual uncertainty estimate. After this, the uncertainty of a given epoch is obtained as , where the sum is calculated over all exposures obtained during a given night. Finally, based on their longterm monitoring of inactive stars, Pepe et al. (2011) inferred a noise level of 0.7 m s^{1} to account for instrumental and stellar noise. After some tests, we found that adding 0.5 m s^{1} in quadrature to the uncertainties of the nightly averages ensured that none of the epochs had uncertainties below the 0.7 m s^{1} level. The typical uncertainties of a single night derived this way were of the order of 0.8 m s^{1}. These corrections are basically only welleducated guesses based on the prior experience with RV data and reported stability of the instrument. Therefore, one must be especially careful not to overinterpret the results derived from them (e.g., powers in periodograms and significance of the signals). In the fully Bayesian approach, we treat the excess noise as a free parameter of the model, therefore the Bayesian estimates of the noise properties should in principle also be more reliable.
First, we reanalysed the nighly binned RVs to see whether we could independently reproduce the results of Mayor et al. (2009a) when HARPSTERRA measurements and our Bayesian methods were used. Assuming an unknown Gaussian noise parameter (e.g. Tuomi et al. 2011) in addition to the estimated measurement uncertainties, the posterior samplings and the corresponding model probabilities easily revealed the three strong signals corresponding to periods of 4.3, 9.6, and 20.4 days. The residual periodogram of the threeKeplerian model revealed additional strong periodicities exceeding the 1% FAP level (Fig. 1, top panel) and we tested more complicated models with up to six Keplerian signals. Especially, we tested whether the additional power present in the threeKeplerian model residuals (Fig. 1, top panel) at periods of 28.6, 34.8, 51.3, and 308 days, peaking above the 10% FAP level, are statistically significant by starting our MCMC samplings at nearby seed periods.
Fig. 1 Leastsquares periodograms of the binned HD 40307 radial velocities for the residuals of the models with three (top) to six (bottom) periodic signals. The analytic 10%, 1%, and 0.1% FAPs are shown as horizontal lines. 

Open with DEXTER 
The global fourKeplerian solution was found to correspond to the three previously known superEarth signals and an additional signal with an MAP period of 320 days. This period was bounded from above and below and its amplitude was strictly positive – in accordance with our detection criteria. The corresponding posterior probability of the fourKeplerian model was 1.9 × 10^{5} times greater than that of a threeKeplerian one, making the 320day signal significant. In addition to this signal, we could identify a 51day periodicity (Fig. 1, second panel) that satisfied the detection criteria as well. Including this fifth signal in the model further increased the model probabilitity by a factor of 6.6 × 10^{6}. We could furthermore identify a sixth signal with our sixKeplerian model, correponding to a period of 34.4 days. However, even though the samplings converged well and the solution looked wellconstrained, the sixKeplerian model was only five times more probable than the fiveKeplerian one and would not be detected using our criteria in Sect. 3.2. Both of the two new significant signals had MAP estimates of their radial velocity amplitudes slightly lower than 1.0 m s^{1} – the signals at 51 and 320 days had amplitudes of 0.70 [0.31, 1.09] m s^{1} and 0.75 [0.38, 1.12] m s^{1}, respectively, where the uncertainties are denoted using the intervals corresponding to the 99% BCSs. We note that the periodogram of sampling does not have strong powers at the periods we detect (see Fig. 1 in Mayor et al. 2009a).
Given the uncertain nature of the signal at 34.4 days and the potential loss of information when using the nightly averaged RVs (artificial reduction of the number of measurements), we performed a complete Bayesian reanalysis of the full dataset (345 RVs), now including the aforementioned moving average approach to model the velocities.
4.2. Analysis using all RV measurements
The analyses of the unbinned data immediately showed the three previously announced signals (Mayor et al. 2009a) with periods of 4.3, 9.6, and 20.4 days. Modelling the data with the superposition of k Keplerian signals and an MA(3) noise model plus the two Gaussian white noise components, our posterior samplings and periodogram analyses identified these signals very rapidly, enabling us to draw statistically representative samples from the corresponding parameter densities.
The residual periodogram of this model (three Keplerians and MA(3) components of the noise removed) revealed some significant powers exceeding the 0.1% and 1% FAP level at 320 and 50.8 days, respectively (Fig. 2, top panel). Samplings of the parameter space of a fourKeplerian model indicated that the global solution contained the 320day periodicity as the fourth signal and yielded a posterior probability for the fourKeplerian model roughly 1.5 × 10^{6} times higher than for the threeKeplerian one. The nature of this signal and its relation to the stellar activity (Sect. 5) is discussed in Sect. 5.3.
Fig. 2 Leastsquares periodograms of all 345 RVs of HD 40307 for the residuals of the models with three (top) to six (bottom) periodic signals together with the analytic 10%, 1%, and 0.1% FAPs. 

Open with DEXTER 
We continued by calculating the periodogram of the residuals of our fourKeplerian model (Fig. 2, second panel) and observed a periodogram power that almost reached the 1% FAP level at a period of 50.8 days. Including this fifth signal further increased the posterior probability of our model by a factor of 6.4 × 10^{5}.
The residuals of the fiveKeplerian model contained a periodicity at 34.7 days (Fig. 2, third panel) exceeding the 1% FAP level. The corresponding sixKeplerian model with this candidate received the highest posterior probability – roughly 5.0 × 10^{8} times higher than that of the fiveKeplerian model. Since the parameters of this sixth candidate were also wellconstrained, we conclude that including the 34.7day signal in the statistical model is fully justified by the data.
We also attempted to sample the parameter space of a sevenKeplerian model but failed to find a clear probability maximum for a seventh signal (see also the residual periodogram of the sixKeplerian model in Fig. 2, bottom panel). Although the periodicity of the seventh signal did not converge to a wellconstrained probability maximum, all periodicities in the sixKeplerian model at 4.3, 9.6, 20.4, 34.7, 50.8, and 320 days were still wellconstrained, i.e. their radial velocity amplitudes were statistically distinguishable from zero and their periods had clear probability maxima. Because we were unable to constrain the parameters of a seventh signal using the posterior samplings, we cannot be sure whether the corresponding Markov chains had converged to the posterior density and cannot reliably calculate an estimate for the posterior probability of the sevenKeplerian model. Therefore we stopped looking for additional signals.
From this analysis, we can state confidently that there are six significant periodicities in the HARPSTERRA radial velocities of HD 40307 when the whole spectral range of HARPS is used. As we show in the next section, one of them has the same period as the chromospheric activity indicator (Sindex) and requires more detailed investigation. The analysis of all 345 RVs indicates that for these data binning appears to be a retrograde step in extracting periodic signals from the RV data. We infer that binning serves to alter measurement uncertainties and damp the significance levels of the periodicities in the data.
5. Stellar activity
We examined the time series of two activity indicators derived from the cross correlation function properties as provided by the HARPSDRS. They are the bisector span (BIS) and the full width at half maximum (FWHM) of the CCF. These indices monitor different features of the average stellar line. Briefly, BIS is a measure of the stellar line asymmetry and should correlate with the RVs 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 mean spectral line width. Its variability (when not instrumental) is usually associated with changes in the convective patterns on the surface of the star. A third index, the socalled Sindex in the Mount Wilson system (Baliunas et al. 1995), is automatically measured by HARPSTERRA on the blazecorrected onedimensional spectra provided by the HARPSDRS. The Sindex is a measure of the flux of the CaII H and K lines (λ_{H} = 3933.664 Å and λ_{K} = 3968.470 Å, respectively) relative to a locally defined continuum (Lovis et al. 2011) and is an indirect measurement of the total chromospheric activity of the star. For simplicity, the analysis of the activity indicators was performed throughout for the 135 nightly averaged values using sequential leastsquares fitting of periodic signals that are each described by a sinewave model (period, amplitude and phase).
5.1. Analysis of the FWHM and BIS
The BIS was remarkably stable (rms ~ 0.5 m s^{1}) and the periodogram of its time series did not show any significant powers. Visual inspection of the time series for the FWHM already shows a very significant trend of 5.3 m s^{1} yr^{1}. The 345 measurements of the FWHM are listed in Table A.1. A sinusoidal fit to this trend suggested a period of 5000 days or more (see top panel in Fig. 3). After removing the trend, two more signals strongly show up in the residuals. The first one was found at 23 days and had an analytic FAP of 0.005%. After fitting a sinusoid to this signal and calculating the residuals, an extremely significant peak appeared at 1170 days with an analytic FAP of 0.002%. After including this in a model with three sinusoids, no additional signals could be seen in the periodogram of the residuals with analytic FAP estimates lower than 10% (Fig. 3, bottom panel). We also show the FWHM values together with the fitted periodic curves in Fig. 4. While the signals in the FWHM were significant, we did not clearly detect their counterparts in the RVs. Given that BIS does not show any obvious signals either, we suspect that the periodicities in the FWHM might be caused by instrumental effects, e.g. tiny changes of the focus inside the spectrograph, rather than intrinsic variability of the stellar lines. The instrumental origin would reconcile the absence of correspondent drifts in the BIS and in the RVs. A similar indication of drifts and sensitivity of the FWHM to instrumental issues has been reported in e.g. Lovis et al. (2011). While this adds some caveats on the longterm stability of the HARPS instrumental profile (and therefore its longterm precision), unless it is found to have similar periods, we see no reason to suspect that any of the signals in the RVs are spuriously induced by changes in the FWHM (intrinsic or instrumental).
Fig. 3 Periodogram series of the signals detected in the FWHM, from most significant to less significant (top to bottom). 

Open with DEXTER 
Fig. 4 Time series of the FWHM activity index. The solid black line represents the best fit to a model containing three sinusoids (periods of 5000, 23 and 1170 days). 

Open with DEXTER 
5.2. Analysis of the Sindex
The Sindex also shows strong coherent variablity. The full set of 345 measurements of the Sindex are provided in Table A.1. Again, this analysis was performed for the nightly binned measurements using leastsquares periodograms and the sequential inclusion of sinusoidal signals. The last two Sindex measurements were well above the average and could not be reproduced by any smooth function. Even if they are representative of a physical process, such outlying points cannot be easily modelled by a series of a few sinusoids because these would add many ambiguities to the interpretation of the results. When these two points were removed, the periodograms looked much cleaner and we could unambiguously identify clear periodicities. Therefore, the analysis discussed here is based on the first 133 nightly averages only. As illustrated in the periodograms of Fig. 5, this analysis showed strong evidence for three signals in the Sindex. From most significant to less significant, we found these to consist of a longperiod quadratic trend (with a period exceeding the data baseline), a signal with P ~ 320 days and a third signal at 43 days. Figure 6 shows the best fit to a parabolic trend together with the two periodicities at 320 and 43 days. As has been reported for several other stars (Dumusque et al. 2011a), we find it likely that the longperiod trend is the signature of the stellar cycle that is not yet fully covered by the baseline of the observations. This trend might also be related to the trend observed in the FWHM, perhaps through an increase of the average magnetic field.
Fig. 5 Periodogram series of the signals detected in the Sindex, from most significant to less significant (top to bottom). 

Open with DEXTER 
Fig. 6 Timeseries of the Sindex. The solid black represents the best fit to a model containing three sinusoids (periods of 4000, 320, and 43 days). Longperiod trend and two signals at 320 and 43 days are clearly detected in the timeseries of the Sindex. 

Open with DEXTER 
Using more epochs, Lovis et al. (2011) reached the conclusion that the stellar magnetic cycle of HD 40307 is about 5000 days – a result that is compatible with ours. However, the nature of the 320day signal is less clear. This period seems to be too long to be caused by rotation because it is more than a factor of three longer than the longest currently known rotation periods which are about 100 days (Irwin et al. 2011). Long periods like this would be difficult to explain by angular momentum evolution (Barnes et al. 2012). Indeed, the third signal in the Sindex (43 days) appears to be a much better match to the expected stellar rotation of HD 40307 (Mayor et al. 2009b). The Sindex also shows a few other tentative periodicities at 26, 63, and 167 days (FAP ~ 2–5%).
HD 40307 is a very quiet star (log ; Mayor et al. 2009a) and it is in the effective temperature range where the correlation between magnetic activity and RV jitter shows the smallest correlation (~4300–5000 K, see Fig. 2 in Dumusque et al. 2011b). Therefore, apart from the two most obvious signals in the Sindices (longperiod trend and 320day period), the other tentative periodicities are unlikely to produce detectable effects on the RVs. Indeed, no RV counterpart was identified at 43 days (or any of the potential fourth signals at periods of 26, 63, or 167 days), not even at low significance. Several studies have demonstrated a clear correlation of the Sindex with observed RV offsets. For example, Dumusque et al. (2011b) performed a comprehensive analysis of simultaneous HARPS RVs and chromospheric emission measurements of G and K dwarfs, showing that the correlation between the chromospheric emission and RV offsets can be clearly traced and quantified as a function of spectral type. Pepe et al. (2011) used this correlation to detrend the RV measurements of HD 85512 (K6 V) and reported the detection of a msini ~ 2 M_{⊕} candidate with an orbital period of ~58 days. More recently, AngladaEscudé & Butler (2012) showed that longperiod signals correlated with Doppler signals can be strongly wavelength dependent even within the limited HARPS wavelength range. Since Keplerian signals cannot be wavelength depedent, this provides a new set of diagnostics to distinguish activityinduced signals from Keplerian signals. Instead of using the detrending technique employed by Pepe et al. (2011), we investigate in the next section if any of the signals are wavelength dependent.
5.3. Wavelengthdependent signals
AngladaEscudé & Butler (2012) showed that there is a correlation between the Sindices and RV measurements of HD 69830 (G8 V) and that an apparent longperiod trend in the RVs completely disappeared when one obtains the velocities using the redmost half of the HARPS spectra only. Following AngladaEscudé & Butler (2012), we call this effect chromatic jitter. Chromatic jitter can cause coherent RV variability (e.g., periodic blueshifts following the magnetic cycle) but can also behave randomly, effectively adding uncertainty and correlated noise to the measurements. At present, we do not sufficiently understand the underlying physics to be able to predict the existence and nature of the wavelength dependence of activityrelated signals. Despite this lack of detailed knowledge, the diminished sensitivity of redder wavelengths to activity has been shown to be an efficient method in verifying or falsifying the existence of proposed candidates. For example, CRIRES/VLT RV measurements (Huélamo et al. 2008; Figueira et al. 2010) in the Hband (~1.5 μm) were used to rule out the existence of a gas giant candidate around TW Hya that was initially detected in the optical RVs by Setiawan et al. (2008). Another example was the brown dwarf candidate around LP 94420. This candidate was first detected at optical wavelengths (2 km s^{1} semiamplitude) and then ruled out by the same group using nIR NIRSPEC/Keck observations (Martín et al. 2006). The chromatic nature of activityinduced signals has also abundant qualitative support from theory and simulations (e.g. Reiners et al. 2010; Barnes et al. 2011). Below we explore the wavelength dependence of the 320day signal (and all the others) to investigate its (their) nature in more detail.
As suggested by AngladaEscudé & Butler (2012), the first qualitative way of assessing the presence of chromatic jitter consists of obtaining radial velocity measurements for all epochs using the stellar spectrum from 680 nm (the longest wavelength available at HARPS) to some blue cutoff wavelength λ_{BC}, and plot the rms of the resulting RVs as a function of λ_{BC}. For a perfectly stable star, the rms should decrease monotonically as more echelle apertures are used to obtain the velocity data. However, AngladaEscudé & Butler (2012) showed that even for quiet G and K dwarfs without reported planets, a minimum in this rms was typically reached when λ_{BC} was chosen to be between 450–550 nm. For M dwarfs, this cutoff was found to be closer to 600 nm. In Fig. 7, we show the rms of the radial velocities of HD 40307 as a function of λ_{BC} and note that there appears to be a shallow minimum at λ_{BC} = 453 nm. Although this minimum is not very deep (no Keplerian solution was subtracted from the RVs), its existence already suggests the presence of wavelengthdependent noise that acts stronger in the blue part of the spectrum.
Fig. 7 Rms as a function of the bluest echelle aperture used. The corresponding λ_{BC} are illustrated at the top of the panel. A shallow minimum is reached at aperture 26 (λ_{BC} = 433 nm), indicating activity induced signals/jitter at the bluest wavelengths. 

Open with DEXTER 
The next test consists of assessing if a signal observed using the full spectrum is also present when examining a subsection of it. For simplicity, we restrict this test to the nightly binned measurements and only use periodogram tools. In agreement with 4.2, when the full spectrum was used, the fourth dominant signal was found at a period of ~320 days. Any activityinduced signals only present at the bluest wavelengths should lose their significance as one shifts λ_{BC} towards the red wavelengths. As illustrated in Fig. 8, this behaviour is clearly found for the 320day signal. Interestingly, when λ_{BC} is set to 430 nm (20th HARPS aperture), the 320day peak becomes as high as a new one emerging around 120 days. When the λ_{BC} is set to 503 nm (40th HARPS aperture), the 120day signal clearly dominates the periodogram. It is noteworthy that the powers of both the 34 and 51day peaks remain mostly constant, which shows that their significance is not seriously affected by the choice of the cutoff wavelength.
Fig. 8 Residual periodogram of the threeplanet solution as a function of the blue cutoff wavelength λ_{BC}. The top periodogram is obtained using the fullspectrum RVs. The bottom periodogram is obtained using the aperture 40 (505.0 nm) as λ_{BC}. All signals discussed in the text are marked with narrow vertical lines. The horizontal line represents the peak value of peak closer to 320 days. When λ_{BC} = 433 nm (central periodogram), the emerging 125day peak is already as high as the activityinduced signal. The 120day peak is an alias of the most likely period of 200 days for planet candidate d. 

Open with DEXTER 
As we show in the next section, the 120day signal is a yearly alias of the most favoured period at 200 days, which also appears conspicuously in the bottom periodogram of Fig. 8. To further illustrate the wavelength dependence of the 320day signal, we show the fitted RV amplitudes of the four candidate periodicities (periods of 34, 51, 200, and 320 days) as a function of λ_{BC} (Fig. 9). These amplitudes are obtained by simultaneous leastsquares fitting of seven circular orbits to the seven periods of interest (periods are also allowed to adjust). While the planet candidates at 34 and 51 days show roughly constant amplitudes, the 320day signal almost completely vanishes when λ_{BC} is redder than 450 nm (K < 0.4 m s^{1}). It is also worth noting that the 200day signal increases its amplitude as the 320day signal disappears. The 200day signal is seriously affected by yearly aliases which, effectively, correlates its amplitude with the signal at 320 days. The 320day period is so suppressed at red wavelengths that for any λ_{BC} redder than 450 nm, we were forced to keep the period fixed at 320 days to estimate its amplitude and to prevent it from converging to a very different periodicity.
Fig. 9 Amplitudes of the four signals under discussion as a function of the blue cutoff λ_{BC}. While all proposed candidates keep a constant (or slightly increase) the derived semiamplitude, the 320day signal is almost nonexistent (the period actually becomes unconstrained) when wavelengths redder than 480 nm are used to derived a sevenplanet orbital solution. 

Open with DEXTER 
At this point, it is fair to ask why the longterm variability in the Sindex (and in the FWHM) does not show in the RV data as well. We can show that it actually does, very prominently, when calculating the RVs using the bluemost part of the spectrum only. To illustrate this, we obtained the RV measurements using the wavelength range from 380 nm to 440 nm (20 bluest echelle apertures). After subtracting the signals of the three planets reported by Mayor et al. (2009a), a parabolic shape is clearly left in the RV residuals. Figure 10 illustrates how well the parabolic trend matches the longperiod trend also observed in the Sindex (compare Figs. 6 and 10). No blue RV counterparts of the other tentative periodicities in the Sindex were identified in the data.
Fig. 10 Residuals of the threeplanet fit when only the bluemost part of the spectra were used to extract the Doppler measurements. A quadratic trend similar to the one detected in the Sindex is the most obvious signal left in these residuals. 

Open with DEXTER 
As a conclusion to this section, we identified the chromatic nature of the ~320day signal. We note that the internal errors of the individual spectra, given in the third column of Table A.1, have a median value of 0.4–0.5 m s^{1}. This number is significantly smaller than the measured intranight scatter and also smaller than the best rms we could receive from any six or sevenplanet fit. This means that systematic noise (stellar and/or instrumental) is a dominant source of uncertainty and that accuracy is not lost by ignoring the Doppler information at the bluest wavelengths. A careful optimisation of λ_{BC} might provide even more significant detections but, given the unknown nature of this chromatic jitter, choosing λ_{BC} = 490.0 nm (37th HARPS aperture) appears to provide a fair compromise. We suspect that the observed behaviour of HD 40307 RVs (also HD 69830, HD 85512, and several other nominally stable stars, see AngladaEscudé & Butler 2012) is a manifestation of the same processes observed in magnetically active stars but occuring at lower RV amplitude levels. The correlation between the Sindex and the occurrence of the RV period at 320 days allows us to establish a very likely connection between this activityrelated signal and a spurious RV signal and rule out the Keplerian nature of the corresponding candidate. We are in the process of analysing the chromatic dependence of stellar jitter using a more representative sample (i.e. the sample analysed in Dumusque et al. 2011a), the results of which will be presented in a future publication. Given the substantial changes in the potential signals in the system (e.g., appearance of the candidate at a period of 125 or 200 days), we perform a complete orbital reanalysis of the redmost half of the spectrum only (490–680 nm) in the next section.
Overview of the significant periodicities found in two activity indicators.
6. Analysis of the radial velocities from the 490–680 nm spectra
First, we performed a rapid analysis of the RVs obtained using the redmost part of the spectra using leastsquares periodograms. As was found to be the case when analysing the RVs from full spectra, three additional periodicities compared to the (Mayor et al. 2009a) solutions could be extracted from the data if all orbital eccentricities were fixed to zero. Two of the signals had the same periods (34 and 51 days) and the third one showed at ~200 days (and on its yearly alias at 125 days). Although a model with circular orbits always looks more appealing aesthetically, we admit that this assumption cannot be justified using a pure leastsquares approach without the use of prior information and/or physical constrains (e.g., orbital dynamics). As before, and in order to obtain a robust detection of all the candidates, we performed a full Bayesian analysis of the 345 red RV measurements using the methods and models described in Sect. 3.2.
The parameter spaces of models with k = 0,...,6 were easy to sample using the adaptive Metropolis algorithm and the Markov chains converged rapidly to the posterior densities. Again, the three candidates presented by Mayor et al. (2009a) with periods of 4.3, 9.6, and 20.4 days were extremely significant and easy to find using our MCMC samplings.
After removing the MA(3) component of the noise and the signals in the threeKeplerian model, the power spectrum of the residuals showed very strong peaks at periods of 50.8, 82.1, and 128.3 days (Fig. 11, top panel) exceeding 0.1% FAP levels. Even though the highest peak corresponds to the 125day period, the MCMC samplings showed that the 50.8day periodicity is the global solution by providing higher posterior probabilities to the fourKeplerian model. We note that when calculating the residual periodogram of the threeKeplerian model, it is necessary to assume a fixed 3Keplerian solution (and noise model) for the data. In contrast, the MCMC automatically accounts for correlations by also adjusting the noise parameters simultaneously to find the most probable areas of the parameter space (not only minimising the leastsquares sum). After several MCMC samplings starting in the vicinity of the different possible periods suggested by the periodogram, this 50.8day signal was always found to correspond to the global fourKeplerian solution and made the fourKeplerian model 3.1 × 10^{12} times more probable than the threeKeplerian one.
Fig. 11 As in Fig. 2 but for the unbinned radial velocities of the 490–680 nm spectrum. 

Open with DEXTER 
The residual periodograms of four and fiveKeplerian models had moderate powers at 34.6, 81, 125, and 197 days (Fig. 11, panels two and three). We performed similar samplings of the parameter spaces of five and sixKeplerian models using the favoured periodogram powers as starting points of the Markov chains. As a result, we found that the period space near the 197 and 34.6day powers provided the global maxima and improved the posterior model probabilities by factors of 3.6 × 10^{11} and 2.5 × 10^{9} when moving from k = 4 to k = 5 and from k = 5 to k = 6, respectively, where k is the number of Keplerian signals in the model. Removing all six signals and the MA(3) components of the noise from the data left the periodogram of the residuals without any considerable powers (Fig. 11, bottom panel). Similarly as for the full spectrum, the period of the seventh signal in a sevenKeplerian model did not converge to a clear probability maximum, thus failing to satisfy our requirement that the period of a signal has to be wellconstrained for a positive detection.
Relative posterior probabilities of models with k = 0,...,6 Keplerian signals (ℳ_{k}) and the periods (P_{s}) of the signals added to the model when increasing the number of signals in the model by one.
To summarise, the posterior probabilities (and the corresponding Bayesian evidences) of models with k = 0,1,2,... Keplerian signals were found to be heavily in favour of k = 6 (Table 3). This was the case although the selected prior probabilities penalise the models more heavily as the number of signals increases. We do not show the results for a sevenKeplerian model in Table 3 because the posterior samplings did not converge to a clear maximum (nor several maxima) in the sevenKeplerian model and therefore we cannot be sure whether the OBMH estimate yielded trustworthy values for the marginal integrals needed to assess the model probabilities. We plot the phasefolded MAP estimated orbits of the six Keplerian signals in Fig. 12 after removing the MA(3) components.
Fig. 12 Phasefolded MAP Keplerian signals of the sixplanet solution with the other five signals and the MA(3) components removed from each panel. 

Open with DEXTER 
Sixplanet solution of HD 40307 radial velocities.
Fig. 13 Distributions estimating the posterior densities of orbital periods (P_{x}), eccentricities (e_{x}), and radial velocity amplitudes (K_{x}) of the six Keplerian signals. The dashed lines show the prior densities for comparison for the orbital eccentricities. The solid curves are Gaussian densities with the same mean (μ) and variance (σ^{2}) as the parameter distributions. Additional statistics, mode, skewness (μ^{3}), and kurtosis (μ^{4}) of the distributions are also shown. 

Open with DEXTER 
One of our detection criteria for Keplerian signals in the radial velocities is that their amplitudes should be statistically significantly different from zero. The MAP estimates and the corresponding 99% BCSs of the sixKeplerian model parameters are shown in Table 4. The BCSs of the RV amplitudes clearly indicate that all six amplitude parameters indeed differed significantly from zero in accordance with the model probabilities that imply the existence of six Keplerian signals in the data. This was clearly the case despite the 34.7 and 198day signals having amplitudes with MAP estimates below 1.0 m s^{1}. Also, the periods of all signals were wellconstrained and had clear MAP values and closeGaussian densities (Table 4; Fig. 13). The distributions of orbital eccentricities are also shown in Fig. 13. In agreement with the quicklook detection using leastsquares periodograms assuming circular orbits, these distributions demonstrate that all signals had eccentricities consistent with closecircular orbits.
We also analysed the RVs derived from two independent blocks of the red part of the spectrum to assess if any of the signals was only present at certain wavelength ranges. The first block comprised the range between 465 and 553 nm (from now on Red1) and the second one covered the range between 553 to 690 nm (Red2). Red1 contains the richest part of a HARPS spectrum of a midK dwarf in terms of Doppler information (higher signaltonoise, but also higher density of lines). The mean uncertainty of Red1 is approximately 0.6 m s^{1}. As a result, we found that the same six signals were directly spotted even with periodograms only. The amount of Doppler information in Red2 is lower and nominal uncertainties were found to be approximately 1.0 m s^{1} on average. Using sequential periodogram subtraction of the signals we could recover the same six candidates except that the 125day alias was preferred against the 200day signal for HD 40307 g. However, the Bayesian samplings including the complete model for the uncertainties still converged to the same sixplanet solution (P_{g} ~ 200 days). Given that there is no a priori reason to prefer one or the other (125 versus 200 days, no activity signals appear to be related to any of them) and because the MCMC samplings insist on favouring the 200day period, we conclude that the solution in Table 4 is the one preferred by the data. The significances and model probabilities using Red1 and Red2 are obviously lower than when using the whole red half of the spectrum. Yet, with the exception of the 34day candidate in Red2 (chains converged to the same solution, but model probability ratio was close to the 1:150 threshold), all signals could be confidently identified even using two subsets of the HARPS data (each contains 1/3 of the full HARPS wavelength coverage). These tests give us further confidence in the reality of the reported candidate signals.
As a final consistency check, we calculated the model probabilities using different information criteria and compared them to the Bayesian model probabilities derived using our approximation for the marginal integrals. First, we estimated the relative posterior probabilities using the Akaike information criterion (AIC; see e.g. Liddle 2007; Burnham & Anderson 2002), which is known to provide only rough approximations for the marginal integrals. The AIC (actually its smallsample version, AICc) yielded posterior probabilities for the different models that did not differ by more than a factor of 5 from the probabilities received using the OBMH estimate (Table 3). Using the AIC in model selection yielded the same qualitative results the OBMH estimate did. A common criticism of the AIC is that does not penalise the more complex models as much as it should in all applications (e.g. Burnham & Anderson 2002, and references therein). To obtain a more reliable comparison, we repeated the model selection procedure using the Bayesian information criterion (BIC), sometimes called the Schwarz information criterion (Schwarz 1978), which is known to penalise more complicated models, e.g. models with more Keplerian signals in this particular case, more heavily than the AIC. Even according to the more restictive BIC, the sixKeplerian model had the highest posterior probability and was 1.3 × 10^{5} times more probable than the fiveKeplerian model. These results provide further confidence in the conclusion that six periodic signals explain the variations in the data.
6.1. Noise properties
In addition to ruling out wavelengthdependent signals, the Bayesian analysis of the RVs at different wavelength ranges allows the comparison of the noise properties. As described in Sect. 3.2, we modelled the radial velocity noise using a thirdorder MA model (three free parameters) together with a Gaussian white noise component (one free parameter). The posterior distributions and MAP values for these parameters quantify the significance and the magnitude of each component of the noise (white vs correlated).
The MA parameters φ_{j}, j = 1,2,3, which describe the correlation between the error of the ith and (i − j)th measurement, are clearly positive for the RVs from full spectra (Table 5). The noise of the measurements appears to correlate positively the strongest with the noise of the previous measurement and the effect decreases for the measurements obtained before that. However, this correlation decreases for the RVs obtained from the redmost part of the spectra (490–680 nm). According to our results, this decrease is balanced by a corresponding increase in the magnitude of the Gaussian white noise component (Table 5). Therefore, the noise in the radial velocities becomes “whiter” when excluding the bluest half of the spectra from the analyses. This result means that in addition to removing longterm activityrelated signals, using redder wavelengths also helps reducing correlations in the noise that might cause biases to any estimates of the Keplerian signals if not accounted for.
Noise parameters obtained when using the RVs from the full spectra and the red part (490–680 nm) of the spectra only as denoted using the MAP estimates and the 99% BCSs.
We note that the six signals in the velocities were rather independent of the exact MA model used. We obtained consistent solutions with MA(2) and MA(4) models as well, though these noise models were found to have slightly lower posterior probabilities. This indicates, as can be seen in Table 5, that the MA components MA(p) with p > 1 do not affect the solution much because their 99% BCSs imply that they are consistent with zero.
7. Orbital stability
We investigated the longterm stability of the planetary system corresponding to our sixKeplerian solution (Table 4). This analysis was performed in three steps. First, we checked the Lagrange stability criteria, which provide a first qualitative assessement of the allowed orbital configurations. Second, we numerically integrated orbits directly drawn from the MCMC samplings and assessed the fraction corresponding to stable orbits at the Myr time scale. As a last test, we analysed the stability of the sixplanet configuration by testing different eccentricities and arguments of the node of each planet.
7.1. Lagrange stability criteria
As a rough initial test, we checked whether the subsequent planets in the system satisfied the approximate Lagrange stability criterion expressed analytically in Barnes & Greenberg (2006) and applied for instance in Tuomi (2012). We plotted the resulting areas of instability in Fig. 14 – if the i − 1th or i + 1th planet has orbital parameters within the shaded areas surrounding the ith planet, the system is likely unstable over the long term. While this criterion is but a rough approximation of reality because it does not take the stabilising and/or destabilising effects of orbital resonances into account and is only valid for two planets at a time, it provides simple guidelines on the prospects of stability in the system.
According to the Lagrange stability criterion and the corresponding analytical thresholds in Fig. 14, the six planet candidates orbiting HD 40307 form a stable system if their orbital eccentricities are close to or below their MAP estimates. Alternatively, it can be stated that if the planets exist and have masses approximately equal to the minimum masses in Table 4, they are all likely on closecircular orbits. Especially, the outer companion is unlikely to have orbtital eccentricity in excess of 0.4 (Fig. 14).
Fig. 14 Approximate Lagrange stability thresholds (shaded areas) between each two planets (red circles) and the MAP orbital parameters of the sixplanet solution (Table 4). 

Open with DEXTER 
Fig. 15 Distributions estimating the orbital eccentricities given measurements alone (grey histogram; mean μ_{m} and standard deviation σ_{m}), showing the values the eccentricities had during orbital integrations of 10^{6} years that could not be shown to be unstable (white histogram; mean μ_{d} and standard deviation σ_{d}), and their combination (red histogram). 

Open with DEXTER 
Fig. 16 Stability analysis of the circular sixplanet solution (S1). The panels show the influence of changing the initial eccentricity (left) and initial mean anomaly (right) on the stability of the complete planetary system. “Red” corresponds to highly unstable orbits while “black” denotes regular motion. Each of the 100 initial conditions per planet were integrated numerically for 1 kyr, for which we then computed the stability criterion D (Eq. (2)). 

Open with DEXTER 
7.2. Numerical integration of MCMC samples
To further test the dynamical stability of the system, we integrated the planetary orbits for 10^{6} years, which is appropriate to exclude “almost all” unstable configurations (Barnes & Quinn 2004). In our integrations, we used the BulirschStoer integrator (Bulirsch & Stoer 1966) that has been used before when integrating planetary orbits (e.g. Tuomi et al. 2009; Tuomi 2011). We found that most of the unstable configurations corresponded to one (or more) of the eccentricities of the companions being higher than its MAP estimate and the system turned out to be unstable in 10^{4} − 10^{5} years. Specifically, we drew a sample from the posterior density of the orbital parameters and used each vector in this sample as an initial state of orbital integrations. We increased the sample size such that there were 100 initial states that could not be shown to correspond to unstable orbital configurations. This subsample was found to correspond to 4% of the total sample.
In Fig. 15 we plot the posterior densities of the planetary eccentricities together with the densities obtained from the subsample of solutions that were not found to be unstable. Figure 15 states essentially the same thing as Fig. 14. That is, it shows that the stable configurations mostly correspond to orbital eccentricities closer to zero than those derived from the Bayesian MCMC analysis. The corresponding updated estimates for orbital eccentricities are shown in Table 4 and denoted using e_{dyn}. While this experiment shows that a nonnegligible fraction of the solutions allowed by the data corresponds to stable configurations, we aimed to explore the source of the dynamical instabilities and the role of potential dynamical resonances in the system as well.
7.3. Sources of (in)stability
Owing to the dense packing of the planets in the HD 40307 system, one still has to rely on systematic numerical integrations to correctly investigate any aspects of its stability. For example, from a direct numerical integration of the MAP sixplanet solution given in Table 4, one quickly finds that the orbital evolution leads to collisional trajectories within a timescale of 10^{3} years.
Fig. 17 Representation of the liquid water habitable zone around HD 40307 as described by Selsis et al. (2007) compared to our solar system. Note the compactness of the orbital configuration of the five inner planet candidates around HD 40307 compared to the planets in the inner solar system. 

Open with DEXTER 
As was already shown above, it is likely that the planets have closecircular orbits. This configuration is also supported by the confidence intervals given for the individual eccentricities in Table 4, which essentially state that circular orbits cannot be ruled out by the data. Therefore, to slightly simplify the dynamical analysis, we performed a new fit to the observations, for which we allowed only circular orbits. For this nominal solution (S1 hereafter) we analysed the impact of changing the initial eccentricity e and initial mean anomaly M of individual planets on the stability of the whole system.
Starting from S1 we changed only the initial e or M of an individual planet. We then integrated the obtained system for 1000 years using the symplectic SABA2 integrator (Laskar & Robutel 2001) with a time step τ = 0.1 days. The stability of the system was afterwards estimated using a frequency analysis (Laskar 1993). For this we divided the solution into two parts and computed the mean motions n for all planets for each part separately. To obtain the measure of stability (D), we calculated the maximum relative change in the difference between both mean motions as (2)For a stable system, the individual mean motions will not vary with time, i.e. . The results are shown in Fig. 16. Here “red” denotes strongly chaotic orbits, while “black” corresponds to very regular ones.
The left panel of Fig. 16 confirms that the stability of the system is crucially dependent on the individual eccentricities of the planets. Only for initial eccentricities very close to zero regular motion can be found (D ~ 10^{9}). Using the right panel of Fig. 16 it is also possible to identify the reasons for dynamical stability in HD 40307. While the system is not sensitive at all for example to the initial mean anomaly of planet g, it crucially depends on M_{e} and M_{f}. In each column of these planets one finds in Fig. 16 initial conditions that lead to very regular motion (D ~ 10^{9} for e.g. M_{f}(t_{0}) ≈ 310°) and also conditions that lead to chaotic motion (D ~ 10^{2} for e.g. M_{f}(t_{0}) ≈ 267°). When analysing the results of the respective integrations, we found that the orbital periods of candidates e and f need to have a ratio of 3:2 to maintain stable motion. This resonant configuration helps to avoid close encounters between the planets due to the libration of the resonance angle θ = −2λ_{e} + 3λ_{f} − ω_{e}. We note that this configuration is also maintained over long time intervals. Integrations up to 1 Myr using these initial conditions did not show any signs of chaotic motion.
In contrast, for highly unstable initial conditions the angle θ circulates, which leads very soon to collisions between some planets or the ejection of one of them.
In conclusion to the discussion on the dynamics of the system we find that – despite the presence of many planet candidates in the system – there exists a substantial number of dynamically stable orbits allowed by the data. In particular, orbits with eccentricities close to zero are dynamically favoured due to the existence of the stabilising 3:2 mean motion resonance between candidates e and f. We note that we have only tested solutions based on the minimum masses of the candidates. It is possible that increasing the masses of all planets by a small factor (e.g., × 1.2) makes the system unstable and could therefore impose interesting upper limits to the masses. Additional RV measurements are necessary to better constrain the eccentricities of the several candidates and derive more robust conclusions on the dynamical features of the system.
8. Habitability and prospective followup
The socalled liquid water habitable zone (HZ) is the area around a star where an Earthlike planet could support liquid water on its surface. For HD 40307, this HZ lies between 0.43 and 0.85 AU (see Fig. 17). This interval was computed using the recipes given in Selsis et al. (2007) and is a function of stellar luminosity and effective temperature. We used the luminosity and temperature estimates given by Ghezzi et al. (2010). With an orbital semimajor axis of 0.60 AU, the planet candidate HD 40307 g receives ~62% of the radiation the Earth receives from the Sun and lies safely within the HZ of its stellar host. Even though the radiation is somewhat low compared to that received by the Earth, we note that the Earth lies actually reasonably close to the inner boundary of the Sun’s HZ (Fig. 17).
As is the case with the previously reported planets in the HZs of nearby stars, additional observational information is necessary – on top of confirming the planetary nature of the three new signals we report in this work – to decide if HD 40307 g indeed can support water on its surface. In the meantime, and given that it would be a rather massive object for a telluric planet, detailed climatic (e.g. Barnes et al. 2012) and planetary interior simulations (e.g. Korenaga 2010; Stein et al. 2011) can be used to estimate its suitability for supporting liquid water and, perhaps, life. Unlike the other previously reported candidate habitable planets around nearby M dwarfs (e.g., GJ 581 and GJ 667C), HD 40307 g is located farther away from the central star and, like the habitable zone planet candidate Kepler 22b (Borucki et al. 2012), it should not suffer from tidal locking either (Kasting et al. 1993; Barnes et al. 2009b).
Given its relatively longperiod orbit, the transit probability of planet candidate HD 40307 g is very low (~0.6% Charbonneau et al. 2007). Also, previous attempts to detect the transits of the inner candidate HD 40307 b (P ≈ 4.3 days) have been unsuccessful (Gillon et al. 2010). Therefore, if coplanarity is a common feature of planetary systems, the fact that HD 40307 b is not transiting its star bodes ill for the transit chances of any outer companion. HD 40307 is a nearby star (13 pc) compared to the typical target of the Kepler mission (e.g., Kepler 22 is located at ~200 pc). Because of this, the starplanet angular separation can be as large as ~46 mas. Spacebased directimaging missions such as Darwin/ESA (Cockell et al. 2009) or NASA/TPF (Lawson et al. 2006) aim to reach the sensitivity to detect Earthsized planets at an inner working angle of ~40 mas. Therefore, HD 40307 g would be the first habitablezone candidate in the superEarth mass regime that could be targeted by the planned direct imaging observatories.
9. Discussion
After examining the wavelength dependence of the signals in the RVs of HD 40307 and their relation to the chromospheric activity indicator (Sindex), we easily detect the three superEarth candidates previously reported by Mayor et al. (2009a) using independent data analysis methods. In addition to these, we also confidently detect three new periodic signals in the RVs of this K2.5 V dwarf star. Our analyses of the HARPS RVs indicates the existence of a system with six superEarth candidates around this star (Table 3; Fig. 11) and our dynamical analyses show that longterm stable orbits are allowed within our credibility sets of the parameters. Our results make HD 40307 one of the few known planetary systems with more than five planet candidates: HD 10180 (Lovis et al. 2011; Tuomi 2012), the Kepler11 sixplanet system (Lissauer et al. 2011), and the solar system.
A handful of potentially habitable planets have been reported to date. GJ 581 d (Mayor et al. 2009b) and HD 88512 b (Pepe et al. 2011) orbit at the extended HZ where thick clouds of CO_{2} or water would be required to sustain wet atmospheres. As is the case with Kepler22 b (Borucki et al. 2012) and GJ 667C c (AngladaEscudé et al. 2012), HD 40307 g orbits well within the more classic and conservative definition of the HZ. With a radius of 2.4 R_{⊕}, Kepler22 b is likely to be a smallscale version of Neptune rather than a rocky planet with a solid surface (Borucki et al. 2012). However, the precise natures of GJ 667C c and HD 40307 g are unknown and they could also be scaleddown versions of Neptunelike planets with thick atmospheres and solid cores. Compared to the candidates orbiting nearby lowmass stars (GJ 581 d, GJ 667C c, HD 85512 b), HD 40307 g is not likely to suffer from tidallocking – which improves its chances of hosting an Earthlike climate. Like for the other HZ candidates, it is not yet possible to determine its physical and geochemical properties and a directimaging mission seems to be the most promising way of obtaining observational information on these properties.
For stable stars, the HARPS data can be readily used to detect Keplerian periodicities in RVs with amplitudes below 1.0 m s^{1} (Pepe et al. 2011). This accuracy can be additionally improved by obtaining RV measurements using optimized datareduction techniques (e.g. HARPSTERRA; AngladaEscudé & Butler 2012) and analysing the wavelength dependence of the signals. Also, detections should not rely on criteria solely based on the analyses of the model residuals – this procedure is prone to biases and the weakest signals in the data can easily escape detection (compare e.g. the results of Lovis et al. 2011; Feroz et al. 2011; Tuomi 2012).
The ability to detect such lowamplitude signals also introduces a practical problem. The orbits of the planets can no longer be readily compared with large numbers of data points using phasefolded curves (Fig. 12) because – while the signals may be significant – their existence cannot be so readily verified by eyesight from such figures. Therefore, when discovering lowmass planets with the radial velocity method, it is necessary to demonstrate the significance of the signals using some other means. The detection criterion we describe, namely, that the radial velocity amplitude is statistically distinguishable from zero, provides an easy way of demonstrating this. The distributions of all amplitude parameters in Fig. 13 clearly satisfy this criterion and show that all signals we report are indeed significantly present in the data.
Additional measurements are needed to verify the existence and better constrain the orbits and minimum masses of the three new planet candidates reported here. Moreover, given that all signals are of planetary origin, detailed dynamical analyses can additionally constrain the orbits and masses by excluding unstable configurations, and especially, constraining the upper limits of the planetary masses that could then be used to determine the inclinations of their orbits. In this sense, and given the relatively dense dynamical packing of the system, the planet candidates around HD 40307 are likely to have masses close to their minimum values.
The outer planet candidate (HD 40307 g) has a minimum mass of ~7 M_{⊕} and orbits the star within the liquid water HZ. Given its priviledged position comfortably within the HZ, it is likely to remain inside it for the full orbital cycle even if its eccentricity is not strictly zero. A more detailed characterization of this candidate is very unlikely using ground based studies because it is very inlikely to transit the star, and a direct imaging mission seems the most promising way of learning more about its possible atmosphere and lifehosting capabilities. The planetary system around HD 40307 has an architecture radically different from that of the solar system (see Fig. 17), which indicates that a wide variety of formation histories might allow the emergence of roughly Earthmass objects in the habitable zones of stars.
Acknowledgments
M. Tuomi is supported by RoPACS (Rocky Planets Around Cool Stars), a Marie Curie Initial Training Network funded by the European Commission’s Seventh Framework Programme. G. AngladaEscudé is supported by the German Ministry of Education and Research under 05A11MG3. E. Gerlach would like to acknowledge the financial support from the DFG research unit FOR584. A. Reiners acknowledges research funding from DFG grant RE1664/91. S. Vogt gratefully acknowledges support from NFS grant AST0307493. This work is based on data obtained from the ESO Science Archive Facility under request number GANGLFGGCE173161. This research has made extensive use of the SIMBAD database, operated at CDS, Strasbourg, France; and the NASA’s Astrophysics Data System. We acknowledge the significant efforts of the HARPSESO team in improving the instrument and its data reduction pipelines, and obtaining the observations that made this work possible.
References
 AngladaEscudeé, G., & Butler, R. P. 2012, ApJS, 200, 15 [NASA ADS] [CrossRef] (In the text)
 AngladaEscudé, G., LópezMorales, M., & Chambers, J. E. 2010, ApJ, 709, 168 [NASA ADS] [CrossRef] (In the text)
 AngladaEscudé, G., Arriagada, P., Vogt, S. S., et al. 2012, ApJ, 751, L16 [NASA ADS] [CrossRef] (In the text)
 Baliunas, S. L., Donahue, R. A., Soon, W. H., et al. 1995, ApJ, 438, 269 [NASA ADS] [CrossRef] (In the text)
 Barnes, S. A. 2007, ApJ, 669, 1167 [NASA ADS] [CrossRef] (In the text)
 Barnes, R., & Greenberg, R. 2006, ApJ, 647, L163 [NASA ADS] [CrossRef] (In the text)
 Barnes, R., & Quinn, T. 2004, ApJ, 611, 494 [NASA ADS] [CrossRef] (In the text)
 Barnes, R., Jackson, B., Raymond, S. N., et al. 2009a, ApJ, 695, 1006 [NASA ADS] [CrossRef] (In the text)
 Barnes, R., Jackson, B., Greenberg, R., & Raymond, S. N. 2009b, ApJ, 700, L30 [NASA ADS] [CrossRef] (In the text)
 Barnes, J. R., Jeffers, S. V., & Jones, H. R. A. 2011, MNRAS, 412, 1599 [NASA ADS] [CrossRef] (In the text)
 Barnes, R., Mullins, K., Goldblatt, C., et al. 2012, Astrobiology, accepted [arXiv:1203.5104] (In the text)
 Batalha, N. M., Borucki, W. J., Bryson, S. T., et al. 2011, ApJ, 729, 27 [NASA ADS] [CrossRef] (In the text)
 Borucki, W. J., Koch, D. G., Batalha, N., et al. 2012, ApJ, 745, 120 [NASA ADS] [CrossRef] (In the text)
 Bulirsch, R., & Stoer, J. 1966, Numer. Math., 8, 1 [CrossRef] (In the text)
 Burnham, K. P., & Anderson, D. R. 2002, Model Selection and Multimodel Inference, 2nd edn. (New York: SpringerVerlag) (In the text)
 Charbonneau, D., Brown, T. M., Burrows, A., & Laughlin, G. 2007, Proc. Protostars and Planets V (University of Arizona Press), 951, 701 [NASA ADS] (In the text)
 Chib, S., & Jeliazkov, I. 2001, J. Am. Stat. Ass., 96, 270 [CrossRef] (In the text)
 Clyde, M. A., Berger, J. O., Bullard, F., et al. 2007, Statistical Challenges in Modern Astronomy IV, eds. G. J. Babu, & E. D. Feigelson, ASP Conf. Ser., 371, 224 (In the text)
 Cockell, C. S., Léger, A., Fridlund, M., et al. 2009, Astrobiology, 9, 1 [NASA ADS] [CrossRef] [PubMed] (In the text)
 Cumming, A. 2004, MNRAS, 354, 1165 [NASA ADS] [CrossRef] (In the text)
 Dumusque, X., Udry, S., Lovis, C., et al. 2010, A&A, 525, A140 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Dumusque, X., Santos, N. C., Udry, S., et al. 2011a, A&A, 527, A82 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Dumusque, X., Lovis, C., Udry, S., & Santos, N. C. 2011b, Proc. IAU Symp., 276, 530 [NASA ADS] (In the text)
 Feroz, F., Balan, S. T., & Hobson, M. P. 2011, MNRAS, 415, 3462 [NASA ADS] [CrossRef] (In the text)
 Figueira, P., Pepe, F., Melo, C. H. F., et al. 2010, A&A, 511, A55 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Ford, E. B., & Gregory, P. C. 2007, Statistical Challenges in Modern Astronomy IV, eds. G. J. Babu, & E. D. Feigelson, ASP Conf. Ser., 371, 189 (In the text)
 Funk, B., Wuchterl, G., Schwarz, R., et al. 2010, A&A, 516, A82 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Ghezzi, L., Cunha, K., Smith, V. V., et al. 2010, ApJ, 720, 1290 [NASA ADS] [CrossRef] (In the text)
 Gray, R. O., Corbally, C. J., Garrison, R. F., et al. 2006, AJ, 132, 161 [NASA ADS] [CrossRef] (In the text)
 Gillon, M., Deming, D., Demory, B.O., et al. 2010, A&A, 518, A25 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Haario, H., Saksman, E., & Tamminen, J. 2001, Bernoulli, 7, 223 [CrossRef] [MathSciNet] (In the text)
 Hastings, W. 1970, Biometrika, 57, 97 [CrossRef] [MathSciNet] (In the text)
 Huélamo, N., Figueira, P., Bonfils, X., et al. 2008, A&A, 489, L9 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Irwin, J., Berta, Z. K., Burke, C. J., et al. 2011, ApJ, 727, 56 [NASA ADS] [CrossRef] (In the text)
 Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Ass., 90, 773 [CrossRef] [MathSciNet] (In the text)
 Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108 [NASA ADS] [CrossRef] [PubMed] (In the text)
 Korenaga, J. 2010, ApJ, 752, 43 [CrossRef] (In the text)
 Laskar, J. 1993, Celest. Mech. Dyn. Astron., 56, 191 [NASA ADS] [CrossRef] (In the text)
 Laskar, J., & Robutel, P. 2001, Celest. Mech. Dyn. Astron., 80, 39 [NASA ADS] [CrossRef] (In the text)
 Lawson, P. R., Ahmed, A., Gappinger, R. O., et al. 2006, Proc. SPIE, 6268, 626828 [CrossRef] (In the text)
 Léger, A., Rouan, D., Schneider, J., et al. 2009, A&A, 506, 287 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] (In the text)
 Liddle, A. R. 2007, MNRAS, 377, L74 [NASA ADS] (In the text)
 Lissauer, J. J., Fabrycky, D. C., Ford, E. B., et al. 2011, Nature, 470, 53 [NASA ADS] [CrossRef] [PubMed] (In the text)
 Lomb, N. R. 1976, Astrophys. Space Sci., 39, 447 [NASA ADS] [CrossRef] (In the text)
 Lovis, C., Ségransan, D., Mayor, M., et al. 2011, A&A, 528, A112 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Martín, E. L., Guenther, E., Zapatero Osorio, M. R., Bouy, H., & Wainscoat, R. 2006, ApJ, 644, L75 [NASA ADS] [CrossRef] (In the text)
 Mayor, M., Pepe, F., Queloz, D., et al. 2003, ESO Messenger, 114, 20 [NASA ADS] (In the text)
 Mayor, M., Udry, S., Lovis, C., et al. 2009a, A&A, 493, 639 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Mayor, M., Bonfils, X., Forveille, T., et al. 2009b, A&A, 507, 487 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., et al. 1953, J. Chem. Phys., 21, 1087 [NASA ADS] [CrossRef] (In the text)
 Pepe, F. 2010, From HARPS/HARPSN to ESPRESSO: pushing the limits further, Astronomy of Exoplanets with Precise Radial Velocities, exoplanets.astro.psu.edu/workshop/ [exoplanets.astro.psu.edu], ed. C. Beichman, et al. (In the text)
 Pepe, F., Lovis, C., Ségransan, D., et al. 2011, A&A, 534, A58 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Reiners, A., Bean, J. L., Huber, K. F., et al. 2010, ApJ, 710, 432 [NASA ADS] [CrossRef] (In the text)
 Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef] (In the text)
 Schwarz, G. E. 1978, Ann. Stat., 6, 461 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
 Selsis, F., Kasting, J. F., Levrard, B., et al. 2007, A&A, 476, 1373 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Setiawan, J., Henning, T., Launhardt, R., et al. 2008, Nature, 451, 38 [NASA ADS] [CrossRef] [PubMed] (In the text)
 Sousa, S. G., Santos, N. C., Mayor, M., et al. 2008, A&A, 487, 373 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] (In the text)
 Stein, C., Finnenkötter, A., Lowman, J. P., & Hansen, U. 2011, GeoRL 38, Issue 21 (In the text)
 Tuomi, M. 2011, A&A, 528, L5 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Tuomi, M. 2012, A&A, 543, A52 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Tuomi, M., & Kotiranta, S. 2009, A&A, 496, L13 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Tuomi, M., Kotiranta, S., & Kaasalainen, M. 2009, A&A, 494, 769 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Tuomi, M., Pinfield, D., & Jones, H. R. A. 2011, A&A, 532, A116 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Tuomi, M., Jones, H. R. A., Jenkins, J. S., et al. 2012a, MNRAS, submitted (In the text)
 Tuomi, M., Jones, H. R. A., Jenkins, J. S., et al. 2012b, A&A, submitted (In the text)
 van Leeuwen, F., 2007, A&A, 474, 653 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Vogt, S. S., Allen, S. L., Bigelow, B. C., et al. 1994, Proc. SPIE, 2198, 362 [NASA ADS] [CrossRef] (In the text)
 Winn, J. N., Matthews, J. M., Dawson, R. I., et al. 2011, ApJ, 737, L18 [NASA ADS] [CrossRef] (In the text)
Online material
Appendix A: Radial velocities and activity indices
HARPSTERRA differential radial velocities of HD 40307 from the full spectra (F) and the red part (490–680 nm) of the spectra (R).
All Tables
Relative posterior probabilities of models with k = 0,...,6 Keplerian signals (ℳ_{k}) and the periods (P_{s}) of the signals added to the model when increasing the number of signals in the model by one.
Noise parameters obtained when using the RVs from the full spectra and the red part (490–680 nm) of the spectra only as denoted using the MAP estimates and the 99% BCSs.
HARPSTERRA differential radial velocities of HD 40307 from the full spectra (F) and the red part (490–680 nm) of the spectra (R).
All Figures
Fig. 1 Leastsquares periodograms of the binned HD 40307 radial velocities for the residuals of the models with three (top) to six (bottom) periodic signals. The analytic 10%, 1%, and 0.1% FAPs are shown as horizontal lines. 

Open with DEXTER  
In the text 
Fig. 2 Leastsquares periodograms of all 345 RVs of HD 40307 for the residuals of the models with three (top) to six (bottom) periodic signals together with the analytic 10%, 1%, and 0.1% FAPs. 

Open with DEXTER  
In the text 
Fig. 3 Periodogram series of the signals detected in the FWHM, from most significant to less significant (top to bottom). 

Open with DEXTER  
In the text 
Fig. 4 Time series of the FWHM activity index. The solid black line represents the best fit to a model containing three sinusoids (periods of 5000, 23 and 1170 days). 

Open with DEXTER  
In the text 
Fig. 5 Periodogram series of the signals detected in the Sindex, from most significant to less significant (top to bottom). 

Open with DEXTER  
In the text 
Fig. 6 Timeseries of the Sindex. The solid black represents the best fit to a model containing three sinusoids (periods of 4000, 320, and 43 days). Longperiod trend and two signals at 320 and 43 days are clearly detected in the timeseries of the Sindex. 

Open with DEXTER  
In the text 
Fig. 7 Rms as a function of the bluest echelle aperture used. The corresponding λ_{BC} are illustrated at the top of the panel. A shallow minimum is reached at aperture 26 (λ_{BC} = 433 nm), indicating activity induced signals/jitter at the bluest wavelengths. 

Open with DEXTER  
In the text 
Fig. 8 Residual periodogram of the threeplanet solution as a function of the blue cutoff wavelength λ_{BC}. The top periodogram is obtained using the fullspectrum RVs. The bottom periodogram is obtained using the aperture 40 (505.0 nm) as λ_{BC}. All signals discussed in the text are marked with narrow vertical lines. The horizontal line represents the peak value of peak closer to 320 days. When λ_{BC} = 433 nm (central periodogram), the emerging 125day peak is already as high as the activityinduced signal. The 120day peak is an alias of the most likely period of 200 days for planet candidate d. 

Open with DEXTER  
In the text 
Fig. 9 Amplitudes of the four signals under discussion as a function of the blue cutoff λ_{BC}. While all proposed candidates keep a constant (or slightly increase) the derived semiamplitude, the 320day signal is almost nonexistent (the period actually becomes unconstrained) when wavelengths redder than 480 nm are used to derived a sevenplanet orbital solution. 

Open with DEXTER  
In the text 
Fig. 10 Residuals of the threeplanet fit when only the bluemost part of the spectra were used to extract the Doppler measurements. A quadratic trend similar to the one detected in the Sindex is the most obvious signal left in these residuals. 

Open with DEXTER  
In the text 
Fig. 11 As in Fig. 2 but for the unbinned radial velocities of the 490–680 nm spectrum. 

Open with DEXTER  
In the text 
Fig. 12 Phasefolded MAP Keplerian signals of the sixplanet solution with the other five signals and the MA(3) components removed from each panel. 

Open with DEXTER  
In the text 
Fig. 13 Distributions estimating the posterior densities of orbital periods (P_{x}), eccentricities (e_{x}), and radial velocity amplitudes (K_{x}) of the six Keplerian signals. The dashed lines show the prior densities for comparison for the orbital eccentricities. The solid curves are Gaussian densities with the same mean (μ) and variance (σ^{2}) as the parameter distributions. Additional statistics, mode, skewness (μ^{3}), and kurtosis (μ^{4}) of the distributions are also shown. 

Open with DEXTER  
In the text 
Fig. 14 Approximate Lagrange stability thresholds (shaded areas) between each two planets (red circles) and the MAP orbital parameters of the sixplanet solution (Table 4). 

Open with DEXTER  
In the text 
Fig. 15 Distributions estimating the orbital eccentricities given measurements alone (grey histogram; mean μ_{m} and standard deviation σ_{m}), showing the values the eccentricities had during orbital integrations of 10^{6} years that could not be shown to be unstable (white histogram; mean μ_{d} and standard deviation σ_{d}), and their combination (red histogram). 

Open with DEXTER  
In the text 
Fig. 16 Stability analysis of the circular sixplanet solution (S1). The panels show the influence of changing the initial eccentricity (left) and initial mean anomaly (right) on the stability of the complete planetary system. “Red” corresponds to highly unstable orbits while “black” denotes regular motion. Each of the 100 initial conditions per planet were integrated numerically for 1 kyr, for which we then computed the stability criterion D (Eq. (2)). 

Open with DEXTER  
In the text 
Fig. 17 Representation of the liquid water habitable zone around HD 40307 as described by Selsis et al. (2007) compared to our solar system. Note the compactness of the orbital configuration of the five inner planet candidates around HD 40307 compared to the planets in the inner solar system. 

Open with DEXTER  
In the text 