A low-mass planet candidate orbiting Proxima Centauri at a distance of 1.5 AU

See allHide authors and affiliations

Science Advances  15 Jan 2020:
Vol. 6, no. 3, eaax7467
DOI: 10.1126/sciadv.aax7467


Our nearest neighbor, Proxima Centauri, hosts a temperate terrestrial planet. We detected in radial velocities evidence of a possible second planet with minimum mass mc sin ic = 5.8 ± 1.9M and orbital period Pc=5.210.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.


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 ≥ 4MJupiter) (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 MJupiter 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 MJupiter. Using the radial velocity (RV) technique, Endl and Kürster (8) excluded the presence of planets with minimum masses greater than 1 MNeptune 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 (, 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).

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.


RV extraction

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.

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.


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.

Spectroscopic activity diagnostics

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.

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).

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 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 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[P1/P2] ∼ 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 Pc ∼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 Tc, conj, each spaced by ~1900 days. The best-fit semi-amplitude of the additional Keplerian orbit is Kc ∼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 eb=0.170.10+0.12 (eb<0.22 at 68.3% level of confidence) and ec=0.410.26+0.34 (ec<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 (Zcirc/Zecc)∼ +0.8]; therefore, the orbits can be assumed being circular according to our data. However, we note that the median of the posteriors for eb equals the maximum a posteriori value, suggesting a real nonzero eccentricity for planet b, and that our estimate eb = 0.17 is not too different from eb = 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 Pc=190082+96 days (corresponding to more than three complete orbits over the time span of the observations), semi-major axis ac = 1.48 ± 0.08 AU, minimum mass mc sin ic = 5.8 ± 1.9 M, and equilibrium temperature Teq=3918+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 Pc 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 Pc 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.

Is the candidate planetary signal actually linked to the stellar activity?

To investigate whether the signal can be attributed to a planetary companion, or it is possibly due to a long-term stellar magnetic activity cycle, we searched ancillary data for evidence of a periodicity close to ~1900 days. The best dataset for our purpose is represented by the ASAS-3–ASAS-4 V-band light curve, first analyzed in (17). In this work, we used a larger dataset, as detailed in the Supplementary Materials, which now covers 6688 days. This very extended ASAS light curve represents an invaluable dataset, because it allows to study the photometric evolution of Proxima activity over almost the whole time span of the RVs. The ASAS photometric time series is shown in Fig. 4 (top). With respect to the data analyzed in (17), the star’s brightness is observed increasing starting around the epoch HJD 2,457,500. We reassessed the average periodicity of the activity cycle identified in (17) with the older light curve (2576 ± 52 days) by performing an MC fit, which takes into account both the rotational and long-period activity cycle modulations. We used a model composed of two sinusoidal functions (one modeling the rotation and one modeling the activity cycle) and of one quadratic term to model the long-term trend seen in Fig. 4. For the rotation period Prot, we used a uniform prior in the range of 80 to 90 days, while for the average activity cycle period Pact. cycle, we used a uniform prior in the large range of 1000 to 6500 days. The best-fit model is represented by the red curve in the top plot of Fig. 4, corresponding to an average period of the activity cycle Pact.cycle=238244+47 days and to a rotation period Prot = 83.12 ±0.07 days. In the bottom panel of Fig. 4, we show the ASAS light curve, after removing the 83-day rotational signal and the quadratic long-term trend, folded at the best-fit period of the activity cycle. Our newly determined value for the mean Pact. cycle is nearly 200 days shorter than that estimated in (17). It differs from our best-fit orbital period for the candidate planet c (Pc=190082+96 days) by 482 ± 107 days, i.e., the two periods differ by 4.5σ. Even if these results could be influenced by the different sampling of the datasets, the currently available ASAS dataset indicates that the period of our candidate planetary signal is significantly distinct from that of the activity cycle; therefore, it presently does not support the interpretation of the ∼1900-day signal in terms of stellar activity. This conclusion certainly needs a more robust confirmation through a photometric follow-up extended over the next years. We also performed a GP regression analysis of the ASAS light curve by adopting the same quasi-periodic kernel used for the RV time series and large priors (Table 1), without including any other term in the model. Our results show that the stellar rotation period θ is well recovered and the timescale of the correlations λ attains a value close to that estimated for the RVs. This suggests that photometry and RVs contain quasi-periodic activity-related signals with similar properties, and this is particularly suggestive by taking into account that the datasets have different length, time span, and sampling.

Fig. 4 Analysis of the light curve of Proxima Centauri.

Top: ASAS light curve of Proxima (gray dots). The red curve represents our best-fit model of the photometric data, which includes two sinusoids (for the rotational and activity cycle modulations) and a quadratic term to take into account the rise in brightness particularly evident after the epoch HJD 2,457,500. Bottom: ASAS light curve, after removing the 83-day rotational signal and the quadratic long-term trend, folded at the best-fit period of the activity cycle. The best-fit sinusoid modeling the activity cycle is represented by the red curve. The epoch corresponding to phase 0 is HJD 2,458,049.79

Collins et al. (24) analyzed 13 years of spectroscopic and photometric observations and derived 82.1 days for the stellar rotation period of Proxima based on Hα activity index, in agreement with the period of 83 days calculated in (4) with Hubble Space Telescope (HST) photometry and confirmed in (17, 25). Robertson et al. (25) did not find anything relevant on longer timescales, likely as a consequence of the stochastic nature of the microflaring activity of Proxima, which dominates the Hα line emission and makes the analysis of this index very complex in the time domain, as already pointed out in (8). We performed an analysis of activity diagnostics extracted from UVES and HARPS spectra to search for evidence of a long-term activity cycle and correlations with the RV time series that could explain the candidate planetary signal alternatively in terms of activity. Concerning indexes derived from chromospheric emission lines, we consider here only that derived from the Hα line (see Materials and Methods), for which the median S/N is 52 compared to a median S/N of 0.5 of the index based on the CaII H + K lines. Figure S5 shows the time series of the Hα index, the GLS periodogram, and correlation diagrams with the RV residuals (after subtracting from the original dataset the spectroscopic signal induced by Proxima b, the instrumental offsets, and a secular acceleration term, thus still including activity-related signals). The main peak in the periodogram occurs at 236 days, and no long-term activity cycle is detected, neither in the original dataset nor in the residuals after pre-whitening the data by subtracting the 236-day signal. The same result is obtained when analyzing data extracted from HARPS spectra only. The correlation between the Hα index and the RV residuals is not significant for both the complete UVES + HARPS and only HARPS datasets (Spearman’s correlation coefficients are ρ = 0.20 and ρ = 0.23, respectively).

Very recently, Pavlenko et al. (26) analyzed the temporal variations of some chromospheric emission lines from the same HARPS optical spectra of our sample. Their analysis focused on the pseudo-equivalent widths and profiles of the emission lines, and they found that all the lines show short timescale variability (at least 10 min). We do not report about long-period activity cycles detected in the time series of the analyzed indexes. We calculated the GLS periodogram of their Hα index dataset (after removing five outliers) and found that the main peak is located at 234 days, with no evidence of significant signals with long periods, i.e., a result very similar to that obtained with our dataset.

In conclusion, our analysis of activity diagnostics, as the long time span photometric and the Hα index datasets, did not reveal the existence of a periodicity similar to that of our 5.2-year candidate planet and of correlations that would explain the long-period signal detected in the newly extracted RV dataset in terms of stellar activity. Therefore, we affirm to the best of our knowledge that the signal may be due to an additional planet in the system, Proxima c. Nonetheless, the lack of correlation between Hα index and RVs does not entirely exclude that the ~1900-day signal is linked to the activity. Still, little is known about the impact of activity cycles on the RV measurements for slowly rotating cool dwarfs as Proxima. Moreover, because the statistical significance of our GP + 2 planets model is not high with the present dataset, and due to the long period and small semi-amplitude of the signal, its real nature needs to be further investigated with additional extensive RV and photometric follow-up, by using alternative models and data analysis methods to treat the stellar activity, and with different detection techniques. How we will show in the next section, a decisive contribution to this regard is expected to come from Gaia astrometric observations.


The semi-amplitude of the candidate’s Doppler signal is equal to that of Proxima b, and it is significant at a 3σ level. By adopting the metric K/N=Kplanet × Nepochs/RMSdataset defined in (27) as the figure of merit to attest the veracity of planetary signals detected in RV measurements (where Nepochs = 279 and RMS = 2.2 m s−1 as measured from the combined UVES/HARPS dataset), we get K/N ∼10, which is above the K/N = 7.5 threshold proposed in (27).

Proxima was observed with the Atacama Large Millimeter/submillimeter Array (ALMA) at a wavelength of 1.3 mm by Anglada et al. (28), who reported the existence of an unknown source at a projected distance of about 1.2 arc sec from the star (equivalent to 1.6 AU). We note here the similarity of our derived orbital semi-major axis of the candidate planet with that determined in (28) for the point-like source in their ALMA images. Nonetheless, Anglada et al. (28) could not rule out the possibility of it being a background galaxy or a transient phenomenon. ALMA imaging could corroborate the existence of Proxima c if the secondary 1.3-mm source is confirmed: In this sense, ALMA follow-up observations will be essential. In (28), the possible existence of a cold dust belt at ∼30 AU, with inclination of 45°, is also mentioned. If Proxima c orbits on the same plane, its real mass would be mc = 8.2 M and both planets in the system would fall in the range of super-Earths, making Proxima the nearest system of multiple super-Earths to the Sun.

Considerations about planet formation and evolution

The existence of Proxima c is highly significant for planet formation models. This planet would be the one with the longest period and a minimum mass in the super-Earth regime presently detected with the RV technique around a low-mass star (fig. S6). It would also be the first at a distance from the parent star much larger than the expected original location of the snowline in the protoplanetary disk, which was within 0.15 AU according to (29) (see the Supplementary Materials for a detailed discussion). Microlensing observations have hinted at the existence of super-Earths at distances of ∼1 AU from 0.1 M stars (e.g., MOA‑2010-BLG-328L b) but with large uncertainties. It is unlikely that Proxima c has been kicked out during a planet system instability from an initial position much closer to the star, because its orbit is consistent with a circular one and because of the absence of more massive planets on shorter orbital distance. The formation of a super-Earth well beyond the snowline challenges formation models according to which the snowline is a sweet spot for the accretion of super-Earths, due to the accumulation of icy solids at that location (30, 31), or it suggests that the protoplanetary disk was much warmer than usually thought (29, 32). The planet is massive enough relative to the central star to have opened a relatively deep gap in the protoplanetary disk and have migrated in type II mode. According to Kanagawa et al. (33), its migration timescale would have been 1 million years (see the Supplementary Materials).

Thus, the planet might have formed at a few times its current distance. This can be linked with a possible inner ring at 1 to 4 AU, proposed in (28) from ALMA 1.3-mm data, although debated in (34). We speculate that this inner ring may be due to dust produced by planetesimals clustered in one of the planet’s inner mean motion resonances. This clustering might have happened during the planet’s inward migration. To this regard, new observations with ALMA will be fundamental for characterizing the system.

The crucial role of Gaia astrometry

A major role in confirming the existence of Proxima c will be played by space-based astrometry. Using the Hipparcos catalog and Gaia DR2, Kervella et al. (7) detected an anomaly (i.e., deviation from a purely linear tangential proper motion) in Proxima’s tangential velocity significant at a 1.8σ level. If confirmed, this is compatible with the existence of a planet with true mass ~10 to 20 M at a distance of ~1 to 2 AU [see figure 14 in (7)]. We note that our estimates for the minimum mass and orbital radius of the candidate planet Proxima c are in good agreement with the ranges calculated in (7). An independent estimate of Proxima’s proper motion variation between the Hipparcos and Gaia DR2 epochs is found in (35), which confirms the existence of an anomaly significant at a 5.3σ level. To further support the presence of the proper motion anomaly, we calculated the parameter ΔQ between the Hipparcos and Gaia proper motions, as defined in equation 10 of (36). We found ΔQ = 19, which indicates significant deviation from linear tangential proper motion between the two epochs. To assess whether this anomaly can be explained by the candidate planet alone, we performed a Markov chain Monte Carlo analysis of the observed proper motion variation, as computed in (35), by fitting for the orbital inclination i and the ascending node longitude Ω, while keeping the orbital parameters derived from RVs fixed. Because of the few measurements available, and of the large difference between the precisions of Hipparcos and Gaia proper motions (nearly one order of magnitude), we cannot yet state whether the observed Δμ is due to Proxima c, pointing out the need for measurements at multiple epochs to possibly reach a conclusive analysis.

Given the target brightness and the expected minimum size of the astrometric signature (αsin ic=16746+47 μas), Gaia alone should clearly detect the astrometric signal of the candidate planet at the end of the 5-year nominal mission, all the more so in case of a true inclination angle significantly less than 90°. Proxima is one of the very few stars in the Sun’s backyard for which Gaia alone might be sensitive to an intermediate separation planetary companion in the super-Earth mass regime. We carried out a detailed numerical experiment combining the available RVs with synthetic, but realistic, Gaia observations of Proxima more than the 5-year nominal mission (see the Supplementary Materials), which spans almost one orbit of planet c. Assuming the orbital parameters derived from the RV analysis, we explored the range of allowed inclinations (and thus actual mass values for Proxima c) compatible with the 1σ upper limits of 0.1 to 0.2 MJ from (7). The results of the analysis (Fig. 5) indicate that the true mass of Proxima c should be measurable with a typical accuracy of 17% with astrometry only, while combining astrometry with RVs increases the accuracy to 5%. From the simulations, we infer that the combined RV + astrometry solution will allow a very accurate reconstruction of the full set of orbital elements of Proxima c. We note that these results should be considered as a best-case scenario, provided the Gaia per-measurement astrometric performance will be that expected (37) on a star with Proxima’s magnitude at G-band (G = 8.9 mag). Given the encouraging results, even in the case of degraded precision, Gaia is expected to play a crucial role for the improved characterization of Proxima c.

Fig. 5 Outcomes of the combined analysis of the astrometric and RV datasets.

Left: True mass of Proxima c versus the sine of the orbital inclination, as obtained from the astrometric simulations. The black line is the simulated exact solution, the blue dots represent the values derived from the Gaia astrometry alone, while the red dots are the values derived by combining the Gaia astrometry with the radial velocities. Right: Fractional deviation of the true mass (defined as the difference between the simulated and retrieved masses for Proxima c divided by the simulated value) versus sine of the orbital inclination.

Final considerations

The precise mass estimate will, in turn, permit to resolve, at least in part, important model degeneracies in predictions of Proxima c’s apparent brightness in reflected visible light due to orbit geometry, companion mass, system age, orbital phase, cloud cover, scattering mechanisms, and degree of polarization. The flux contrast between Proxima c and the parent star for thermal emission and for reflected light is between 10−8 and 10−9, depending on the geometric albedo. These values are beyond the capabilities of presently available instruments for direct imaging, but, given the apparent maximum separation of ∼1 arc sec, future high-contrast imaging instrumentation [similar to that on the European Extremely Large Telescope (38), among several other ground and space-based facilities] will be capable of observing the object for a large fraction of its orbit. The combined RV + astrometry solution will therefore make possible detailed ephemerides predictions to optimize the planning and interpretation of follow-up/characterization measurements of Proxima c with such instrumentation. Although Proxima c represents a very challenging target for the combination SPHERE + ESPRESSO envisioned in (39) for detecting the atmosphere of Proxima b with high-contrast imaging, it will be the most appealing target for similar aggregates of instruments mounted on extremely large telescopes.


Supplementary material for this article is available at

Fig. S1. Stacked BGLS periodograms for subsets of RV residuals.

Fig. S2. Posterior distributions of parameters and hyperparameters relative to an MC exploratory run where the global model is composed of a GP quasi-periodic kernel and two planetary orbital equations, with eccentricities treated as free parameters.

Fig. S3. Posterior distributions for all the free parameters of the model adopted in this work, which includes a GP quasi-periodic kernel to fit the stellar activity term in the RV, and two circular orbital equations (see Table 1).

Fig. S4. Stellar activity signal as found in the RV through a GP regression using a quasi-periodic kernel.

Fig. S5. Analysis of Ha activity index.

Fig. S6. Minimum mass versus orbital semi-major axis of confirmed planets orbiting low-mass stars (M < 0:6 M) discovered with the RV technique.

Table S1. Other models used for fitting the RVs, as discussed in the text.

References (4051)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We made use of data from the European Space Agency (ESA) mission Gaia (, processed by the Gaia Data Processing and Analysis Consortium (DPAC, Funding for the DPAC has been provided by national institutions, in particular, the institutions participating in the Gaia Multilateral Agreement. We acknowledge the computing centers of INAF-Osservatorio Astronomico di Trieste/Osservatorio Astrofisico di Catania, under the coordination of the CHIPP project, for the availability of computing resources and support. We acknowledge B. Wargelin for useful discussion on the ASAS photometry; N. Lanza, A. Mortier, M. Pinamonti, F. Bauer, M. MacGregor, and P. Kalas for interesting and useful discussions on various aspects of this work; and M. Tuomi, E. Palle, A. Reiners, M. Zechmeister, S. Dreizler, Y. Tsapras, J. Barnes, A. Ofir, S. Jeffers, J. Morin, R. Nelson, G. Coleman, I. Ribas, V. J. Sanchez Bejar, L. Tal-Or, J. B. P. Strachan, C. Haswell, and Z. M. Berdiãs for their contribution to the RD campaign. F.D.S. acknowledges Artemis Del Sordo for unique conversations. Last, we are grateful to T. Milian, F. Lechner, and M. Brega for having contributed to this work with their cheerfulness, from somewhere close to Proxima. Quoting Henry David Thoreau: “If you have built castles in the air, your work need not be lost; that is where they should be. Now put the foundations under them.” Funding: M.D. acknowledges financial support from Progetto Premiale 2015 FRONTIERA funding scheme of the Italian Ministry of Education, University, and Research. P.G. acknowledges financial support from the Italian Space Agency (ASI) under contract 2014-025-R.1.2015 to INAF. D.B. acknowledges financial support from INAF and Agenzia Spaziale Italiana (ASI grant no. 2014-025-R.1.2015) for the 2016 PhD fellowship programme of INAF. H.R.A.J. is supported by UK STFC grant ST/R006598/1. J.S.S. and P.A.P.R. acknowledge support by FONDECYT grant 1161218 and partial support from CONICYT project Basal AFB-170002. M.J.L.-G., N.M., C.R.-L., E.R., G.A., and J.F.G. acknowledge financial support from Spanish MINECO AYA2016-79425-C3-3-P, AYA2017-84390-C2-1-R (co-funded by FEDER), ESP2017-87143R, and ESP2017-87676-C05-02-R and the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709). G.A.-E. research is funded via grant ST/P000592/1 “Astronomy Research at Queen Mary” granted by the Science and Technology Facilities Council/United Kingdom. Author contributions: M.D. conceived this work and performed the analyses of RVs and photometry, detecting and characterizing the candidate Doppler signal. He led the preparation of the manuscript. F.D.S. conceived this work and contributed to the setup of the analyses and interpretation of the results. He contributed to writing the paper. G.A.-E. is the principal investigator of the RD campaign. He extracted the HARPS RVs, collaborated to the analyses of the RVs, and participated in the manuscript writing. P.G. and A.S. conceived and performed the astrometric analyses and contributed to writing the paper. A.M. wrote the section concerning the formation and orbital evolution of Proxima’s system. G.P. provided new photometric data of Proxima from the ASAS survey. D.B. contributed to the astrometric analyses. R.P.B., H.R.A.J., and F.F. worked on the extraction of the UVES RVs and contributed to the manuscript writing. F.-J.H., J.S.J., M.J.L.-G., N.M., P.A.P.R., E.R., and C.R.-L. contributed with photometric follow-up observations of Proxima Centauri within the RD collaboration. P.J.A., G.A., and J.F.G. contributed to the manuscript writing with relevant discussion about searching for the candidate counterpart in ALMA data and preparing follow-up observations. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. All the datasets analyzed in this work (radial velocities, Hα activity index time series, and ASAS photometry) are additionally made available at the Web address

Stay Connected to Science Advances

Navigate This Article