EDP Sciences
Free access
Issue
A&A
Volume 549, January 2013
Article Number A48
Number of page(s) 23
Section Planets and planetary systems
DOI http://dx.doi.org/10.1051/0004-6361/201220268
Published online 17 December 2012

© ESO, 2012

1. Introduction

Current high-precision 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 low-mass planets orbiting nearby stars. During recent years, radial velocity (RV) planet searches have revealed several systems of super-Earths and/or Neptune-mass 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 super-Earths 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 low-mass planets may be found in highly compact multiple systems that are still stable in long-term, 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 scaled-up versions of the Earths is not entirely clear (Barnes et al. 2009a). Their masses, between those of Earth and Neptune, suggest they are Neptune-like proto-gas 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 super-Earths 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 super-Earth mass objects can have rocky compositions.

In this article we re-analyse the 345 HARPS spectra publicly available through the ESO archive using a newly developed software tool called HARPS-TERRA (template-enhanced radial velocity re-analysis application; Anglada-Escudé & Butler 2012). Instead of the classic cross-correlation function method (CCF) implemented by the standard HARPS-ESO data reduction software (HARPS-DRS), we derive Doppler measurements by least-squares matching of each observed spectrum to a high signal-to-noise ratio template built from the same observations. A description of the method and the implementation details are given in Anglada-Escudé & 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 HARPS-DRS. As an example, it allows us to re-obtain 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- surface-related 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 self-consistent 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 (re-analysis of the spectra and Bayesian inference) allows the detection and verification of these low-amplitude 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 320-day signal is replaced by a signal of a super-Earth-mass candidate with a 200-day 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    Lstar / L =  −0.639 ± 0.060; Ghezzi et al. 2010) than the Sun. The star is quiescent (; Mayor et al. 2009a) and relatively metal-poor with  [Fe/H]  =  −0.31 ± 0.03 (Sousa et al. 2008). It also lacks massive planetary companions, which makes it an ideal target for high-precision RV surveys aiming at finding low-mass planets. According to the calibration of Barnes (2007), HD 40307 likely has an age similar to that of the Sun (~4.5 Gyr).

Table 1

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 time-scales 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 zero-mean and known variance (nominal uncertainties in the RVs) and another with zero-mean 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 mi as (1)where rk(ti) is the superposition of k Keplerian signals at epoch ti and γ is the reference velocity. The random variable ϵi is the Gaussian white noise component of the noise model with zero-mean 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 time-scale. Specifically, using hours as units of time, the exponential term (that is always  < 1 because ti > ti − 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 Metropolis-Hastings 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 one-block Metropolis-Hastings (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 k-Keplerian signals.

The prior probability densities in our analyses were essentially uniform densities. As in Tuomi (2012), we adopted the priors of the RV amplitude π(Ki) = U(0,aRV), reference velocity π(γ) = U( − aRV,aRV), and jitter π(σJ) = U(0,aRV), where U(a,b) denotes a uniform density in the interval  [a,b] . Since the observed peak-to-peak difference in the raw RVs is lower than 10 m s-1, the hyperparameter aRV was conservatively selected to have a value of 20 m s-1. The priors of the longitude of pericentre (ω) and the mean anomaly (M0) 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 scale-invariant parameter. The prior of this parameter was set uniform such that the two cut-off periodicities were Tmin and Tmax. These hyperparameters were selected as Tmin = 1.0 days and Tmax = 10Tobs 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 (Tobs), because signals in excess of that can be detected in RV data (Tuomi et al. 2009) and because there might be long-period 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 semi-Gaussian 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. time-shift 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 k-Keplerian signals is as follows. First, we require that the posterior probability of a k-Keplerian model is at least 150 times greater than that of the k − 1-Keplerian 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 well-constrained 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 least-squares periodograms (Lomb 1976; Scargle 1982) to probe the next most significant periods left in the data. In particular, we used the least-squares 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 F-ratio statistic (or power) of the fit. While strong powers likely indicate the existence of a periodic signal (though strong powers may be caused by sampling-related 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 multi-Keplerian fits due to strong correlations and aliases between clearly detected signals and yet-undetected lower amplitude companions (e.g. Anglada-Escudé 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 goodness-of-fit 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 false-alarm 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 long-period trend using our new RVs and we did not detect evidence for a trend in the new CCF RVs obtained using the HARPS-DRS. As shown in Anglada-Escudé & 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 HARPS-DRS v3.5 based on the release notes1 (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 [JD-2 450 000]. This feature can effectively decrease the phase-coverage 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 HARPS-TERRA 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 HARPS-TERRA 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 long-term 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 well-educated guesses based on the prior experience with RV data and reported stability of the instrument. Therefore, one must be especially careful not to over-interpret 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 re-analysed the nighly binned RVs to see whether we could independently reproduce the results of Mayor et al. (2009a) when HARPS-TERRA 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 three-Keplerian 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 three-Keplerian 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.

thumbnail Fig. 1

Least-squares 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 four-Keplerian solution was found to correspond to the three previously known super-Earth 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 four-Keplerian model was 1.9 × 105 times greater than that of a three-Keplerian one, making the 320-day signal significant. In addition to this signal, we could identify a 51-day 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 × 106. We could furthermore identify a sixth signal with our six-Keplerian model, correponding to a period of 34.4 days. However, even though the samplings converged well and the solution looked well-constrained, the six-Keplerian model was only five times more probable than the five-Keplerian 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 four-Keplerian model indicated that the global solution contained the 320-day periodicity as the fourth signal and yielded a posterior probability for the four-Keplerian model roughly 1.5  ×  106 times higher than for the three-Keplerian one. The nature of this signal and its relation to the stellar activity (Sect. 5) is discussed in Sect. 5.3.

thumbnail Fig. 2

Least-squares 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 four-Keplerian 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 × 105.

The residuals of the five-Keplerian model contained a periodicity at 34.7 days (Fig. 2, third panel) exceeding the 1% FAP level. The corresponding six-Keplerian model with this candidate received the highest posterior probability – roughly 5.0 × 108 times higher than that of the five-Keplerian model. Since the parameters of this sixth candidate were also well-constrained, we conclude that including the 34.7-day signal in the statistical model is fully justified by the data.

We also attempted to sample the parameter space of a seven-Keplerian model but failed to find a clear probability maximum for a seventh signal (see also the residual periodogram of the six-Keplerian model in Fig. 2, bottom panel). Although the periodicity of the seventh signal did not converge to a well-constrained probability maximum, all periodicities in the six-Keplerian model at 4.3, 9.6, 20.4, 34.7, 50.8, and 320 days were still well-constrained, 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 seven-Keplerian model. Therefore we stopped looking for additional signals.

From this analysis, we can state confidently that there are six significant periodicities in the HARPS-TERRA 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 (S-index) 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 HARPS-DRS. 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 so-called S-index in the Mount Wilson system (Baliunas et al. 1995), is automatically measured by HARPS-TERRA on the blaze-corrected one-dimensional spectra provided by the HARPS-DRS. The S-index 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 least-squares fitting of periodic signals that are each described by a sine-wave 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 long-term stability of the HARPS instrumental profile (and therefore its long-term 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).

thumbnail Fig. 3

Periodogram series of the signals detected in the FWHM, from most significant to less significant (top to bottom).

Open with DEXTER

thumbnail 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 S-index

The S-index also shows strong coherent variablity. The full set of 345 measurements of the S-index are provided in Table A.1. Again, this analysis was performed for the nightly binned measurements using least-squares periodograms and the sequential inclusion of sinusoidal signals. The last two S-index 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 S-index. From most significant to less significant, we found these to consist of a long-period 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 long-period 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.

thumbnail Fig. 5

Periodogram series of the signals detected in the S-index, from most significant to less significant (top to bottom).

Open with DEXTER

thumbnail Fig. 6

Time-series of the S-index. The solid black represents the best fit to a model containing three sinusoids (periods of 4000, 320, and 43 days). Long-period trend and two signals at 320 and 43 days are clearly detected in the time-series of the S-index.

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 320-day 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 S-index (43 days) appears to be a much better match to the expected stellar rotation of HD 40307 (Mayor et al. 2009b). The S-index 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 S-indices (long-period trend and 320-day 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 S-index 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, Anglada-Escudé & Butler (2012) showed that long-period 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 activity-induced 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. Wavelength-dependent signals

Anglada-Escudé & Butler (2012) showed that there is a correlation between the S-indices and RV measurements of HD 69830 (G8 V) and that an apparent long-period trend in the RVs completely disappeared when one obtains the velocities using the redmost half of the HARPS spectra only. Following Anglada-Escudé & 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 activity-related 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 H-band (~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 944-20. This candidate was first detected at optical wavelengths (2 km s-1 semi-amplitude) and then ruled out by the same group using nIR NIRSPEC/Keck observations (Martín et al. 2006). The chromatic nature of activity-induced 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 320-day signal (and all the others) to investigate its (their) nature in more detail.

As suggested by Anglada-Escudé & 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 cut-off 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, Anglada-Escudé & 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 cut-off 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 wavelength-dependent noise that acts stronger in the blue part of the spectrum.

thumbnail 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 activity-induced 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 320-day signal. Interestingly, when λBC is set to 430 nm (20th HARPS aperture), the 320-day peak becomes as high as a new one emerging around 120 days. When the λBC is set to 503 nm (40th HARPS aperture), the 120-day signal clearly dominates the periodogram. It is noteworthy that the powers of both the 34- and 51-day peaks remain mostly constant, which shows that their significance is not seriously affected by the choice of the cut-off wavelength.

thumbnail Fig. 8

Residual periodogram of the three-planet solution as a function of the blue cut-off wavelength λBC. The top periodogram is obtained using the full-spectrum 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 125-day peak is already as high as the activity-induced signal. The 120-day 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 120-day 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 320-day 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 least-squares 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 320-day signal almost completely vanishes when λBC is redder than 450 nm (K < 0.4 m s-1). It is also worth noting that the 200-day signal increases its amplitude as the 320-day signal disappears. The 200-day signal is seriously affected by yearly aliases which, effectively, correlates its amplitude with the signal at 320 days. The 320-day 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.

thumbnail Fig. 9

Amplitudes of the four signals under discussion as a function of the blue cut-off λBC. While all proposed candidates keep a constant (or slightly increase) the derived semi-amplitude, the 320-day signal is almost nonexistent (the period actually becomes unconstrained) when wavelengths redder than 480 nm are used to derived a seven-planet orbital solution.

Open with DEXTER

At this point, it is fair to ask why the long-term variability in the S-index (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 long-period trend also observed in the S-index (compare Figs. 6 and 10). No blue RV counterparts of the other tentative periodicities in the S-index were identified in the data.

thumbnail Fig. 10

Residuals of the three-planet 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 S-index 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  ~320-day 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 seven-planet 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 Anglada-Escudé & 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 S-index and the occurrence of the RV period at 320 days allows us to establish a very likely connection between this activity-related 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 re-analysis of the redmost half of the spectrum only (490–680 nm) in the next section.

Table 2

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 least-squares 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 least-squares 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 three-Keplerian 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 125-day period, the MCMC samplings showed that the 50.8-day periodicity is the global solution by providing higher posterior probabilities to the four-Keplerian model. We note that when calculating the residual periodogram of the three-Keplerian model, it is necessary to assume a fixed 3-Keplerian 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 least-squares sum). After several MCMC samplings starting in the vicinity of the different possible periods suggested by the periodogram, this 50.8-day signal was always found to correspond to the global four-Keplerian solution and made the four-Keplerian model 3.1 × 1012 times more probable than the three-Keplerian one.

thumbnail 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 five-Keplerian 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 six-Keplerian 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.6-day powers provided the global maxima and improved the posterior model probabilities by factors of 3.6 × 1011 and 2.5  ×  109 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 seven-Keplerian model did not converge to a clear probability maximum, thus failing to satisfy our requirement that the period of a signal has to be well-constrained for a positive detection.

Table 3

Relative posterior probabilities of models with k = 0,...,6 Keplerian signals (ℳk) and the periods (Ps) 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 seven-Keplerian model in Table 3 because the posterior samplings did not converge to a clear maximum (nor several maxima) in the seven-Keplerian 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 phase-folded MAP estimated orbits of the six Keplerian signals in Fig. 12 after removing the MA(3) components.

thumbnail Fig. 12

Phase-folded MAP Keplerian signals of the six-planet solution with the other five signals and the MA(3) components removed from each panel.

Open with DEXTER

Table 4

Six-planet solution of HD 40307 radial velocities.

thumbnail Fig. 13

Distributions estimating the posterior densities of orbital periods (Px), eccentricities (ex), and radial velocity amplitudes (Kx) 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 six-Keplerian 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 198-day signals having amplitudes with MAP estimates below 1.0 m s-1. Also, the periods of all signals were well-constrained and had clear MAP values and close-Gaussian densities (Table 4; Fig. 13). The distributions of orbital eccentricities are also shown in Fig. 13. In agreement with the quicklook detection using least-squares periodograms assuming circular orbits, these distributions demonstrate that all signals had eccentricities consistent with close-circular 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 mid-K dwarf in terms of Doppler information (higher signal-to-noise, 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 125-day alias was preferred against the 200-day signal for HD 40307 g. However, the Bayesian samplings including the complete model for the uncertainties still converged to the same six-planet solution (Pg ~ 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 200-day 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 34-day 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 small-sample 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 six-Keplerian model had the highest posterior probability and was 1.3 × 105 times more probable than the five-Keplerian 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 wavelength-dependent 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 third-order 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 long-term activity-related 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.

Table 5

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 long-term stability of the planetary system corresponding to our six-Keplerian 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 six-planet 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 close-circular orbits. Especially, the outer companion is unlikely to have orbtital eccentricity in excess of 0.4 (Fig. 14).

thumbnail Fig. 14

Approximate Lagrange stability thresholds (shaded areas) between each two planets (red circles) and the MAP orbital parameters of the six-planet solution (Table 4).

Open with DEXTER

thumbnail 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 106 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

thumbnail Fig. 16

Stability analysis of the circular six-planet 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 106 years, which is appropriate to exclude “almost all” unstable configurations (Barnes & Quinn 2004). In our integrations, we used the Bulirsch-Stoer 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 104 − 105 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 edyn. While this experiment shows that a non-negligible 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 six-planet solution given in Table 4, one quickly finds that the orbital evolution leads to collisional trajectories within a timescale of 103 years.

thumbnail 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 close-circular 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 Me and Mf. 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. Mf(t0) ≈ 310°) and also conditions that lead to chaotic motion (D ~ 10-2 for e.g. Mf(t0) ≈ 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 follow-up

The so-called liquid water habitable zone (HZ) is the area around a star where an Earth-like 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 semi-major 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 long-period 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 star-planet angular separation can be as large as  ~46 mas. Space-based direct-imaging missions such as Darwin/ESA (Cockell et al. 2009) or NASA/TPF (Lawson et al. 2006) aim to reach the sensitivity to detect Earth-sized planets at an inner working angle of  ~40 mas. Therefore, HD 40307 g would be the first habitable-zone candidate in the super-Earth 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 (S-index), we easily detect the three super-Earth 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 super-Earth candidates around this star (Table 3; Fig. 11) and our dynamical analyses show that long-term 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 Kepler-11 six-planet 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 CO2 or water would be required to sustain wet atmospheres. As is the case with Kepler-22 b (Borucki et al. 2012) and GJ 667C c (Anglada-Escudé 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, Kepler-22 b is likely to be a small-scale 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 scaled-down versions of Neptune-like planets with thick atmospheres and solid cores. Compared to the candidates orbiting nearby low-mass stars (GJ 581 d, GJ 667C c, HD 85512 b), HD 40307 g is not likely to suffer from tidal-locking – which improves its chances of hosting an Earth-like climate. Like for the other HZ candidates, it is not yet possible to determine its physical and geochemical properties and a direct-imaging 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 data-reduction techniques (e.g. HARPS-TERRA; Anglada-Escudé & 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 low-amplitude signals also introduces a practical problem. The orbits of the planets can no longer be readily compared with large numbers of data points using phase-folded 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 low-mass 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 life-hosting 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 Earth-mass 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. Anglada-Escudé 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/9-1. S. Vogt gratefully acknowledges support from NFS grant AST-0307493. 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 HARPS-ESO team in improving the instrument and its data reduction pipelines, and obtaining the observations that made this work possible.

References

Online material

Appendix A: Radial velocities and activity indices

Table A.1

HARPS-TERRA differential radial velocities of HD 40307 from the full spectra (F) and the red part (490–680 nm) of the spectra (R).

All Tables

Table 1

Stellar properties of HD 40307.

Table 2

Overview of the significant periodicities found in two activity indicators.

Table 3

Relative posterior probabilities of models with k = 0,...,6 Keplerian signals (ℳk) and the periods (Ps) of the signals added to the model when increasing the number of signals in the model by one.

Table 4

Six-planet solution of HD 40307 radial velocities.

Table 5

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.

Table A.1

HARPS-TERRA differential radial velocities of HD 40307 from the full spectra (F) and the red part (490–680 nm) of the spectra (R).

All Figures

thumbnail Fig. 1

Least-squares 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
thumbnail Fig. 2

Least-squares 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
thumbnail 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
thumbnail 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
thumbnail Fig. 5

Periodogram series of the signals detected in the S-index, from most significant to less significant (top to bottom).

Open with DEXTER
In the text
thumbnail Fig. 6

Time-series of the S-index. The solid black represents the best fit to a model containing three sinusoids (periods of 4000, 320, and 43 days). Long-period trend and two signals at 320 and 43 days are clearly detected in the time-series of the S-index.

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 8

Residual periodogram of the three-planet solution as a function of the blue cut-off wavelength λBC. The top periodogram is obtained using the full-spectrum 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 125-day peak is already as high as the activity-induced signal. The 120-day peak is an alias of the most likely period of 200 days for planet candidate d.

Open with DEXTER
In the text
thumbnail Fig. 9

Amplitudes of the four signals under discussion as a function of the blue cut-off λBC. While all proposed candidates keep a constant (or slightly increase) the derived semi-amplitude, the 320-day signal is almost nonexistent (the period actually becomes unconstrained) when wavelengths redder than 480 nm are used to derived a seven-planet orbital solution.

Open with DEXTER
In the text
thumbnail Fig. 10

Residuals of the three-planet 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 S-index is the most obvious signal left in these residuals.

Open with DEXTER
In the text
thumbnail Fig. 11

As in Fig. 2 but for the unbinned radial velocities of the 490–680 nm spectrum.

Open with DEXTER
In the text
thumbnail Fig. 12

Phase-folded MAP Keplerian signals of the six-planet solution with the other five signals and the MA(3) components removed from each panel.

Open with DEXTER
In the text
thumbnail Fig. 13

Distributions estimating the posterior densities of orbital periods (Px), eccentricities (ex), and radial velocity amplitudes (Kx) 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
thumbnail Fig. 14

Approximate Lagrange stability thresholds (shaded areas) between each two planets (red circles) and the MAP orbital parameters of the six-planet solution (Table 4).

Open with DEXTER
In the text
thumbnail 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 106 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
thumbnail Fig. 16

Stability analysis of the circular six-planet 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
thumbnail 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