Our nearest neighbor, Proxima Centauri, hosts a temperate terrestrial planet. We detected in radial velocities evidence of a possible second planet with minimum mass m c sin i c = 5.8 ± 1.9M ⊕ and orbital period P c = 5.21 − 0.22 + 0.26 years. The analysis of photometric data and spectro-scopic activity diagnostics does not explain the signal in terms of a stellar activity cycle, but follow-up is required in the coming years for confirming its planetary origin. We show that the existence of the planet can be ascertained, and its true mass can be determined with high accuracy, by combining Gaia astrometry and radial velocities. Proxima c could become a prime target for follow-up and characterization with next-generation direct imaging instrumentation due to the large maximum angular separation of ~1 arc second from the parent star. The candidate planet represents a challenge for the models of super-Earth formation and evolution.

In this work, we present the results from the analysis of the extended RV dataset carried out within the framework outlined in ( 11 ) and aimed at searching for additional low-mass planetary companions to Proxima. The conclusions of this study are supported by the analysis of spectroscopic activity diagnostics and of a photometric light curve with a baseline longer than the time span of the RV dataset.

Over more than 15 years, our nearest stellar neighbor Proxima Centauri (GJ 551; hereafter Proxima) has been observed with different techniques aimed at the detection of planetary companions. Proxima is an M5.5V star 1.3008 ± 0.0006 pc away from the Sun ( 1 ); therefore, it is an ideal target for astrometric ( 2 ) and direct imaging ( 3 ) searches. To date, these methods excluded the existence of Jupiter mass planets from 0.8 astronomical unit (AU) to farther than 5 AU (>2 AU for masses m ≥ 4M Jupiter ) ( 2 , 3 , 4 ). Holman and Wiegert ( 5 ) predicted a maximum stable orbital radius of 1700 AU for planets orbiting Proxima, because the star orbits the double system αCen AB, as has been demonstrated with a high degree of confidence in ( 6 , 7 ). More recently, Kervella et al. ( 7 ) set a 1σ upper limit of 0.3 M Jupiter to potential companions of Proxima up to 10 AU by analyzing its proper motion taken from the Gaia Data Release 2 (DR2) and excluded the presence of planets between 10 and 50 AU in the mass range 0.3 to 8 M Jupiter . Using the radial velocity (RV) technique, Endl and Kürster ( 8 ) excluded the presence of planets with minimum masses greater than 1 M Neptune out to 1 AU, while they could have detected planets with minimum masses greater than 8.5 M ⊕ and orbital periods out to 100 days. It was because of the RV technique that the temperate, low-mass planet Proxima b, orbiting at a distance of ∼0.05 AU, was discovered ( 9 ). M dwarfs have high occurrence rates of small planets (1.0 to 2.8 R ⊕ ), 3.5 times more than main-sequence FGK stars ( 10 ); therefore, systems with multiple, small, low-mass planets are expected to be common around them. With their RV dataset, including measurements with the HARPS (High Accuracy Radial Velocity Planet Searcher) and UVES (Ultraviolet and Visual Echelle Spectrograph) spectrographs, Anglada-Escudé et al. ( 9 ) could not rule out the presence of an additional super-Earth in the system with orbital periods longer than that of Proxima b and with Doppler semi-amplitude smaller than 3 m s −1 . For the sake of an independent confirmation of the existence of Proxima b and to search for additional low-mass companions, Damasso and Del Sordo ( 11 ) analyzed the same RV measurements using a model that treats the imprint of the stellar activity in a different way than the method adopted in ( 9 ). This analysis uncovered correlated variability in the RV data ascribable to the stellar activity and modulated over the known stellar rotation period. By treating the stellar activity signal with a quasi-periodic model, they could not unambiguously detect a low-mass companion with an orbital period longer than that of Proxima b. The search for additional planets in this system did not stop, and it was at the origin of the Red Dots (RD) initiative ( https://reddots.space/ ), which also focused on other nearby stars, and recently led to the discovery of a candidate super-Earth orbiting Barnard’s star close to the snowline ( 12 ). Because of the RD campaign, additional RVs of Proxima were collected with the HARPS spectrograph, extending the time span by 549 days with respect to the dataset analyzed in ( 9 ).

Although we used the same recipe for both UVES and HARPS, we noted that an offset between the two datasets still exists by comparing the Hα values taken at a similar epoch (BJD = 2,453,207) and at two consecutive epochs (HARPS, BJD = 2,453,812; UVES, BJD = 2,453,813). It is reasonable to expect an instrumental offset when combining data extracted with the same method from spectra collected with different instruments. To produce a complete UVES + HARPS dataset free from offset, we subtracted the average value 1.26 from the UVES Hα index dataset, as determined from measurements at the epochs indicated above.

We excluded outliers potentially due to powerful flares through a 3σ clipping of the data, resulting in two epochs removed from the UVES dataset and three epochs removed from the HARPS dataset, therefore with a very limited impact on the analysis. Then, we binned the Hα index extracted from the HARPS spectra on a nightly basis. Last, because there is not one-to-one correspondence between the Hα index and RV datasets (some of the latter missing because the corresponding spectra were discarded by the TERRA pipeline), we searched for and selected those epochs that are in common between the two time series to perform a correct correlation study (this caused 10 RVs to be excluded from the correlation analysis).

In addition to photometry, we inspected activity indicators extracted from spectra as the full width at half maximum (only for HARPS) and those based on the chromospheric CaII H + K and Hα emission lines. Because the S/N corresponding to CaII H + K (as measured from HARPS spectra) is generally less than 1 (median value, 0.5), we selected only the Hα index for a reliable analysis. A key point to address to search for long-term, periodic modulations due to activity (to be eventually compared with any significant long-period signal found in the RVs) is to deal with indexes extracted homogeneously both from the UVES and HARPS spectra, covering the whole time span of the observations. To do so, we used the UVES activity indexes already published in ( 9 ), and then we extracted the Hα index from HARPS spectra following the recipe used in ( 9 ) by adapting the code ACTIN ( 18 ). The time series of the Hα index measured from HARPS spectra is listed in table S3.

In this work, we used a long time span dataset of ASAS-3 and ASAS-4 V-band observations ( 16 ). The light curve was an extension of that shown in ( Fig. 3 ) ( 17 ), with a time span increased by 1343 days. We used magnitudes measured with the photometric aperture-labeled MAG2 (appropriate to brightness and crowding of the field) in the ASAS data file and considered only high-quality data (flags A and B). Last, we binned the data on a nightly basis. Their dispersion is 0.044 mag, and the median uncertainty is 0.046 mag. The data are listed in table S2.

The RV dataset extracted from the UVES spectra, collected between 2000 and 2007, is an improved version of that used in ( 9 ), obtained after changes have been applied to the pipeline for processing all the spectra homogeneously. We reanalyzed the entire dataset starting from the raw images and using the associated calibration frames. We reduced the raw images to one-dimensional spectra with our custom reduction package and generated velocities from our precision velocity package. After this full reduction, the root mean square (RMS) of the nightly binned RVs of Proxima has been reduced from 2.30 to 2.02 m/s ( 15 ). The final dataset that we analyzed consists of 202 HARPS RVs (88 before the fiber upgrade) and 77 nightly binned UVES RVs, with a time span of 6392 days.

The RVs from HARPS spectra were extracted using the TERRA (Template-Enhanced Radial velocity Re-analysis Application) pipeline ( 13 ) and represent an updated dataset with respect to that published in ( 9 ). As in previous works, all observations taken before (various programs) and after [including the Pale Red Dot 2016 ( 9 ) and Red Dots 2017 campaigns] the HARPS fiber upgrade in May 2015 were treated as coming from a separate instrument to account for the reported offset introduced by the fiber change. HARPS/ESO (European Southern Observatory) reduced spectra have a known “residual” systematic effect with a ~1-year periodicity caused by a small pixel size difference every 512 pixels on the detector, often called the “stitching problem,” coupled to the barycentric motion of the Earth, which implies that some spectral lines go across these pixels ( 14 ). As detailed in ( 12 ), we masked ±40 pixels around each of these 512 stitches and rerun the RV velocity measurements for both pre- and post-fiber upgrade datasets. Despite the fact that this removes about 10% of the useful Doppler pixels, a bit higher random noise is desirable compared to systematic excursions beating with the yearly sampling. All barycentric corrections were applied as in ( 12 ), and secular geometric acceleration was also removed from the final RVs using the known astrometry of the star (DR2). Because we are interested in testing for the presence of longer period signals, we computed nightly weighted means and added 1 m s −1 in quadrature to the formal errors given by the pipeline to account for some unrealistically small uncertainties in some high signal-to-noise (S/N) spectra.

RESULTS

We analyzed the enlarged RV dataset spanning ~17 years by performing Monte Carlo (MC) analyses in a Bayesian framework using models based on Gaussian process (GP) regression, as described in detail in the Supplementary Materials. Initially, our model included only the orbital equation of the planet Proxima b combined with the GP term describing the stellar activity contribution to the RVs. The best-fit values of the one-planet model parameters are shown in table S1. Then, we subtracted from the complete RV time series the best-fit solution for planet b (eccentric orbit), a secular acceleration term, and the RV offsets (thus, without removing any activity-related signal), and we analyzed the generalized Lomb-Scargle (GLS) periodogram (19) of the RV residuals. We found a clear peak with the highest power at P∼1907 days, with a false alarm probability of 0.01% as derived by a bootstrap (with replacement) analysis of 10,000 randomly generated RV samples (Fig. 1).

Fig. 1 Frequency analysis of the RV residuals. Top: RV residuals after subtracting from the original dataset the spectroscopic signal induced by Proxima b, the instrumental offsets, and a secular acceleration term, as fitted by a global model including a GP quasi-periodic term and only the eccentric orbital equation for Proxima b. The residuals still include a stellar activity term. The red line corresponds to the best-fit sinusoid as derived with GLS (P = 1907 days). Middle and bottom: GLS and BGLS periodograms of the residuals. For the GLS periodogram, we calculated the false alarm probability thresholds, indicated by the dashed horizontal lines, through a bootstrap analysis. For clarity, the inset plot shows a zoom-in view of the low-frequency region, with the highest peak at P = 1907 days marked by a vertical dotted line. The second highest peak in both periodograms occurs at P ~ 307 days, which is the 1-year alias of the candidate planetary signal. Bottom: Window function of the RV time series. The inset plot shows a zoom-in view of the low-frequency region, with the period P = 1907 days marked by a dotted vertical line.

Evolution of the long-period signal over time Before investigating the nature of this signal in detail, we checked its evolution and stability over time by analyzing the stacked Bayesian GLS (SBGLS) (20) periodogram of the one-planet RV residuals (we show the BGLS periodogram of the full dataset in Fig. 1). SBGLS (available online at the Web page https://anneliesmortier.wordpress.com/sbgls/) allows us to check the variability and incoherence with time of a signal—both properties expected to be detected if it is due to the stellar activity—as it calculates BGLS periodograms for subsets of RV data by adding sequentially one data point per time, until the whole dataset is analyzed. We show the SBGLS periodogram for the RV residuals in Fig. 2. It can be seen that the long-period signal is emerging after ~170 RV measurements. We note that the signal appearing at P ~ 307 days is the 1-year alias of the ~1900-day signal, whose significance decreases after N ~ 260 RV measurements, while that of the ~1900-day signal increases. From the logarithm of the posterior probabilities, we can calculate the relative probability of the two peaks in the BGLS periodogram for the full RV dataset. We find that the period of the candidate planet signal is at least ∼1013 times (log[P 1 /P 2 ] ∼ 13) more probable than the following highest peak. The bottom panel of Fig. 2 shows the evolution, with the increasing number of RVs over time, of the values of the orbital period, semi-amplitude, and phase of the candidate planetary signal, obtained from a least squares fit (we remark that the stellar activity contribution has not been removed from these RVs). We see, in particular, that the semi-amplitude and phase of the long-period signal remain constant within the error bars, indicating that it is stable and coherent over time for N > 180 measurements. Fig. 2 Stability and coherence of the long-period signal in the RVs. Top: SBGLS periodograms of the RV residuals (same data as in Fig. 1). Middle: Zoomed-in view of the top plot, starting from the 150th RV measurement. The vertical dashed line marks the orbital period P ~ 1900 days of the candidate planet to guide the eye. Bottom: Evolution of the orbital period, semi-amplitude, and phase of the candidate planet signal with increasing number of RV points, as calculated by GLS through a least squares fit. We also investigated the question of whether just one dataset (HARPS or UVES) mainly contributes to the appearance of the ∼1900-day signal and checked how the periodograms change by removing groups of data from the analysis. We calculated the SBGLS periodograms for four subsets of data, and the results are shown in fig. S1. Although both datasets cover the time span of the proposed signal, we note that after excluding the UVES dataset (panel A), the candidate planetary signal does not emerge clearly with HARPS data only. This result could be influenced by the fact that the UVES dataset alone covers a time span of ~2500 days, which is longer than the ~1900-day signal under investigation, and by the higher precision of the UVES RVs compared to that of HARPS (median σ RV = 0.6 and 1.3 m/s, respectively). However, given the relative nonuniformity of the HARPS datasets, we investigated the sensitivity of the signal to the different HARPS datasets and find that the signal is detected (fig. S1) without including HARPS data from 2016 (panel B, where the ~1900-day period is the most significant) and, with higher probability, excluding data from 2017 (panel C, the probability being similar to that of the 1-year alias at 307 days) or those collected in the time interval 2011–2014 (panel D). Together, these results indicate that the signal with period of ~1900 days appears significant in our data.

RV modeling with activity and planetary signals We investigated the nature of this signal by introducing a second orbital equation in the global model (GP + Keplerians) and running a new analysis. As a first exploratory run, we adopted large uniform priors on the orbital period of the second Keplerian (20 to 6500 days, the upper limit corresponding nearly to the time span of the data) and on the GP hyperparameters λ (0 to 10,000 days) and θ (70 to 100 days), representing the evolutionary timescale of the stellar active regions and the stellar rotation period, respectively. This analysis showed that the majority of the samples piled up around the orbital period P c ∼1900 days (fig. S2), confirming the outcomes of the GLS/BGLS periodograms; i.e., the existence in the data of a second coherent signal as shown by the existence of well-localized solutions for the time of inferior conjunction T c, conj , each spaced by ~1900 days. The best-fit semi-amplitude of the additional Keplerian orbit is K c ∼1.3 m s−1. By assuming that this signal is real and can be modeled reliably by a Keplerian, we inferred the final set of best-fit values of the fitted and derived system parameters through a second run of MC analyses, this time using restricted intervals for some of the priors (which are still treated as uniform/noninformative) and testing two different models (i.e., both planets on circular or eccentric orbits). We get e b = 0.17 − 0.10 + 0.12 (e b <0.22 at 68.3% level of confidence) and e c = 0.41 − 0.26 + 0.34 (e c <0.58 at 68.3% level of confidence) for the eccentricities of planets b and c, respectively. Both have a level of significance less than the 2.45σ threshold derived in (21); thus, the results do not allow us to constrain the eccentricities. Moreover, the comparison of the Bayesian evidences do not favor the eccentric model [ ln (Z circ /Z ecc )∼ +0.8]; therefore, the orbits can be assumed being circular according to our data. However, we note that the median of the posteriors for e b equals the maximum a posteriori value, suggesting a real nonzero eccentricity for planet b, and that our estimate e b = 0.17 is not too different from e b = 0.25 given in (22). The posterior distributions of the free parameters of our model with two circular planets are shown in fig. S3. Our adopted best-fit solution is summarized in Table 1, and Fig. 3 shows the phase-folded RV curves for planet b and candidate planet c. We note that our fit resulted in small values of the uncorrelated jitter for each instrument, and instrumental offsets are consistent with zero. The results for the model with two planets on eccentric orbits are summarized in table S1. According to the adopted model, the candidate planet Proxima c has orbital period P c = 1900 − 82 + 96 days (corresponding to more than three complete orbits over the time span of the observations), semi-major axis a c = 1.48 ± 0.08 AU, minimum mass m c sin i c = 5.8 ± 1.9 M ⊕ , and equilibrium temperature T eq = 39 − 18 + 16 K. The difference between the Bayesian evidences of the two-planet circular model and that including only one circular planet (for which we used the same priors for all the parameters in common) is Δln Z = 1.6, corresponding to a weak-to-moderate evidence in favor of the two-planet model according to the scale given in (23). We note that the Bayesian evidence does not favor the two-planet circular model over the one-planet circular model when a much larger prior on P c is adopted [𝒰(20–6500) days]. In that case, Δln Z∼ −5.8. The sensitivity of the statistical evidence from the choice of the prior range for P c implies that the data place weak constraints on the model evidences and that additional observations are necessary over the next years to cover more orbits of the candidate companion and make the Bayesian statistics less dependent from the prior range in a GP framework. However, as discussed above, the existence of the ~1900-day signal in our present dataset appears justified. We show in fig. S4 the best-fit solution for the stellar activity term as fitted by the GP regression. Table 1 Results of the GP regression analysis applied to RVs, including two circular orbital equations, and to ASAS photometry. View this table: Fig. 3 Phase-folded spectroscopic orbits for Proxima b and c. Top: RV curves of Proxima b and of the candidate planet Proxima c, phase-folded to the orbital periods listed in Table 1. The red curves represent the best-fit orbital solutions, and the red points are phase-binned RV values. Bottom: Distributions of the number of measurements along the planets’ orbits. We also tested a more complex model including an additional circular Keplerian orbit, with the orbital period of a possible third companion explored up to 1600 days (using uniform priors in both linear and logarithmic scales). We did not find evidence for any significant additional signal in the data. We only note that the posterior for the orbital period for the third Keplerian gives hints for a signal at ∼240 days, which we also detected in the periodogram of the RV residuals (after removing only the signals of the two planets; see fig. S4) and that of the Hα activity indicator (see the next section). Moreover, we also used a different model based on the sum of sinusoids to treat the stellar activity. Details and results are described in the Supplementary Materials.