## Abstract

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.9*M*_{⊕} and orbital period Pc=5.21+0.26−0.22 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.

## INTRODUCTION

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

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.

## MATERIALS AND METHODS

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

### Photometry

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.

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

### 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 ∼10^{13} 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.

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 eb=0.17+0.12−0.10 (*e*_{b}<0.22 at 68.3% level of confidence) and ec=0.41+0.34−0.26 (*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 Pc=1900+96−82 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 Teq=39+16−18 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.

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 *P*_{rot}, we used a uniform prior in the range of 80 to 90 days, while for the average activity cycle period *P*_{act. 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=2382+47−44 days and to a rotation period *P*_{rot} = 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 *P*_{act. 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=1900+96−82 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.

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.

## DISCUSSION

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 *N*_{epochs} = 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 *m*_{c} = 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=167+47−46 μ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 *M*_{J} 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.

### 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 MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/3/eaax7467/DC1

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.

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.

Fonte: *Science Advances*