Απόσπασμα:
Αρχική Δημοσίευση από SenseiG
«Το άστρο HD 40307 είναι ένα εξαιρετικά ήσυχο παλαιό άστρο-νάνος, όποτε δεν υπάρχει κάποιος λόγος ώστε ένας τέτοιος πλανήτης να μην μπορεί να διαθέτει «γήινη» ατμόσφαιρά», τόνισε o Anglad-Escude.
Η έρευνα των αστρονόμων δημοσιεύεται εδώ. |
Καλά που έβαλες και την πηγή,
γιατί διαφορετικά θα αμφισβητούσαμε την γνησιότητα του άρθρου,
του οποίου παραθέτω τον πρόλογο παρακάτω, για κάποιους δύσπιστους ... 
Habitable-zone super-Earth candidate in a six-planet system
around the K2.5V star HD 40307
Mikko Tuomi
⋆1,2, Guillem Anglada-Escud´e⋆⋆3, Enrico Gerlach4, Hugh R. A. Jones1, Ansgar Reiners3,
Eugenio J. Rivera
5, Steven S. Vogt5, and R. Paul Butler6
1
University of Hertfordshire, Centre for Astrophysics Research, Science and Technology Research Institute, College Lane, AL10
9AB, Hatfield, UK
2
University of Turku, Tuorla Observatory, Department of Physics and Astronomy, V¨ais¨al¨antie 20, FI-21500, Piikki¨o, Finland
3
Universit¨at G¨ottingen, Institut f¨ur Astrophysik, Friedrich-Hund-Platz 1, 37077 G¨ottingen, Germany
4
Lohrmann Observatory, Technical University Dresden, D-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 XX.XX.2012
/ Accepted XX.XX.XXXX
ABSTRACT
Context.
The K2.5 dwarf HD 40307 has been reported to host three super-Earths. The system lacks massive planets and is therefore
a potential candidate for having additional low-mass planetary companions.
Aims.
We re-derive 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 low-amplitude 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 di
ffering 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 S-index 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
M sin i ∼ 7 M⊕ (HD 40307 g).
Conclusions.
We show that Doppler measurements can be filtered for activity-induced 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 su
ffer from tidal locking. Also, at an angular separation of ∼ 46 mas, HD 40307 g
would be a primary target for a future space-based direct-imaging mission.
Key words.
Methods: Statistical, Numerical – Techniques: Radial velocities – Stars: Individual: HD 40307
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., 2011a; 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 lowmass
planets may be found in highly compact multiple systems
⋆
e-mail: guillem.anglada@gmail.com
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´eger
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´e& Butler, 2012).
Instead of the classic cross-correlation function method (CCF)
implemented by the standard HARPS-ESO data reduction software
(HARPS-DRS), we derive Dopplermeasurements by leastsquares
matching of each observed spectrum to a high signal-tonoise
ratio template built from the same observations. A description
of the method and the implementation details are given in
Anglada-Escud´e & Butler (2012). In addition to an increase in
precision (especially for K andMdwarfs), this method allows us
to performadditional 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 Section 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
e
ffects, 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 (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 (Section 2) and describe the statistical modelling
of the observations and the data analysis techniques we used
(Section 3). In Section 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 Section 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 (Section 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
∼ 7M⊕ orbiting
within the liquid water habitable zone of HD 40307. The analysis
of the dynamical stability of the system (Section 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 Section 8.We give some concluding remarks and
discuss the prospects of future work in Section 9.
K2.5 V star is a nearby dwarf with a Hipparcos parallax of 77.95
2. Stellar properties of HD 40307
We list the basic stellar properties of HD 40307 in Table 1. This
??
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 (log
R′
HK
< −4.99;
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
Table 1.
Stellar properties of HD 40307.
Parameter Estimate Reference
Spectral Type K2.5 V Gray et al. (2006)
log
R′
HK
-4.99 Mayor et al. (2009a)
π
[mas] 77.95??0.53 van Leeuwen (2007)
log L
star/L⊙ -0.639??0.060 Ghezzi et al. (2010)
log
g 4.47??0.16 Sousa et al. (2008)
M
star [M⊙] 0.77??0.05 Sousa et al. (2008)
T
eff [K] 4956??50 Ghezzi et al. (2010)
[Fe
/H] -0.31??0.03 Sousa et al. (2008)
v
sin i [kms−1] <1 Mayor et al. (2009a)
P
rot [days] ∼ 48 Mayor et al. (2009a)
Age [Gyr]
∼ 4.5 Barnes (2007)
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).
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 anMA 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 unknownmanner.We
limit our analysis to MA models of third order (MA(3) models)
because higher order choices did not improve the noise model
significantly. E
ffectively, 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 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 e
ffects 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
m
i
= rk (ti) + γ + ǫi +
p
X
j
=1
φ
jhmi−j − rk(ti−j) − γi exp ti−j − ti, (1)
where
rk(ti) 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
zero-mean and variance
σ2 = σ2i
+
σ2J
, where
σ2i
is the (fixed)
nominal uncertainty of the
ith measurements and σ2J
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 Equation (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
t
i
> ti −j) vanishes in few hours.
To demonstrate the impact of binning on this data set,
Section 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 algorithmis 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 di
fferent models. We used the oneblock
Metropolis-Hastings (OBMH) method (Chib & Jeliazkov,
2001; Clyde et al., 2007) to calculate the relative posterior probabilities
of models with di
ffering numbers of Keplerian signals
(e.g. Tuomi & Kotiranta, 2009; Tuomi, 2011; Tuomi et al., 2011;
Tuomi, 2012). We performed several samplings using di
fferent
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|Mk), where m is the measurement
vector and
Mk 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 di
fference in the raw RVs is lower
than 10 m s
−1, the hyperparameter aRV was conservatively selected
to have a value of 20 ms
−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-o
ff periodicities were Tmin and
T
max
. 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
π(ei) ∝ N(0, σ2e
) 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
themif 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 theMA 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 di
ffering numbers of Keplerian signals. Instead, following
Tuomi (2012), we set them as
P(Mk) = 2P(Mk+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 di
fferent from zero. Third, we also require that the
periods of all signals are well-constrained fromabove 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 o
ffset 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´e, 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 e
ffectively 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´e & Butler (2012) (Fig. 3), changes
in the continuum flux accross each echelle order (also called
blaze function) induce RV shifts of several ms
−1 if not properly
accounted for. This e
ffect was reported to affect HARPS
measurements in Pepe (2010) and appears to have been fixed
by HARPS-DRS 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−1yr−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-2450000].This feature can
e
ffectively 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 ˆ
ve 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 ms
−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
σ−1
e
= PNe
i
(σei
)
−1, 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 ms
−1 to account for instrumental and stellar
noise. After some tests, we found that adding 0.5ms
−1 in quadrature
to the uncertainties of the nightly averages ensured that none
of the epochs had uncertainties below the 0.7 ms
−1 level. The
typical uncertainties of a single night derived this way were of
the order of 0.8 ms
−1. These corrections are basically only welleducated
guesses based on the prior experiencewith RV data and
reported stability of the instrument. Therefore, onemust 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.
1
www.eso.org/sci/facilities/lasilla/instruments/harps/tools/drs.html
5
10
15
20
0
5
10
15
20
Power
0
5
10
15
20
Power
0
Power
10 100 1000
5
10
15
20
Period [days]
0
Power
Nightly binned RVs
10%
1%
0.1%
10%
10%
10%
1%
1%
1%
0.1%
0.1%
0.1%
51-d
34-d
320-d
4-th signal
5-th signal
6-th signal
7-th signal
200-d
200-d
200-d
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.
(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-Keplerianmodel 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.
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 Section 3.2. Both
of the two new significant signals had MAP estimates of their
radial velocity amplitudes slightly lower than 1.0 ms
−1 – the signals
at 51 and 320 days had amplitudes of 0.70 [0.31, 1.09] ms
−1
and 0.75 [0.38, 1.12] ms
−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
(Section 5) is discussed in Section 5.3.
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 wellconstrained
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
5
10
15
20
0
5
10
15
20
Power
0
5
10
15
20
Power
0
Power
10 100 1000
5
10
15
20
Period [days]
0
Power
Full spectrum RVs
10%
1%
0.1%
10%
10%
10%
1%
1%
1%
0.1%
0.1%
0.1%
51-d
34-d
320-d
4-th signal
5-th signal
6-th signal
7-th signal
200-d
200-d
200-d
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.
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.
probability maxima. Because we were unable to constrain the
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 fullwidth
at half-maximum (FWHM) of the CCF. These indices
monitor di
fferent 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 o
ffsets 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 theMountWilson 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., 2011a) 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 ms−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 ms
−1 yr−1. The 345 measurements
of the FWHM are listed in Table A.2. 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 e
ffects, 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. (2011a).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).
5.2. Analysis of the S-index
of 345 measurements of the S-index are provided in Table A.2.
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
The S-index also shows strong coherent variablity. The full set
5
10
15
20
25
30
0
5
10
15
20
25
30
Power
0
5
10
15
20
25
30
Power
0
Power
10 100 1000
5
10
15
20
25
30
Period [days]
0
Power
FWHM
10%
1%
0.1%
10%
10%
10%
1%
1%
1%
0.1%
0.1%
0.1%
5000+ d
23-d
1170-d
1st signal
2nd signal
3rd signal
4th signal
Fig. 3.
Periodogram series of the signals detected in the FWHM, from
most significant to less significant (top to bottom).
53000 53500 54000 54500 55000
-20
-10
0
10
20
30
JD-2400000 [days]
-30
FWHM - 5910 [m s
-1]
Fig. 4.
Habitable-zone super-Earth candidate in a six-planet system
around the K2.5V star HD 40307
Time series of the FWHM activity index. The solid black line
represents the best fit to a model containing three sinusoids (periods of
Mikko Tuomi
⋆1,2, Guillem Anglada-Escud´e⋆⋆3, Enrico Gerlach4, Hugh R. A. Jones1, Ansgar Reiners3,
Eugenio J. Rivera5, Steven S. Vogt5, and R. Paul Butler6
1
University of Hertfordshire, Centre for Astrophysics Research, Science and Technology Research Institute, College Lane, AL10
9AB, Hatfield, UK
2
University of Turku, Tuorla Observatory, Department of Physics and Astronomy, V¨ais¨al¨antie 20, FI-21500, Piikki¨o, Finland
3
Universit¨at G¨ottingen, Institut f¨ur Astrophysik, Friedrich-Hund-Platz 1, 37077 G¨ottingen, Germany
4
Lohrmann Observatory, Technical University Dresden, D-01062 Dresden, Germany
5
UCO/Lick Observatory, Department of Astronomy and Astrophysics, University of California at Santa Cruz, Santa Cruz, CA 95064,
USA
ABSTRACT
Context.
Aims.
Methods.
6
Department of Terrestrial Magnetism, Carnegie Institute of Washington, Washington, DC 20015, USA
Received XX.XX.2012 / Accepted XX.XX.XXXX
The K2.5 dwarf HD 40307 has been reported to host three super-Earths. The system lacks massive planets and is therefore
a potential candidate for having additional low-mass planetary companions.
We re-derive 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 low-amplitude signals.
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 S-index are strongly wavelength dependent.
Results.
Conclusions.
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 M sin i ∼ 7 M⊕ (HD 40307 g).
We show that Doppler measurements can be filtered for activity-induced 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 su ffer from tidal locking. Also, at an angular separation of ∼ 46 mas, HD 40307 g
would be a primary target for a future space-based direct-imaging mission.
Key words.
1. Introduction
Methods: Statistical, Numerical – Techniques: Radial velocities – Stars: Individual: HD 40307
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., 2011a; 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 lowmass
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´eger
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´e& Butler, 2012).
Instead of the classic cross-correlation function method (CCF)
implemented by the standard HARPS-ESO data reduction software
(HARPS-DRS), we derive Dopplermeasurements by leastsquares
matching of each observed spectrum to a high signal-tonoise
ratio template built from the same observations. A description
of the method and the implementation details are given in
Anglada-Escud´e & Butler (2012). In addition to an increase in
precision (especially for K andMdwarfs), this method allows us
to performadditional 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 Section 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
e ffects, 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 (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 (Section 2) and describe the statistical modelling
of the observations and the data analysis techniques we used
(Section 3). In Section 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 Section 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 (Section 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 ∼ 7M⊕ orbiting
within the liquid water habitable zone of HD 40307. The analysis
of the dynamical stability of the system (Section 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 Section 8.We give some concluding remarks and
discuss the prospects of future work in Section 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 (log R′
HK
< −4.99;
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
Table 1.
Stellar properties of HD 40307.
Parameter Estimate Reference
Spectral Type K2.5 V Gray et al. (2006)
log R′
HK
-4.99 Mayor et al. (2009a)
π
v
P
[mas] 77.95??0.53 van Leeuwen (2007)
log Lstar/L⊙ -0.639??0.060 Ghezzi et al. (2010)
log g 4.47??0.16 Sousa et al. (2008)
M star [M⊙] 0.77??0.05 Sousa et al. (2008)
T eff [K] 4956??50 Ghezzi et al. (2010)
[Fe /H] -0.31??0.03 Sousa et al. (2008)
sin i [kms−1] <1 Mayor et al. (2009a)
rot [days] ∼ 48 Mayor et al. (2009a)
Age [Gyr] ∼ 4.5 Barnes (2007)
3. Statistical analyses
3.1. Statistical models
m
i
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).
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 anMA 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 unknownmanner.We
limit our analysis to MA models of third order (MA(3) models)
because higher order choices did not improve the noise model
significantly. E
ffectively, 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 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 e ffects 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
= rk(ti) + γ + ǫi +
p
X
j
=1
φ
jhmi−j − rk(ti−j) − γi exp ti−j − ti, (1)
where rk(ti) 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
zero-mean and variance σ2 = σ2i
σ2J
, where
σ2i
is the (fixed)
nominal uncertainty of the
ith measurements and σ2J
is a free
parameter describing the magnitude of the jitter component.
Finally, the free parameters of the MA(
p) model are denoted as
φ
t
i
j, j = 1, ..., p – they describe the amount of correlation between
the noise of the ith and i − jth measurements. The exponential
term in Equation (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−j) vanishes in few hours.
To demonstrate the impact of binning on this data set,
Section 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 algorithmis 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 di fferent models. We used the oneblock
Metropolis-Hastings (OBMH) method (Chib & Jeliazkov,
2001; Clyde et al., 2007) to calculate the relative posterior probabilities
of models with di ffering numbers of Keplerian signals
(e.g. Tuomi & Kotiranta, 2009; Tuomi, 2011; Tuomi et al., 2011;
Tuomi, 2012). We performed several samplings using di fferent
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|Mk), where m is the measurement
vector and Mk 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
(γ) = U(−aRV, aRV), and jitter π(σJ) = U(0, aRV), where
(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 ms −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-o ff periodicities were Tmin and
T
max
. 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 π(ei) ∝ N(0, σ2e
) 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
themif they explain the data better than lower ones (Tuomi,
2012; Tuomi et al., 2012a).
For the MA components φj, we selected uniform priors as
π
3.3. Periodogram analysis
(φj) = U(−1, 1), for all j = 1, 2, 3. This choice was made to
ensure that theMA 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 di ffering numbers of Keplerian signals. Instead, following
Tuomi (2012), we set them as P(Mk) = 2P(Mk+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 di fferent from zero. Third, we also require that the
periods of all signals are well-constrained fromabove 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).
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 o
ffset 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´e, 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 e ffectively 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´e & Butler (2012) (Fig. 3), changes
in the continuum flux accross each echelle order (also called
blaze function) induce RV shifts of several ms −1 if not properly
accounted for. This e ffect was reported to affect HARPS
measurements in Pepe (2010) and appears to have been fixed
by HARPS-DRS 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−1yr−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-2450000].This feature can
e ffectively 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
e
i
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 ˆ
ve 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 ms −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 σ−1
= PNe
(σei
)
−1, 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 ms −1 to account for instrumental and stellar
noise. After some tests, we found that adding 0.5ms −1 in quadrature
to the uncertainties of the nightly averages ensured that none
of the epochs had uncertainties below the 0.7 ms −1 level. The
typical uncertainties of a single night derived this way were of
the order of 0.8 ms −1. These corrections are basically only welleducated
guesses based on the prior experiencewith RV data and
reported stability of the instrument. Therefore, onemust 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.
1
www.eso.org/sci/facilities/lasilla/instruments/harps/tools/drs.html
[www.eso.org]
0
5
10
15
20
Power
0
5
10
15
20
Power
0
5
10
15
20
Power
10 100 1000
Period [days]
0
5
10
15
20
Power
Nightly binned RVs
10%
1%
0.1%
10%
10%
10%
1%
1%
1%
0.1%
0.1%
0.1%
320-d
51-d
34-d
4-th signal
5-th signal
6-th signal
7-th signal
200-d
200-d
200-d
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.
(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-Keplerianmodel 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.
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 Section 3.2. Both
of the two new significant signals had MAP estimates of their
radial velocity amplitudes slightly lower than 1.0 ms −1 – the signals
at 51 and 320 days had amplitudes of 0.70 [0.31, 1.09] ms −1
and 0.75 [0.38, 1.12] ms
−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
(Section 5) is discussed in Section 5.3.
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
0
5
10
15
20
Power
0
5
10
15
20
Power
0
5
10
15
20
Power
10 100 1000
Period [days]
0
5
10
15
20
Power
Full spectrum RVs
10%
1%
0.1%
10%
10%
10%
1%
1%
1%
0.1%
0.1%
0.1%
320-d
51-d
34-d
4-th signal
5-th signal
6-th signal
7-th signal
200-d
200-d
200-d
Fig. 2.
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
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 wellconstrained
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
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.
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 fullwidth
at half-maximum (FWHM) of the CCF. These indices
monitor di
fferent 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 o ffsets 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 theMountWilson 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., 2011a) 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 ms−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 ms −1 yr−1. The 345 measurements
of the FWHM are listed in Table A.2. 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 e ffects, 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. (2011a).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).
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.2.
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
0
5
10
15
20
25
30
Power
0
5
10
15
20
25
30
Power
0
5
10
15
20
25
30
Power
10 100 1000
Period [days]
0
5
10
15
20
25
30
Power
FWHM
10%
1%
0.1%
10%
10%
10%
1%
1%
1%
0.1%
0.1%
0.1%
5000+ d
23-d
1170-d
1st signal
2nd signal
3rd signal
4th signal
Fig. 3.
Periodogram series of the signals detected in the FWHM, from
most significant to less significant (top to bottom).
53000 53500 54000 54500 55000
JD-2400000 [days]
-30
-20
-10
0
10
20
30
FWHM - 5910 [m s
-1]
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
.