EDP Sciences
Free Access
Issue
A&A
Volume 599, March 2017
Article Number A126
Number of page(s) 7
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201630050
Published online 13 March 2017

© ESO, 2017

1. Introduction

The discovery of a potentially terrestrial planet in an approximately 11-day orbit around the M dwarf Proxima Centauri, announced by Anglada-Escudé et al. (2016; hereafter AE16), represents a major breakthrough in the extrasolar-planet science. Proxima is the closest known star to the Sun, and the planet Proxima b might have an equilibrium surface temperature suitable for liquid water to exist. Intrinsic stellar activity can mask planetary signals measured in radial velocity time series, especially for M dwarf systems. This stellar activity signal can be of the order of a few  m s-1, larger than the semi-amplitude of the planetary signal previously measured for the Proxima system. To date, several techniques have been developed to mitigate this effect. Dumusque et al. (2017) devised a blind test, based on simulated RV measurements including realistic stellar activity noise (Dumusque 2016), to compare the performances of different methods. To derive orbital and physical properties of Proxima b from high-precision RV measurements it is therefore necessary a proper treatment of the stellar activity (the “noise”) contribution. To do so, AE16 used a Bayesian framework consisting of a moving average term and linear correlations with activity indexes. The technique used by AE16 was tested on the RV datasets of Dumusque et al. (2017) and found to be very effective for recovering signals induced by low-mass planets. A different Bayesian technique applied to the test which performed similarly well is based on modelling the correlated noise using Gaussian processes (GPs) (see e.g. Rasmussen & Williams 2006; and Roberts et al. 2012, for a general description of GPs). GPs are a powerful tool to mitigate the correlated noise in a set of measurements, such as the stellar activity signature in RV data. This technique has been used successfully in several recent works dealing with real RV data (Haywood et al. 2014; Grunblatt et al. 2015; Rajpaul et al. 2015; Affer et al. 2016; Faria et al. 2016; Lopez-Morales et al. 2016).

In this article we have applied the GP technique to the radial velocity observations of Proxima published in AE16, taken with the UVES and HARPS spectrographs. HARPS data are divided in a pre-2016 and 2016 (PRD: Pale Red Dot campaign) subsamples. We aim to confirm and characterise the signal of Proxima b in an independent way through an alternative and robust modelling of the stellar noise. We also consider the existence of a second planetary signal, by testing the hypothesis in a Bayesian framework.

2. Description of the model and technique

We model the correlated noise by adopting a quasi-periodic covariance function, similarly to what was done, for example, in Affer et al. (2016). This function is described by parameters, called hyperparameters, which try to model in a simple fashion some of the physical phenomena undelying the stellar noise. The quasi-periodic kernel is described by the covariance matrix (1)where t and t indicate two different epochs. This kernel is composed by a periodic term coupled to an exponential decay one, in order to model a recurrent signal linked to stellar rotation and taking into account the size-evolution of finite-lifetime active regions. Such approach is therefore particularly suitable to model the stellar noise on short-term timescales, as we consider that modulated by the stellar rotation period. The hyperparameter h of the covariance function represents the amplitude of the correlations; θ is related to the rotation period of the star; w is the length scale of the periodic component, linked to the size evolution of the active regions; λ is the correlation decay timescale, that we assume to be related to the active regions lifetime. Here, σinstr,RV(t) is the RV internal error at time t for a given instrument; σinstr,jit is the additional uncorrelated “jitter” term, one for each instrument, that we add in quadrature to the internal errors to take into account instrumental effects and other 1 m s-1 noise sources included neither in σinstr,RV(t) nor in our stellar activity framework; δt,t is the Dirac delta function.

In the GP framework, the log-likelihood function to be maximised by the Markov chain Monte Carlo (MCMC) procedure is (2)where n is the number of the data points, K is the covariance matrix built from the covariance function in Eq. (1), and represents the RV residuals, obtained by subtracting the Keplerian signal(s) from the original RV dataset.

The general form for the models that we tested in this work is given by the equation (3)where nplanet = 1,2; ν is function of time t, time of the inferior conjuntion Tc,j, orbital period Pj, eccenticity e and argument of periapsis ωj; γinstr is the RV offset, one for each instrument; is the secular acceleration; GP is the stellar noise modelled with the Gaussian process. Instead of fitting separately ej and ωj, we use the auxiliary parameters and to reduce the covariance between ej and ωj.

Our analysis is performed in two steps through two separate MCMC runs. First, we investigate the existence of the signal attributed to Proxima b by including one Keplerian in the model together with the GP noise term. Through our GP analysis, we can constrain the values of our model hyperparameters, and thus constrain the stellar properties we assume these hyperparameters are related to. Here we examine the existence of a second planetary signal in the RV data, whose presence could not be ruled out at long orbital periods, as suggested by AE16. Finally, we evaluate the Bayesian evidences (marginal likelihoods) for the two models to assess whether the present RV dataset supports the two-planet scenario over the other. We use the publicly available emcee algorithm (Foreman-Mackey et al. 2013) to perform the MCMC analyses of the RV data, and the publicly available GEORGE Python library to perform the GP fitting within the MCMC framework (Ambikasaran et al. 2014). We used 150 random walkers for each MCMC run, and to derive the parameter posterior distributions we performed a burn-in by removing from the samples those with lnℒ lower than the median of the whole lnℒ dataset. The strategy of our MCMC analysis is to put as little constraint as possible to the model parameters. Following such a philosophy, the starting points of the walkers in the parameter space were randomly chosen from normal distributions with σ sufficiently large to allow a broad exploration of the parameter space around guess values. We also used uninformative priors for almost all the jump parameters. The MCMC convergence is checked using the diagnostic proposed by Ford (2006), and the best fit values and uncertainties for each jump parameter are calculated as the median of the marginal posterior distributions and the 16% and 84% quantiles. We quantify the relative statistical significance of the two tested models by deriving their Bayesian evidences and , according to Bayes’ theorem. The selection of the model better describing the data is made by calculating the ratio /. In fact, by assuming that there is not any a priori reason to favour one model over the other, the -ratio represents the Bayes factor, that is, the figure of merit encoding all the support the data give to one model over the other (see, e.g. Díaz et al. 2016, for a detailed discussion in the framework of the extrasolar-planets detection). We estimate the natural logarithm of each evidence by using two analytic approximations due to Chib & Jeliazkov (2001) and Perrakis et al. (2014; hereafter indicated as the CJ01 and Perr14 estimators). For the statistical interpretation we follow here a conventional empirical scale, the so-called Jeffreys’ scale, which states that the model with the highest is strongly favoured over the other when , while denotes moderate evidence, weak evidence, and corresponds to inconclusive evidence.

Table 1

Prior probability distributions for the one-planet model parameters.

thumbnail Fig. 1

Distributions of total MCMC values for the jump parameters of the one-planet model as function of the natural logarithm of the likelihood function.

Open with DEXTER

thumbnail Fig. 2

Distributions as in Fig. 1 after constraining the GP hyperparameter θ in the range [85, 89] days. The vertical red lines represent the median values (solid) and the 16th and 84th percentiles (dotted) of these distributions, which are listed in Table 2.

Open with DEXTER

3. One-planet model

In order to minimally constrain our model, we have chosen Bayesian prior probability distributions as described in Table 1. All are uniform (uninformative) priors, except for the orbital eccentricity: for this we adopted a Beta distribution, because it is well-known that the noise tends otherwise to favour higher eccentricities (e.g. see Appendix B.4 in Gregory 2010, April 2016 supplement1). We considered the GP periodic hyperparameter θ and the orbital period P as scale invariant parameters, for which an uninformative prior is one that is uniform in lnP. We have chosen a very large range for the GP hyperparameter w, which usually is constrained to values less than one when there is clear evidence of a stellar-rotation-linked periodic component in the correlated noise (e.g. Lopez-Morales et al. 2016). Our MCMC run stopped at ~200 000 steps: the median of lnℒ calculated over the 150 chains maintained a nearly constant value over ~140 000 steps. After the burn-in of the first 100 000 steps, we have performed an additional burn-in in the lnℒ space by deriving the marginal posterior distributions of samples with lnℒ > median(lnℒ) (see Fig. 1). The posterior distribution of the GP hyperparameter θ shows a concentrated and symmetrical sub-sample between 85 and 89 days with the highest values of lnℒ, which is the range where the known rotation period of Proxima falls in. Considering that θ is expected to represent the stellar rotation period, this evidence is particularly relevant because it reveals that the correlated noise does contain a periodic term which is modulated by the rotation of the star. The rotation period of Proxima has been robustly estimated from ASAS V-band photometry over an interval spanning up to 15 yr (e.g. see Fig. 3 in Robertson et al. 2016, or Fig. 3 in Wargelin et al. 2017). In particular, Wargelin et al. (2017) present the Lomb-Scargle periodogram of the ASAS-3 and ASAS-4 photometry covering 15 yr, which shows a very significant and sharp peak at 83.1 days. The same authors show evidence for differential rotation in the light curve by analysing periodograms of different observing seasons, finding significant peaks ranging from Prot = 77 to Prot = 90 days. Based on this robust result from photometry, we can naturally restrict our analysis to the 697 510 posterior samples for which θ falls in the range 85–89 days, ensuring the periodic features of our model are tied to the stellar rotation. The corresponding posterior distributions for all the model parameters are showed in Fig. 2. For the GP hyperparameter w we get a symmetric distribution with all the values less than one, which is what is expected, as described, for example, in Lopez-Morales et al. (2016). Moreover, the distribution of the GP hyperparameter λ has a median of 311 days and is nearly symmetric (excluding the points with λ< 200 days that are relatively few compared to the whole sample, as indicated by the value of the 16th percentile). This long evolutionary timescale appears particularly interesting when compared with photometric results. Wargelin et al. (2017) show the evidence for the rotational phasing to remain remarkably constant over the ~15 yr timespan at a fixed Prot = 83.1 days, and for the modulation amplitude not to change significantly over two years, which could be indicative of the presence of persistent active longitudes on the stellar photosphere. For instance, the light curve corresponding to epochs between 2010 and 2012, which overlap with the timespan of the RV data, can be well described by a sinusoidal function with constant semi-amplitude and phase (see Fig. 3 in Wargelin et al. 2017), suggesting that the photospheric structures responsible for the correlations observed in the RV residuals could indeed have a decay timescale of the order of approximately one year. The best-fit estimates for all the model parameters are showed in Table 2, along with values for the planet orbital solution reported by AE16. The residuals of our global model have an rms of 1.3 . In Fig. 3 we show the RV residuals, after removing the stellar signal described by the GP, folded at our best-fit estimate of the orbital period (Fig. 4 shows the same plot for each dataset taken separately). It can be seen that a circular model well describes the data, in agreement with our finding that the eccentricity is consistent with zero within ~1.4σ, whilst usually at least a 2.5σ level is required for claiming a significant eccentricity. Figure 5 shows the RV residuals, after subtracting the Keplerian solution, with our best-fit model for the correlated stellar noise superposed. We note that the model describes the second half of the HARPS-PRD dataset particularly well.

Table 2

Percentiles of the distributions in Fig. 2, referred to the cut for θ in the range 85–89 days, compared with values found by AE16.

thumbnail Fig. 3

Radial velocity residuals, after subtracting the best-fit GP correlated noise model, folded at the best-fit orbital period P = 11.1855 days. Different symbols and colours are used for each dataset. Grey dots represent average values for 20 phase bins. Superposed are the best-fit eccentric model (continuous red curve) and the simpler circular model (dashed red curve).

Open with DEXTER

AE16 show (extended data Fig. 2 therein) that their likelihood-ratio periodogram of the HARPS pre-2016 RV dataset has the highest peak at P = 215 days (p-value = 1%), but they could not conclude about the nature of this signal. To investigate in more detail the possible cause of its existence, we analyse the same dataset with the generalised Lomb-Scargle algorithm (GLS; Zechmeister & Kürster 2009). We find that the highest power density peak appears at P ~ 57 days (upper plot in Fig. 6), while at P = 215 days there is no significant power (p-value > 1%, estimated through a bootstrap with re-sampling analysis on 10 000 fake datasets). Then, we run GLS on the RV residuals of our global one-planet plus noise model. We do not find power excess in the periodogram at low frequencies (lower plot in Fig. 6), indicating that our correlated noise model has suppressed any long-term modulation in the data. The impact of the GP model on these residuals is not trivial to be characterised, and could prevent the detection of additional planetary signals in the periodogram. Therefore, we explore an additional two-planet plus stellar activity model to test whether this signal is more effectively modelled as a planet or stellar activity.

thumbnail Fig. 4

Radial velocity residuals for the one-planet model, after subtracting the best-fit GP correlated noise model, folded at the best-fit orbital period P = 11.1855 days. Each dataset is showed here separately for clarity. Superposed are the best-fit eccentric model (continuous red curve) and the simpler circular model (dashed red curve).

Open with DEXTER

thumbnail Fig. 5

Radial velocity residuals time series (black dots), after subtracting our best-fit orbital solution for Proxima b. The blue line with grey shaded 1σ regions represents our best-fit GP quasi-periodic model for the correlated stellar noise. The upper plot shows the complete dataset, while the two plots in the second row show selected epochs, to easier visualize the agreement between the model and the data.

Open with DEXTER

thumbnail Fig. 6

Upper plot: GLS periodogram of the HARPS pre-2016 RV data. The vertical red dashed line marks P = 215 days, the period found by AE16 in their likelihood-ratio periodogram, for which they could not give a physical interpretation. The p-value levels have been derived through a bootstrap analysis. Lower plot: GLS periodogram of the complete RV datates residuals after subtracting the GP plus one-planet model. The inset plot shows the same periodogram around P = 215 days, for more clarity.

Open with DEXTER

4. Two-planet model

We tested the two-planet scenario using Gaussian priors for all the GP hyperparameters, except for h, based on the physically reliable results obtained for the stellar noise representation for the one-planet model. In modelling the RVs we ignored the mutual gravitational perturbations between the two planets. Taking further advantage of the one-planet model results to speed up the analysis, for Proxima b we fixed the upper limit for RV semi-amplitude to 3  m s-1, that of the orbital period to 20 days and the eccentricity to zero. We explored the parameter space for the orbital period of the second Keplerian starting from 100 days, on the basis that no clear modulation with a period equal or less than that is observed in the HARPS-PRD RV dataset (see AE16). We ran both a model where the orbit of the outer planet is assumed circular, and one where the eccentricity is treated as a free parameter. The list of the adopted priors is showed in Table 3. The marginal posterior distributions for the simpler circular model are showed in Fig. 7. For this case, we note that the orbital period of the outer planet appears to be multi-modal, and the semi-amplitude of the Keplerian is below the average value of the internal RV uncertainties (σRV⟩ = 0.94 m s-1). Moreover, the rms of the residuals remains the same between the one and two-planet model (1.3 and 1.2 m s-1, respectively). These evidences alone do not favour the existence of a second Keplerian signal in the data. A similar conclusion can be drawn for the model with the eccentricity as free parameter. Here, the eccentricity appears to be insignificant within 2σ.

Table 3

Prior probability distributions for the two-planet model parameters.

Table 4

Two-planet model parameter estimates, obtained by fixing the eccentricity of the second Keplerian to zero, and also keeping it as a free parameter.

5. Summary and conclusions

We used a Gaussian process framework to allay the stellar correlated noise in the radial velocity timeseries of Proxima Centauri published by Anglada-Escudé et al. (2016). In our fitting procedure we considered only radial velocity measurements, without including other ancillary data. We then compared the outputs of our procedure with those coming from photometry, to provide a consistent physical interpretation of the results. We adopted a quasi-periodic kernel to describe the correlated noise: this is a widely used function that represents the covariance between measurements at different epochs in terms of parameters that can be related to some physical properties of the star. The philosophy behind our approach – but not the analysis tools – is similar to that of Faria et al. (2016), which used GPs to model the stellar noise in RVs of the host-star CoRoT-7 and were able to retrieve planets CoRoT-7 b and CoRoT-7 c without analysing any other additional dataset.

In our study, we confirm the detection of a coherent signal well described by a Keplerian orbit equation that can be attributed to the planet Proxima b. We find system parameters to be in good agreement with those of Anglada-Escudé et al. (2016). We note that our estimates of the uncorrelated jitter terms for each independent RV dataset are lower than in Anglada-Escudé et al. (2016, Table 2), suggesting that our model removes the stellar noise more efficiently.

As a major outcome, the correlated noise has indeed a periodic component modulated by the stellar rotation period, and we can describe it through a simple GP-based regression. The best-fit estimates of all the hyperparameters assume realistic physical values, as the stellar rotation period or the typical lifetime of active regions, because we have a measured or deduced counterpart from an independent photometric dataset. Therefore our model is physically trustworthy as well as simpler, because it makes use of fewer free parameters. This is likely the reason because we find smaller errors for the Doppler semi-amplitude and the orbital period. We note that Anglada-Escudé et al. (2016) do not mention the existence in their RV periodograms of a significant signal with a frequency close to the stellar rotation period, but in light of our results the rotation of Proxima could be considered the main cause of the clear trend visible in the high-sampling HARPS-PRD dataset.

Table 5

Bayesian evidences for the tested models.

thumbnail Fig. 7

Marginal posterior distributions for the two-planet model (circular option).

Open with DEXTER

Once constrained the stellar noise on the results found for the one-planet model, we tested the presence of a second companion by adding a Keplerian to the model. Table 5 shows the Bayesian evidences for the considered models from two different estimators. Both strongly favour the one-planet scenario, with ) between 10 and 18. Our analysis thus shows that, with our proposed treatment of the stellar noise, the present RV dataset does not show unambiguous evidence for the existence of an additional, low-mass planetary companion to Proxima on an outer orbit with an orbital period of between 100 to 6000 days.


Acknowledgments

This work was conceived and partially carried out during the “Approaching the Stellar Astrophysical Limits of Exoplanet Detection: Getting to 10 cm s-1” workshop held at the Aspen Center for Physics on Sept. 2016. The Aspen Center for Physics is supported by National Science Foundation grant PHY-1066293. M.D. was granted by the Simons Foundation to attend the workshop, and he acknowledges funding from Progetto Premiale INAF “Way to Other Worlds” (WoW), decreto 35/2014. F.D.S. acknowledges the Swedish Research Council International Postdoc fellowship for support. We thank A. Sozzetti, A.S. Bonomo and R. Haywood for valuable discussions, and the referee for useful comments. We finally thank M. Boldi and C. De Sica for some pleasant hours spent in Aspen, twenty years in the past.

References

All Tables

Table 1

Prior probability distributions for the one-planet model parameters.

Table 2

Percentiles of the distributions in Fig. 2, referred to the cut for θ in the range 85–89 days, compared with values found by AE16.

Table 3

Prior probability distributions for the two-planet model parameters.

Table 4

Two-planet model parameter estimates, obtained by fixing the eccentricity of the second Keplerian to zero, and also keeping it as a free parameter.

Table 5

Bayesian evidences for the tested models.

All Figures

thumbnail Fig. 1

Distributions of total MCMC values for the jump parameters of the one-planet model as function of the natural logarithm of the likelihood function.

Open with DEXTER
In the text
thumbnail Fig. 2

Distributions as in Fig. 1 after constraining the GP hyperparameter θ in the range [85, 89] days. The vertical red lines represent the median values (solid) and the 16th and 84th percentiles (dotted) of these distributions, which are listed in Table 2.

Open with DEXTER
In the text
thumbnail Fig. 3

Radial velocity residuals, after subtracting the best-fit GP correlated noise model, folded at the best-fit orbital period P = 11.1855 days. Different symbols and colours are used for each dataset. Grey dots represent average values for 20 phase bins. Superposed are the best-fit eccentric model (continuous red curve) and the simpler circular model (dashed red curve).

Open with DEXTER
In the text
thumbnail Fig. 4

Radial velocity residuals for the one-planet model, after subtracting the best-fit GP correlated noise model, folded at the best-fit orbital period P = 11.1855 days. Each dataset is showed here separately for clarity. Superposed are the best-fit eccentric model (continuous red curve) and the simpler circular model (dashed red curve).

Open with DEXTER
In the text
thumbnail Fig. 5

Radial velocity residuals time series (black dots), after subtracting our best-fit orbital solution for Proxima b. The blue line with grey shaded 1σ regions represents our best-fit GP quasi-periodic model for the correlated stellar noise. The upper plot shows the complete dataset, while the two plots in the second row show selected epochs, to easier visualize the agreement between the model and the data.

Open with DEXTER
In the text
thumbnail Fig. 6

Upper plot: GLS periodogram of the HARPS pre-2016 RV data. The vertical red dashed line marks P = 215 days, the period found by AE16 in their likelihood-ratio periodogram, for which they could not give a physical interpretation. The p-value levels have been derived through a bootstrap analysis. Lower plot: GLS periodogram of the complete RV datates residuals after subtracting the GP plus one-planet model. The inset plot shows the same periodogram around P = 215 days, for more clarity.

Open with DEXTER
In the text
thumbnail Fig. 7

Marginal posterior distributions for the two-planet model (circular option).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.