© ESO, 2013

1. Introduction

Observing the periodic Doppler signatures of planetary companions in the stellar spectra is currently the most efficient method of finding exoplanets around nearby stars. This radial velocity (RV) method can be readily applied to M dwarfs that are the most populous stars in the Solar neighbourhood and are known to be hosts to several planetary systems, e.g. GJ 581 (Bonfils et al. 2005; Udry et al. 2007; Mayor et al. 2009), GJ 667C (Anglada-Escudé et al. 2012; Delfosse et al. 2013; Bonfils et al. 2013), GJ 676A (Forveille et al. 2011; Anglada-Escudé & Tuomi 2012), and GJ 876 (Delfosse et al. 1998; Marcy et al. 1998, 2001; Rivera et al. 2005, 2010). Because of the several on-going high-precision RV surveys, this population of exoplanet systems with multiple planets can readily be expected to increase as the observational baselines extend and as the amount of high-precision data increases and enables the detections of signals with lower and lower amplitudes.

In this article, we apply Bayesian multi-planet detection techniques (e.g. Tuomi et al. 2013a,b) to analyse new RV measurents obtained with the HARPS-TERRA software (Anglada-Escudé & Butler 2012) from public HARPS spectra of the nearby M dwarf GJ 163¹. Since a number of subjective choices are required by any Bayesian statistical data analyses, we discuss how this subjectivity – i.e. prior probability densities or prior models – affects our detection sensitivity and reliability. In particular, we analyse the same dataset by assuming different prior choices and noise models in order to obtain as objective results as possible.

Since the nature of RV noise of main sequence stars is not very well known and some stars show clear autocorrelation in their velocity noise (Tuomi et al. 2013b; Tuomi & Jenkins 2012; Baluev 2012), we analyse the GJ 163 velocities using different noise models, i.e. Gaussian white noise and a simple red noise model. We perform the analyses using posterior samplings (following Tuomi 2012; Tuomi & Jones 2012; Tuomi et al. 2013a) and calculate estimates for the Bayesian evidences using the truncated posterior mixture estimate (Tuomi & Jones 2012). As a result, we report our best estimates for the number and properties of planetary signals in the RV data of GJ 163 using the different prior and noise model choices. As a final validation procedure, we also investigate whether these signals correspond to any clear periodic features in the activity data of the star.

Prior probability densities are an integral part of any statistical analyses based on the Bayes’ rule of conditional probabilities. However, most, if not all, frequentist methods can in fact be derived from Bayesian ones by making certain simplifying assumptions on the shapes and natures of the prior probability densities and likelihood functions of the model parameters and measurements. For instance, all the maximum likelihood methods (e.g., the popular multi-variate least-squares minimisation methods) can be derived from the Bayesian ones by adopting uniform prior probability densities for the model parameters over some suitable range. Similarly, all the confidence level tests based on χ² statistics can be derived from the Bayes’ rule and Bayesian evidence ratio tests by assuming Gaussian likelihood functions and uniform prior densities. Because of its generality, only Bayesian framework enables to study the effects these subjective choices might have on the obtained results. All statistical analysis methods contain prior information (uniform prior is a prior as well) but only Bayesian ones enable taking them into account in a logically consistent manner.

When searching for planetary signatures in noisy data using the Doppler method, the choice of prior probability densities has not been discussed very extensively in the literature (some examples can be found in Baluev 2012; Tuomi & Jones 2012). Ford & Gregory (2007) defined a set of rather uninformative prior densities for analyses of radial velocity data using the typical statistical models consisting of Keplerian signals and white noise. These priors were modified slightly in Tuomi (2012) and Tuomi et al. (2013a) but the effects different prior choices have on the obtained results have not been studied. We study these effects in this article by comparing the effects of prior models on the obtained results when analysing the RV data of GJ 163. We discuss and motivate our prior choices in Sect. 2 and present the full analysis results of GJ 163 HARPS-TERRA velocities in Sect. 5.

2. Prior choice

The prior probability densities of model parameters represent the other component of Bayesian statistical tools in addition to the likelihood functions. The need to compare different likelihood functions, or likelihood models, is understood to be an important feature of statistical analyses of noisy data when the attempt is to find out the processes producing the observations, but the prior probability densities play a significant role as well and different prior models warrant comparisons in order to find the most trustworthy descriptions for the data.

While it is typically the case that the likelihood “overwhelms” priors in the sense that the likelihood function sets much stronger constraints to the posterior density than the prior, this is not necessarily such a good idea when there are not many more measurements than free parameters and when the noisy data does not constrain the parameters much. In such cases, it might be necessary to use stronger prior constraints, i.e. informative priors, to obtain any sensible solutions at all. Uninformative prior choices are generally preferred because such prior beliefs have as little effect as possible on the obtained results. In general, this is important because the relevance of a result is usually strongly tied to how constraining the prior choices are in practice, i.e. what are the assumptions (implicit or explicit) enabling one to reach a significant conclusion. In the next subsection, we will show that even uninformative priors can have, in principle, considerable effects on the obtained results and the inferred conclusions.

2.1. Uninformative priors

When performing statistical analyses of noisy data, the goal is typically to obtain as objective results as possible. In practice, the prior beliefs of a scientist performing data analyses should not be allowed to affect the obtained results – especially if they differ radically from the prior beliefs of fellow scientists. Priors that all scientists can consider rather objective are usually called uninformative priors because they aim at describing a maximum amount of ignorance on the system that the data is assumed to describe.

Uniform probability densities are a common example of such uninformative priors. Generally, for parameter θ, they can be written as π(θ) = c for all θ ∈ Ω, where Ω is the parameter space that is usually assumed to be a bounded subset of the real line to make the density a proper probability density and c some positive constant. However, because analysis results cannot depend on the chosen unitary system and should remain independent after any linear transformation of the parameters, it is always possible to choose c = 1.

Uniform priors cannot be considered uninformative for all the parameters, especially so, when searching for periodic signals in noisy data corresponding to signatures of extra-solar planets. For instance, it is unrealistical to assume a priori that an RV data set has an equal probability of containing a periodic signal of planetary origin between 0–100 days as between 100–200 days because the former contains much more orbits that are stable than the latter. Therefore, period (P) is a scale-parameter and an invariant prior for such a parameter is the Jeffreys’ prior that can be written as π(P) ∝ P^-1. This prior corresponds to a uniform prior probability density in the log-period space, which is the rationale behind the common choice of log P as a parameter in the statistical model instead of P.

2.2. Transformation of parameters

Under a nonlinear transformation from parameter vector θ to vector θ′, the corresponding changes in prior densities have to be accounted for in order to obtain consistent results using both parameterisations. Generally, a transformation θ → θ′ results in a change in the prior probability density defined as $\mathit{\pi }\mathrm{\left(}{\mathit{\theta }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\mathit{\pi }\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\left|\frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}{\mathit{\theta }}^{\mathrm{\prime }}}\right|\mathit{,}$(1)given that the transformation θ′ = f(θ) has an inverse transformation and $\left|\frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}{\mathit{\theta }}^{\mathrm{\prime }}}\right|$ is the Jacobian of this transformation. It is easy to see that a linear transformation of the form θ′ = aθ + b for π(θ) = c yields π(θ′) = c′. This means that under a linear transformation, i.e. a change of unit system, constant priors remain unchanged. However, this is not the case in general, which means that for any non-linear transformation between θ and θ′, the results might be rather different if uniform priors were used in both parameterisations.

2.3. Candidate priors

In this subsection, we define the candidate prior probability densities we use to analyse the velocities of GJ 163.

As a set of reference priors, we use those presented in Tuomi (2012) because this choice is a rather uninformative one and does not constrain the model parameters strongly (e.g. Anglada-Escudé & Tuomi 2012; Tuomi 2012; Tuomi et al. 2013a). The functional forms of the priors and the limits of the respective parameter spaces of the model parameters are shown in Table 1. Because all the parameters are real numbers, we denote the parameter spaces as intervals of the real line. The parameters of the Keplerian signals in Table 1 are the radial velocity amplitude (K), longitude of pericentre (ω), orbital eccentricity (e), mean anomaly (M₀), and the logarithm of the orbital period (l = log P). In addition to these Keplerian parameters, we also include a constant reference velocity (γ), a radial velocity jitter (σ_J; amount of excess white noise in the data on top of the estimated instrument uncertainties), and a correlation coefficient between the noise of subsequent epochs (φ) as defined in e.g. Tuomi et al. (2013a) and Tuomi et al. (2013b).

We choose the hyperparameters of each prior distribution in Table 1 as follows.

2.3.1. Semi-amplitude

The maximum semi-amplitude is set to K_max = 20 ms^-1. This choice is motivated by the fact that the standard deviation of the RVs is 6.8 ms^-1 and, therefore, we do not expect to find signals that have amplitudes in excess of this maximum amplitude.

2.3.2. Eccentricity

One of the most disputed and controversial priors choices in the statistical analysis of the Keplerian problem is the prior choice for the eccentricity. In this study, we consider three cases by using three choices of the hyperparameter σ_e (Table 1). As in previous works (e.g. Tuomi 2012), our reference choice is σ_e = 0.3 but we also obtain results with values of 0.2 and 0.1.

The low value of σ_e = 0.1 can be justified in the following way. The hyperparameter σ_e controls how we expect the probability of the eccentricity to decrease as it approaches unity. Our reference prior might not be of very practical use when the RV data-set is relatively small, it has uncertainties comparable to magnitude to the Keplerian signals it might contain, and its sampling cadence is very uneven (e.g. long gaps). In such cases, too uninformative assumptions, i.e. too high values of σ_e, would allow posterior probability densities to have roughly equally high values over large subsets of the parameter space making the interpretation of the obtained results very difficult, if possible at all. In the Keplerian problems, it is a common trick to restrict the parameter space of low amplitude candidates to strictly circular orbits (e.g. Lovis et al. 2011; Pepe et al. 2011; Dumusque et al. 2012), which decreases the number of free parameters in the model by 2 × N_planets. This assumption also avoids the strong non-linear regime of highly eccentric orbits where the implicit assumptions of point estimate methods (e.g., multivariate Gaussian posteriors) do not hold anymore. This choice corresponds to a delta-function prior density of the form π(e) = δ(e) that gives all the prebability density to e = 0. We considered this choice to be too limiting and selected a slightly less informative one, i.e. $\mathrm{𝒩}\mathrm{\left(}\mathrm{0}\mathit{,}{\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{2}}\mathrm{\right)}$ with σ_e = 0.1. This hyperparameter results in a prior density that is much narrower than the observed eccentricity distribution of low-mass planets (Fig. 1) but not as limiting as an eccentricity fixed to zero.

To investigate how the current statistical properties of the observed exoplanet candidates support our choice of the eccentricity prior, we obtained data from The Extrasolar Planets Encyclopaedia² and selected a sub-sample of all listed exoplanets with minimum masses lower than 0.1 M_Jup. As obtained in Dec. 6th, 2012, this sample contained 113 planet candidates for which eccentricity has been estimated. For 22 of them the reported estimate was exactly zero, which suggests that their eccentricities have been fixed to zero in the statistical analyses, or that only an upper limit was provided at the time of publication. Nevertheless, we assume that the eccentricities of these 22 candidates are indeed close to zero and put them in the [0 − 0.1] bin in Fig. 1. The obtained eccentricity distribution is illustrated in Fig. 1 where we also add a representation of our three test cases (eccentricity priors with σ_e = 0.1, 0.2 and 0.3). This figure shows that – for the sample of low-mass planets – the eccentricity prior with σ_e = 0.2 appears to be a reasonably good representation of the observed population, and therefore, it also poses an interesting choice for a prior density to be investigated.

Fig. 1 Eccentricity distribution (red circles) of low-mass (mp < 0.1  MJup) exoplanet candidates and the relative eccentricity priors with σe = 0.3 (solid line), σe = 0.2 (dash-dotted line), and σe = 0.1 (dashed line). The priors have been scaled to have maxima of 0.37 to enable visual comparison with the distribution.

We note that the eccentricities in the literature are likely overestimated in general for the simple reason that eccentricity is a positive number and cannot be lower than zero. This means that if the eccentricity estimate of a low-eccentricity orbit is biased, it it necessarily more likely to be overestimated than underestimated because the former would be more likely than the latter (Zakamska et al. 2011). Eccentricities are also overestimated systematically when the signal amplitude is not much larger than the measurement noise and when the number of observations w.r.t. the model parameters is low (O’Toole et al. 2008; Shen & Turner 2008). Therefore, our eccentricity prior might be a reasonable choice even for σ_e = 0.1.

2.3.3. Jitter

Because we might know that the target star is a quiescent and a rather inactive one (as is the case of GJ 163, see Sect. 4), we assign an informative prior density to the stellar jitter. A fixed jitter parameter is sometimes chosen based on reasonable, yet subjective, considerations such as maximum demonstrated precision of an instrument, or observed typical variability of other similar stars without planets. As for the eccentricity prior on the circular orbital case, it is a common practice to assume a rather restrictive and informative δ-function prior density for the jitter, i.e. to fix the jitter parameter to some a priori estimated value. To relax this condition and still be able to use our prior knowledge of the sample, we chose $\mathit{\pi }\mathrm{\left(}{\mathit{\sigma }}_{\mathrm{J}}\mathrm{\right)}\mathrm{\propto }\mathrm{𝒩}\mathrm{\left(}{\mathit{\mu }}_{\mathit{\sigma }}\mathit{,}{\mathit{\sigma }}_{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\right)}$ and selected the hyperparameters as μ_σ = 1 ms^-1 and σ_σ =1 ms^-1. These estimates are based on a sample of 27 M dwarfs targeted by HARPS that was found to have on average 1.9 ms^-1 variation (after removing the signals of all known planet candidates) and only two of them had variations in excess of 4 ms^-1 (da Silva et al. 2012). Because effectively we add this parameter describing the stellar jitter in quadrature to the estimated instrument uncertainties (Table 3) in our model, we consider that our prior choice represents the observed properties of M dwarfs reasonably well.

2.3.4. Period

The minimum and maximum periods in the parameter space are chosen as P₀ = 1 day and P_max = 2T_obs, where 1 day is the typical separation between consecutive observations and T_obs is the baseline of the observations³.

Especially when dealing with periodic signals, the period is the most important parameter in terms of confidently discovering a signal in the first place. However, because the parameter we actually use in our posterior samplings is not the period as such but its logarithm, we must note that the actual prior in the period space is then obtained by using the Jacobian of the transformation (Eq. (1)). After simple algebra it can be seen that our reference prior choice (uniform in log P) corresponds to an implicit prior of the period parameter of the form π(P) ∝ P^-1. This type of dependence might not seem to be very satisfactory one as a prior because it assigns higher probabilities for shorter periods than for longer ones. However, it corresponds to a scale-invariant prior, i.e. the Jeffreys’ prior, that remains the same regardless of the unit system of the period parameter. A uniform prior in frequency corresponds to an implicit prior of π(P) ∝ P^-2 and is essentially the one used when analysing time series in the frequency space using periodogram based methods (Baluev 2012). Therefore, we compare the three prior choices that seem natural to the problem: uniform in log P, uniform in P, and uniform in P^-1.

As already mentioned, our sampling strategy uses the logarithm of the orbital period as a parameter in the posterior samplings because it is a scale-invariant parameter. Therefore, we must express the priors uniform in P and in P^-1 in the log-orbital period space. Simple application of Eq. (1) shows that for π(P) = 1 it follows that π(l) ∝ exp(l). Similarly, for π(P^-1) = 1 it follows that π(l) ∝ exp( − l). Clearly, the differences between these priors and a uniform prior in l are considerable and certainly can be expected to have an effect on the obtained results. Yet, without using physical constrains, it cannot be said that any of them is more wrong than the others. While our reference prior of choice is the uniform one in log P, we study the effect of the alternative ones on the interpretation of the results obtained from the velocity data in Sect. 5.

3. GJ 163

GJ 163 is a nearby M3.5 V dwarf (Koen et al. 2010) with a Hipparcos parallax of 66.69 ± 1.82 mas (van Leeuwen 2007), which implies a distance of ~14.9 ± 0.4 pc. We use this distance with J, H and K photometry from 2MASS (Cutri et al. 2003) to derive a mass of 0.40 ± 0.02  M_⊙ using the mass-luminosity relation given by Delfosse et al. (2000). Since the uncertainty in the distance is rather low, the estimate of the mass is mostly limited by the precision of the calibration (~5% for M_⋆ < 0.5  M_⊙). Using the photometric metallicity calibrations from Johnson & Apps (2009) and Schlaufman & Laughlin (2010) and the V and K magntidues of the star, we obtain a slightly super-solar metallicity of [Fe/H] = 0.1 ± 0.1 (mean of the two calibrations), where the uncertainty is again intrinsic uncertainty in the empirical relations. This is the metallicity range where the evolutionary models produce a better agreement with observations. Using Chabrier & Baraffe (1997) and assuming an age between 2 and 10 Gyr, we obtain a luminosity of 0.0196  L_⊙ and a corresponding effective temperature of 3500 ± 100 K.

The parallax and proper motion of GJ 163 imply a UVW velocity vector of (− 70.1, − 75.5, 0.51) km s^-1 which has rather high components for a typical thin disk star, meaning that the star is more likely a member of the thick disk (see Fig. 3 in Bensby et al. 2003). Compared to other M dwarfs of similar spectral class, GJ 163 has a rather low S-index (0.61) – between the very stable and planet prolific GJ 581 (0.46) and the slightly more “jittery” but also planet prolific GJ 876 (0.82). The relative heights of the CaII K lines of some typical M dwarfs are shown in Fig. 2. This relatively low emission in CaII indicates that GJ 163 has low activity levels and suggests that it is rather old (>2 Gyr).

Fig. 2 CaII H line of GJ 163 (purple) as obtained from co-adding the 55 available HARPS spectra. The same lines of stars with similar spectral type are shown for comparison (offset in wavelength). The height of the line was scaled to the mean number of counts per pixel in this particular HARPS echelle order. Spectral types were taken from Bonfils et al. (2013). As discussed in the text, GJ 163 is among the least chromospherically active stars in its spectral range.

Fig. 3 HARPS-TERRA radial velocities of GJ 163 with the data median removed.

Given the precision in the distance and in the photometric colours, the estimation of the total luminosity we provide is dominated by systematic model uncertainties at the level of ~5% (Boyajian et al. 2012). While empirical calibrations based on interferometric radius measurements of some M dwarfs exist, they still result in large scatter (~20–30%) at effective temperatures below 4000 K (Boyajian et al. 2012). Therefore, and besides all caveats associated with using models, we believe that Chabrier & Baraffe (1997) provides a more reliable estimate of the stellar luminosity and we assign a 5% uncertainty to it.

3.1. Observations

The radial velocities underlying the current work are based on spectra obtained with the HARPS spectrograph during different observing programs over the recent years⁴. Visual inspection of the RV measurements shown in Fig. 3 shows that the star shows Doppler variability in excess of the typical uncertainties (~1.1 m s^-1) from the very first observations (obtained in 2003). A more intensive set of high-cadence observations was obtained in 2009 in the context of the HARPS-Totems program (PI. Bonfils) whose aim was to detect short period super-Earth- and Neptune-mass planet candidates with high probability of transit. According to the HARPS-ESO archive⁵, the star is still being monitored but observations made after August 2010 have not been made available.

The HARPS spectra are taken with typical exposure times of 900 seconds and the median signal-to-noise ratio (S/N) per pixel at 6100 Å is ~35. While this might not seem very high, spectra of M-dwarfs contain considerable amounts of Doppler information in the redder end of HARPS. This information is optimally extracted by the HARPS-TERRA software which consist on matching the extracted spectrum as provided by the HARPS data reduction software (DRS) to a very high S/N template made by coadding all the observations. The HARPS-DRS spectra are generated from a standardarized calibration set of observations that include dark, flat fielding and calibration lamp exposures (Th-Ar lamp) obtained at the beginning of each night. These sets of calibrations together with a carefully designed calibration plan have ensured the wavelength solution consistency down to a precision of 1 ms^-1 or better over the years (e.g. Lovis & Pepe 2007). Given that there were a total of 55 spectra available, the combined template spectrum has a S/N of ~250 at 6100 Å contributing negligibly to the error budget of each individual measurement. The high precision achievable with the HARPS-TERRA processing has been demonstrated in Anglada-Escudé & Butler (2012) and its comparisons to the HARPS-DRS derived velocities (obtained via cross correlation with a weighted binary mask) have generally shown that the HARPS-TERRA velocities have greater precision for M dwarfs (Anglada-Escudé et al. 2012, 2013; Anglada-Escudé & Tuomi 2012). In these studies, it has also been shown that some of the HARPS M-dwarf sample stars are as stable (even more stable, with rms < 1 m s^-1 over several years) than some other primary targets from the HARPS-GTO sample (mid-K to late G-dwarfs, such as HD 88512 or Tau Ceti Pepe et al. 2011). The Doppler time-series, as well as measurements of activity indices, are given in Table 3 (see Lovis et al. 2011; Anglada-Escudé & Butler 2012, for further details).

4. Signals in activity indices

In addition to random noise, stellar activity can also generate apparent Doppler periodicities that can be confused with true Keplerian signals. Before going into the analysis of the Doppler data, we first investigate whether there are periodic variations in the three representative activity indices of GJ 163 (see e.g. Anglada-Escudé & Tuomi 2012). These indices are the S-index, the line bisector (BIS) and the full-with at half-maximum (FWHM) of the cross correlation function.

The S-index is proportional to flux coming from the chromosphere of the star. It is measured as the amount of emission at the Ca II H+K doublet (393.3664 nm for the K and 396.8470 nm for the H line) compared to a locally defined continuum. The recipes for computing the S-index from HARPS spectra are given in Lovis et al. (2011) and Anglada-Escudé & Butler (2012) describe how HARPS-TERRA obtains them from the HARPS-DRS products. Such chromospheric emission is closely related to the intensity of the stellar magnetic field and larger values imply higher activity levels. The S-index can show periodic variations due to the stellar magnetic cycle (e.g. 11-years magnetic cycle of the Sun), episodic events of higher activity (e.g. flares) and due to the presence of active regions on the rotating surface of the star. The FWHM and the BIS monitor changes in the mean spectral line profiles and should correlate with the presence of spots. Temperature contrast spots and/or magnetic spots are known to show correlation with spurious Doppler signals (e.g. Queloz et al. 2001; Reiners et al. 2013).

Our general strategy is as follows. If a strong periodicity is identified in any of the indices and such periodicity could be related to any of the Doppler signals (compatible period or related through first order aliases), we add a linear correlation term to the model of the Doppler data and perform new samplings of the parameter space. If the data were better described by the correlation term rather than a genuine Doppler signal, the overall model probability would increase and the planet signal in question should disappear (or its amplitude would be altered substantially).

For GJ 163 in particular, none of the activity-indices shows any hints of periodic variability at all. This lack of activity is also supported by the fact that GJ 163 has a mean value of the S-index comparable to the ones measured on the planet prolific M dwarfs GJ 581 and GJ 667C. As a note, we neglected 4 FWHM and BIS measurements due to manifestly incorrect estimates produced by the HARPS-ESO data-reduction software (see Table 3). All 55 measurements of the S-index were successfully extracted by the HARPS-TERRA software and used in the analysis. We show the periodograms of the three activity-indices in Fig. 4.

Fig. 4 Log-likelihood periodograms of HARPS activity indicators of GJ 163: BIS (top panel), FWHM (middle panel), and S-index (bottom panel).

To search for signals in the indices, we used log-likelihood periodograms, which are computationally substantially less intensive and require much less “supervision” than Bayesian Markov chains. These periodograms represent the improvement of the likelihood of the model (log-likelihood periodograms, Baluev 2009) against the null hypothesis for each test period. Compared to more classic periodograms based solely on χ² minimization (e.g. Cumming 2004), our log-likelihood periodogram also adjusts the white noise component of the noise. As any generic periodogram approach, if the model including a periodic signal (sinusoid) substantially improves the merit statistic (likelihood function in this case), a peak over the 1% false alarm probability threshold should emerge. For GJ 163 in particular, none of the indices shows any consistent periodicity at all. For each index, we also attempted to search for a second signal but did not find any significant improvements to the model. In conclusion, and even if some correlations between Doppler variability and the activity indices remain, the variability of activity indices seems purely random and, therefore, a Doppler model with only stochastic noise terms (white and/or red noise) should provide a sufficient description of the data. The absence of signals in the indices supports the interpretation of any periodicity in the Doppler data as a planet candidate.

5. Bayesian analysis of GJ 163 velocities

We analysed the 55 HARPS-TERRA velocities using posterior samplings and by calculating estimates for model probabilities as in e.g. Tuomi (2012), Tuomi & Jones (2012), and Tuomi et al. (2013a). We report the solutions by using maximum a posteriori (MAP) estimates and 99% Bayesian credibility sets (BCSs). In the following subsections, we explore the consequences of analysing the velocities with different noise models, different prior densities, and models containing different numbers of Keplerian signals.

5.1. Reference priors, white noise

We started the analyses of GJ 163 velocities by using the reference priors as defined in Table 1. We used the common white noise model and assumed that any planets orbiting the star do not interact with one another on the time-scale of the observational baseline.

With these assumptions, we estimated the Bayesian evidence for a reference model with k = 0 and searched for periodicities in the data by using a model with k = 1. Our Markov chains converged to a clear signal at a period of 8.63 days. Accounting for this signal increased the model probability to a value 7.3 × 10⁵ times greater than for the model with k = 0.

The two-Keplerian model increased the model probabilities considerably by a factor of 8.0 × 10¹⁰, implying the presence of another periodicity in the data. The second periodicity was found at 225 days and we could constrain this emerging signal according to our detection criteria shown in Tuomi (2012). Samplings of the parameter space of a three-Keplerian model revealed a third significant periodicity at 25.6 days and the corresponding model was found to have a posterior probability of roughly 4.1 × 10⁸ times that of the model with k = 2. However, the MAP periods corresponding to this model (k = 3) were found to be 8.63, 25.6, and 567 days and although we found a local maximum in the period space at 225 days, it does not correspond to a global solution of the three-Keplerian model anymore. We searched for additional periodicities in the data as well, but our Markov chains did not converge to a well-constrained solution for the four-Keplerian model.

To demonstrate the robustness of our solution with k = 3, we performed temperate samplings of the parameter space of the model such that instead of obtaining samples from the posterior density (π) as such, we used π^β as a posterior density, where β is a parameter in the interval [0,1] describing the “temperature” of the sampling. In particular, we chose β = 0.5, and plotted the obtained Markov chain in Fig. 5. In this figure, we show the posterior density corresponding to the temperate sampling as a function of the longest periodicity. It can be seen that the 567-day period corresponds to the global maximum but that there are local maxima as well at periods of roughly 220, 450, and 2000 days. However, none of the local maxima exceeds 0.1% level of the global maximum (solid horizontal line in Fig. 5), which indicates that the 567-day periodicity is clearly the preferred solution of this model.

Fig. 5 Log-posterior density as a function of the longest periodicity in the GJ 163 data based on temperate sampling. The red arrow indicates the highest value in the sample roughly corresponding to the MAP estimate and the horizontal lines show the probability thresholds 10% (dotted line), 1% (dashed line), and 0.1% (solid line) of the maximum.

5.2. Reference priors, white and correlated noise

As pointed out by Baluev (2012), Tuomi et al. (2013b), and Tuomi & Jenkins (2012), RV variations might appear due to another noise component that can be described as having a red colour. This noise component can be accounted for by adding a correlation term in the model of the ith measurement with the previous ones. We applied a moving average (MA) noise model (Tuomi et al. 2013a,b; Tuomi & Jenkins 2012) and observed that even with a k = 0 model, the Bayesian evidences implied that the model including this type of red noise was considerably better than the one with only white noise. We have listed the resulting log-Bayesian evidences (lnP(m)) of pure white noise and moving average models with k = 0,...,3 in Table 4 and plotted the corresponding phase-folded Keplerian signals in Fig. 6.

Fig. 6 Phase-folded signals of the three-Keplerian model with the other three signals subtracted from each panel.

A noteworthy feature in Table 4 is that the MAP solution of the MA model with k = 2 differs from that of the pure white noise model by having a second periodic signal at a period of 25.6 days instead of 225 days. Therefore, it is likely that the apparent periodicity at 225 days can be interpreted as being caused by a combination of yearly aliasing with the 567-day signal (because 1/365 ≈ 1/225−1/567) and/or noise correlations rather than as a genuine Keplerian signal. Indeed, the two-Keplerian model with an MA component has a much higher posterior probability than that with pure white noise. This implies that, once two Keplerian signals are included in the model, the RV noise of subsequent epochs is not independent for k = 2. For k = 3 the situation is reversed and the white noise model is in fact better than the MA model. We believe this is the case because the MA model is slightly overparameterised due to the rather low number of measurements (55) with respect to the number of free parameters (18) and thus penalised by Occam’s razor – the MAP estimate of the correlation parameter φ is clearly positive (0.83) but negligible correlations are also allowed (with a 99% BCS of [0.01, 1]), which implies that the MA component is not necessary for the three-Keplerian model. As was the case with the white noise model, we could not find a significant fourth periodic signal in the data using posterior samplings of a four-Keplerian model and MA noise.

Using the reference prior, we would conclude that there are confidently three Keplerian signals in the data. We have listed the best solution in Table 5. We note that this solution corresponds to an outer planet with an eccentric orbit (e = 0.46), which casts doubt on the stability of the corresponding system in long term. For this reason, we test eccentricity priors that penalise high eccentricities more than the reference prior.

5.3. Alternative priors: σ_e = 0.1 and 0.2, white noise

Using a more restrictive forms for the prior density of eccentricity, i.e. priors with σ_e = 0.2 and σ_e = 0.1 (see Table 1), we repeated the analyses of the GJ 163 velocities with the pure white noise model.

According to our posterior samplings and estimations of Bayesian evidences, when using models with k = 0,...,5 and assuming that all the excess noise in the data is white, the model probabilities increased as k approached four but we could not estimate the Bayesian evidences reliably for k = 5 because we could not spot a fifth periodicity in the data reliably according to our criteria. However, we observed a broad probability maximum in the period space at a period of roughly 1500 days with a low MAP amplitude of 2.7 ms^-1. Specifically, a fifth periodicity could not be constrained from above and below in the period space and we therefore conclude that this prior choice and white noise model favours the existence of four periodic signals in the data. We obtained the same qualitative result for both eccentricity priors.

We have listed the log-Bayesian evidences of models with k = 0,...,4 in Table 6. According to these results, the eccentricity prior with σ_e = 0.2 is much better prior model for the data because it allows reasonably eccentric solutions that are favoured by the data for k = 2,3, and 4. The reason is that the eccentricity of the second strongest signal with a period of 567 days has MAP estimates of approximately 0.45 for models k = 2 and k = 3, respectively. These eccentricities are penalised much more severely with the prior with σ_e = 0.1 which results in decreased Bayesian evidences because the prior actually conflicts with the likelihood function that favours higher eccentricities.

Fig. 7 Phase-folded signals of the four-Keplerian model with the other three signals subtracted from each panel.

For a white noise model with k = 4 the favoured eccentricities of the signal at a period of 567 days and a fourth signal emerging at 125 days are not very high but they still have MAP estimates of roughly 0.2. Yet, eccentricities as high as 0.5 cannot be ruled out either, which is penalised by the eccentricity prior with σ_e = 0.1 and results in a decreased Bayesian evidence for this more restrictive prior (Table 6). Therefore, assuming close-to-circular orbits and that the RV noise is white, we would conclude that the is evidence in favour of four signals in the GJ 163 data. We plotted the phase-folded signals of the four-Keplerian model in Fig. 7 and the remaining RV residuals after subtracting the MAP signals in Fig. 8 to demonstrate that our model also reproduces the data visually.

5.4. Alternative priors: σ_e = 0.1 and 0.2, white and correlated noise

The situation changes considerably when taking into account the possible correlations in the data using the MA model. While there is not much difference in the performance of these two eccentricity priors (Table 6), a fourth periodicity cannot be found according to our detection criteria.

In absolute terms, the model that has the highest posterior probability given the GJ 163 data is the white noise model with σ_e = 0.2. While the MA model is much better for k = 0,1, and 2 with any of the eccentricity priors, the pure white noise model is favoured by the data for k = 3 and is even better for k = 4 (Table 6). The likely reason for this result is that the MA parameter φ is actually consistent with zero for k = 3 (Fig. 9), despite that its MAP estimate is 0.83, which makes the model overparameterised and decreases its Bayesian evidence estimate.

We report the best solution corresponding to the white noise model and four Keplerian signals in Table 7 and Fig. 10. The corresponding probabilities of this model with the moderate eccentricity prior with σ_e = 0.2 are 4.8 × 10^-28, 3.2 × 10^-22, 7.1 × 10^-12, 0.026, and 0.974 for k = 0,1,2,3, and 4, respectively. We note that the four-Keplerian model does not have a posterior probability that exceeds the detection threshold (by being at least 150 times greater than that of the three-Keplerian model). While the probability of the model with k = 4 is only ~38 times greater than that of the model with k = 3, the other two detection conditions are satisfied – all periods are well-constrained and the amplitudes are statistically different from zero, as can be seen in Fig. 10. Therefore, we present a tentative four-Keplerian solution as the preferred one.

We note that using a proper informative prior for the jitter parameter (σ_J) did not change the results significantly for any of the statistical models. The jitter prior with μ_σ = 1 ms^-1 and σ_σ = 1 ms^-1 penalised models with k = 0 and k = 1 slightly but does not have a significant effect on neither the Bayesian evidences nor the obtained posterior densities. We believe this happens because the jitter prior does not represent the data very well for k = 0 and 1 because the two corresponding models contain at least two or three signals that increase the variability of the data the models interpret as noise. This can be seen as additional support for the existence of more than one signals in the data. However, conclusions can only be based on the Bayesian evidences and the corresponding model probabilities, unless there was a very strong a priori reason to believe that the excess jitter cannot exceed 1–2 ms^-1 level considerably, which is not the case here.

5.5. Period prior comparison

Finally, we tested if the results were affected significantly by the chosen period prior, i.e. whether the prior was constructed by assuming a uniform density in the period space, orbital frequency space, or log-period space. We tested this using the three-Keplerian model and an MA noise model and calculated the log-Bayesian evidences using each period prior. We obtained log-Bayesian evidences that were within 0.5 from one another, which cannot be considered a significant difference in any model selection problem because it corresponds to probabilities between 0.23–0.45 for the three models. The likely reason is that in the vicinity of the period maxima, all these period priors are roughly constant and constant coefficients do not have an effect on the estimates of Bayesian evidences. Also, according to our tests with different values of k, the different period priors did not affect the ability to detect signals in the data in the first place compared to the reference prior that was uniform in log-period.

Fig. 8 Residuals of the four-Keplerian orbital solution shown in Fig. 7.

Fig. 9 Distribution of the correlation parameter for an MA noise model and k = 3 with a hyperparameter σe = 0.2.

Fig. 10 Marginal distributions of the orbital periods (Px), eccentricities (ex), and RV amplitudes (Kx) of the four-Keplerian solution in Table 7. The solid curves are Gaussian densities with the same mean and variance as the marginal distributions.

6. Dynamical feasibility of the system

As a final validation of the signals as planet candidates, we performed a dynamical analysis on representative samples of parameters drawn from the posterior densities. According to the integrations we performed by using the Bulirsch-Stoer algorithm (Bulirsch & Stoer 1966), all parameter vectors drawn from the posterior density of the model with k = 3 (samples compatible with Table 7) corresponded to stable planetary systems on a time-scale of 10⁶ years or longer. We repeated the same experiment with samples from the posterior density of the model with k = 4 (represented by Table 7). Again, all vectors corresponded to stable orbital configurations beyond 10⁶ years and, therefore, we could not use any dynamical argument to decide which solution was the most favored one from a physical point of view.

According to our numerical integrations of the orbits, the orital elements remained in bounded areas of the parameter space without corresponding to collisions, orbital crossings, or escapes from the system (beyond 10 AU), in accordance with the Lagrange stability criteria. Following Tuomi (2012) and Tuomi et al. (2013b), we demonstrate this by plotting the approximated Lagrange stability boundaries of the solutions in Tables 5 and 7 in Fig. 11. This figure illustrates that our three- and four-planet solutions correspond to planetary systems with sufficient orbital spacing to enable long-term stability.

Fig. 11 Approximated Lagrange stability boundaries indicating the parameter space around the MAP estimates (red circles) where there are no stable orbits (shaded areas). Three- (top) and four-planet (bottom) solutions.

Finally, we tested the chaotic behaviour of our solutions by using the frequency analysis method of Laskar (1993) as applied in e.g. Correia et al. (2010) and Tuomi et al. (2013a). According to our results, the relative variation in the mean motions of the planets (D-index) is at most 10^-4 for the most eccentric solutions of the three-Keplerian model, indicating that such orbits might show some chaotic behaviour, whereas typically we found this variation to be roughly 10^-6 or even less, which indicates very regular and thus stable motion.

7. Discussion and conclusions

Based on our analyses of the GJ 163 velocity data from the HARPS spectra, the absence of any clear periodicity in the activity indices, and the dynamical feasibility of the orbital solutions, we conclude that GJ 163 has at least three planet candidates orbiting it. This conclusion is independent of noise models and prior densities and is therefore very strongly supported by the data. When using a white noise model and slightly more limiting eccentricity priors, we obtained a solution for the four-Keplerian model that was well constrained in the parameter space (Table 7 and Fig. 10) but did not exceed the detection threshold of being 150 times more probable than the three-Keplerian model. We interpret this result as suggestive but inconclusive evidence for a fourth planet candidate in the system. Given that two of the three detection conditions of Tuomi (2012) are satisfied for this fourth candidate, even a small amount of additional measurements (~20 new observations) might settle the issue in favour of or against the fourth candidate. Additional measurements are also strongly encouraged to investigate whether a very broad local probability maximum we observed in the period space at roughly 1500 days corresponds to yet another candidate.

Having presented results based on different models and priors, we discuss briefly how one should proceed in similar analyses of RV data from HARPS and other instruments. First, it is crucial to investigate whether autocorrelation in a suitable timescale can explain some of the variation in the data in comparison to models where this variation is described using a Keplerian model. If autocorrelation – in a timescale of few dozen days (Tuomi et al. 2013b; Tuomi & Jenkins 2012; Baluev 2012) – is a good description of the data and the inclusion of Keplerian signals instead of autocorrelation does not improve the model much, one cannot safely claim that there is evidence for additional planetary signals in the data. However, if the signals satisfy the detection criteria of Tuomi (2012) even with an MA model, or a similar model that accounts for the correlated noise, it is likely that the signals are real and possibly of planetary origin if they do not have activity-related counterparts. It is always possible that any given statistical model is not an adequate description of the velocity variations (e.g. Tuomi et al. 2011) and, as a consequence, some signals may still be spurious artefacts of poor modelling. Detecting low-mass planets around nearby stars is not only a matter of finding significant signals – we can only start calling them planetary candidates if these signals are 1) reasonably independent of the exact choice of a noise model (given some good candidate noise models); 2) do not have counterparts in the stellar activity data; and 3) correspond to planetary systems that are physically viable.

The same applies to the choice of priors. We have studied the consequences of assuming different informative priors for the orbital eccentricity. According to our results, the exact choice of this prior does not affect the results much. The only exception in case of GJ 163 velocities is the evidence in favour of a fourth companion when orbital eccentricities are kept a priori close to zero.

We also investigated the effect of uninformative prior densities of the period parameter that correspond to uniform priors on different but still plausible parameterisations. For GJ 163 data, the precise choice of the prior did not affect the results in a significant way.

We conclude that the detections of planets in RV data does not seem to be very sensitive to prior choices in practice as long as the chosen priors are not too informative w.r.t. the likelihood functions, i.e. delta-function priors or other very narrow probability densities. As a general rule, we find it is always advisable to repeat the statistical analyses by assuming a few different prior choices to investigate under which hypotheses solutions corresponding to low-amplitude planets are supported by the data (e.g. does one really need to assume close-to-circular orbits to obtain physically viable solutions?).

According to our results (see also Tuomi & Jenkins 2012), it appears that accounting for correlations on the time-scale of ~10 days might disable the detectability of some signals in RV data. When the significance of a signal depends on whether one uses a correlated noise term in the statistical model, it might not be possible to tell if the signal is a genuine one unless the chosen noise model that enables the detection produces a much better solution than the alternative noise models. Regarding GJ 163, the detection of a possible companion at 125 days is not possible if an MA component is included in the noise model (the corresponding signal is completely absorbed by the correlation term and the Markov chains do not converge to a fourth periodicity). On the other hand, a simpler white noise model enables a convergence to a clear four-Keplerian solution and provides a substantially higher model probability. For this data set, it looks like the MA model is over-parameterised and that the inclusion of a fourth periodic signal would be our preferred model choice. Therefore, the best noise model depends on the number of signals and the preferred number of signals depends on the selected noise model. When this happens, the only way to tell what is the most probable situation in reality is to examine all realistically possible combinations (white noise only, red and white noise, k = 0,1,...) and derive the most informed result from there.

The liquid-water habitable zone (HZ) of GJ 163 is located between roughly 0.123 and 0.275 AU (assuming 0% clouds in the equations of Selsis et al. 2007). This suggests that GJ 163 c, with a semi-major axis of 0.126 AU and a likely circular orbit, is located inside the stellar HZ throughout most of its orbital cycle. The minimum mass of GJ 163 c is m_psini = 8.7 [5.3, 12.3] M_⊕ and its expected true value (~13.0  M_⊕) lies therefore roughly in the Neptune-mass regime. If the planet were found to transit in front of the star, the minimum mass would be its true mass, and its size could be obtained and nature elucidated (scaled-down version of Neptune rather or a massive rocky planet with a solid surface). It is therefore a primary target for transit follow-up observations. Given the uncertainties and detailed modelling necessary to account for all the unknowns, assessing the habitability of GJ 163 c in detail is beyond the scope of this work.

The detection of three (possibly four) planet candidates adds GJ 163 to the emerging population of diverse planetary systems around low-mass stars. The architecture of the GJ 163 system resembles a scaled-down version of the Solar System in the sense that the system consists of two super-Earth mass objects with minimum masses of 9.9 and 8.7 M_⊕ and short orbital periods of 8.6 and 25.6 days, respectively; a possible cool Neptune with a minimum mass of 14.0 M_⊕ and an orbital period of 125 days; and a more massive sub-Saturnian outer planet with an orbital period of 572 days.

Online material

Acknowledgments

M. Tuomi acknowledges D. Pinfield and 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 Federal Ministry of Education and Research under 05A11MG3. Based on data obtained from the ESO Science Archive Facility under request number ANGLADA36104. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. The authors acknowledge the significant efforts of the HARPS-ESO team in improving the instrument and its data reduction pipelines that made this work possible. We also acknowledge the efforts of all the individuals that have been involved in observing the target star with HARPS.

References

All Tables

All Figures

	Fig. 1 Eccentricity distribution (red circles) of low-mass (mp < 0.1  MJup) exoplanet candidates and the relative eccentricity priors with σe = 0.3 (solid line), σe = 0.2 (dash-dotted line), and σe = 0.1 (dashed line). The priors have been scaled to have maxima of 0.37 to enable visual comparison with the distribution.
In the text

Fig. 2 CaII H line of GJ 163 (purple) as obtained from co-adding the 55 available HARPS spectra. The same lines of stars with similar spectral type are shown for comparison (offset in wavelength). The height of the line was scaled to the mean number of counts per pixel in this particular HARPS echelle order. Spectral types were taken from Bonfils et al. (2013). As discussed in the text, GJ 163 is among the least chromospherically active stars in its spectral range.

In the text

	Fig. 3 HARPS-TERRA radial velocities of GJ 163 with the data median removed.
In the text

	Fig. 4 Log-likelihood periodograms of HARPS activity indicators of GJ 163: BIS (top panel), FWHM (middle panel), and S-index (bottom panel).
In the text

	Fig. 5 Log-posterior density as a function of the longest periodicity in the GJ 163 data based on temperate sampling. The red arrow indicates the highest value in the sample roughly corresponding to the MAP estimate and the horizontal lines show the probability thresholds 10% (dotted line), 1% (dashed line), and 0.1% (solid line) of the maximum.
In the text

	Fig. 6 Phase-folded signals of the three-Keplerian model with the other three signals subtracted from each panel.
In the text

	Fig. 7 Phase-folded signals of the four-Keplerian model with the other three signals subtracted from each panel.
In the text

	Fig. 8 Residuals of the four-Keplerian orbital solution shown in Fig. 7.
In the text

	Fig. 9 Distribution of the correlation parameter for an MA noise model and k = 3 with a hyperparameter σe = 0.2.
In the text

	Fig. 10 Marginal distributions of the orbital periods (Px), eccentricities (ex), and RV amplitudes (Kx) of the four-Keplerian solution in Table 7. The solid curves are Gaussian densities with the same mean and variance as the marginal distributions.
In the text

	Fig. 11 Approximated Lagrange stability boundaries indicating the parameter space around the MAP estimates (red circles) where there are no stable orbits (shaded areas). Three- (top) and four-planet (bottom) solutions.
In the text