Issue 
A&A
Volume 658, February 2022



Article Number  A115  
Number of page(s)  16  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202142337  
Published online  10 February 2022 
A candidate shortperiod subEarth orbiting Proxima Centauri^{★}
^{1}
Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas,
4150762
Porto,
Portugal
email: joao.faria@astro.up.pt
^{2}
Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua Campo Alegre,
4169007
Porto,
Portugal
^{3}
Instituto de Astrofísica de Canarias,
38205
La Laguna,
Tenerife,
Spain
^{4}
Departamento de Astrofísica, Universidad de La Laguna,
38206
La Laguna,
Tenerife,
Spain
^{5}
European Southern Observatory,
Alonso de Cordova 3107,
Vitacura,
Santiago,
Chile
^{6}
INAF – Osservatorio Astrofisico di Torino,
Via Osservatorio 20,
10025
Pino Torinese,
Italy
^{7}
Département d’astronomie de l’Université de Genève,
Chemin Pegasi 51,
1290
Versoix,
Switzerland
^{8}
Consejo Superior de Investigaciones Científicas,
28006
Madrid,
Spain
^{9}
INAF – Osservatorio Astronomico di Trieste,
Via Tiepolo 11,
34143
Trieste,
Italy
^{10}
Physics Institute of University of Bern,
Gesellschaftsstrasse 6,
3012
Bern,
Switzerland
^{11}
Department of Physics, and Institute for Research on Exoplanets, Université de Montréal,
Montréal,
H3T 1J4,
Canada
^{12}
Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências da Universidade de Lisboa,
Edifício C8, Campo Grande,
1749016
Lisbon,
Portugal
^{13}
Faculdade de Ciências da Universidade de Lisboa (Departamento de Física),
Edifício C8,
1749016
Lisboa,
Portugal
^{14}
Institute for Fundamental Physics of the Universe, IFPU,
Via Beirut 2,
34151
Grignano,
Trieste,
Italy
^{15}
Scuola Normale Superiore, Piazza dei Cavalieri,
7 56126
Pisa,
Italy
^{16}
Centro de Astrobiología (CAB, CSICINTA), Depto. de Astrofísica, ESAC campus,
28692
Villanueva de la Cañada (Madrid),
Spain
^{17}
European Southern Observatory,
KarlSchwarzschildStr. 2,
85748
Garching bei München,
Germany
^{18}
Centro de Astrofísica da Universidade do Porto,
Rua das Estrelas,
4150762
Porto,
Portugal
^{19}
INAF – Osservatorio Astronomico di Palermo,
Piazza del Parlamento 1,
90134
Palermo,
Italy
^{20}
INAF – Osservatorio Astronomico di Brera,
Via E. Bianchi 46,
23807
Merate,
Italy
^{21}
Centro de Astrobiología (CSICINTA),
Crta. Ajalvir km 4,
28850
Torrejón de Ardoz,
Madrid,
Spain
Received:
30
September
2021
Accepted:
24
December
2021
Context. Proxima Centauri is the closest star to the Sun. This small, lowmass, mid M dwarf is known to host an Earthmass exoplanet with an orbital period of 11.2 days within the habitable zone, as well as a longperiod planet candidate with an orbital period of close to 5 yr.
Aims. We report on the analysis of a large set of observations taken with the ESPRESSO spectrograph at the VLT aimed at a thorough evaluation of the presence of a third lowmass planetary companion, which started emerging during a previous campaign.
Methods. Radial velocities (RVs) were calculated using both a crosscorrelation function (CCF) and a template matching approach. The RV analysis includes a component to model Proxima’s activity using a Gaussian process (GP). We use the CCF’s full width at half maximum to help constrain the GP, and we study other simultaneous observables as activity indicators in order to assess the nature of any potential RV signals.
Results. We detect a signal at 5.12 ± 0.04 days with a semiamplitude of 39 ± 7 cm s^{−1}. The analysis of subsets of the ESPRESSO data, the activity indicators, and chromatic RVs suggest that this signal is not caused by stellar variability but instead by a planetary companion with a minimum mass of 0.26 ± 0.05 M_{⊕} (about twice the mass of Mars) orbiting at 0.029 au from the star. The orbital eccentricity is well constrained and compatible with a circular orbit.
Key words: techniques: radial velocities / stars: activity / stars: individual: Proxima
© ESO 2022
1 Introduction
The last couple of decades have seen a fast increase in the number of known exoplanets orbiting solartype stars. Several dedicated space missions (e.g. CoRoT: Baglin et al. 2006, Kepler: Borucki 2016, TESS: Ricker et al. 2010, CHEOPS: Broeg et al. 2013) and groundbased instruments (e.g. HARPS: Mayor et al. 2003, CARMENES: Quirrenbach et al. 2010, HARPSN: Cosentino et al. 2012, HPF: Mahadevan et al. 2012 and soon NIRPS: Wildi et al. 2017 and SPIRou: Artigau et al. 2014) are or will soon be pushing the detection limits towards smaller and lowermass planets, and they already allow for the detailed study of the atmospheres and internal compositions of individual planets (e.g. Ehrenreich et al. 2020). In this context, M dwarfs continue to be prime targets for exoplanet searches due to their smaller sizes and masses, resulting in higheramplitude transit and radial velocity (RV) signals from the planets.
Coincidentally, several of the M dwarfs closest to the Sun are known or suspected to host a planetary system (e.g. Ribas et al. 2018; Bonfils et al. 2018; Jeffers et al. 2020; Díaz et al. 2019; Tuomi et al. 2019; LilloBox et al. 2020). In particular, the discovery of an Earthmass planet orbiting Proxima Centauri (AngladaEscudé et al. 2016), our closest stellar neighbour, was one of the most significant results in the field, in part because the planet orbits inside the habitable zone (HZ) of the star (e.g. Kopparapu et al. 2013). More recently, a second candidate superEarth was detected at a longer orbital period of 5.21 yr (Damasso et al. 2020) but has so far eluded a clear detection with direct imaging or astrometry (Gratton et al. 2020; Kervella et al. 2020; Benedict & McArthur 2020).
The detection of Proxima b and similar planets orbiting M dwarfs demonstrates that highprecision RV measurements allow us to find Earthmass planets that orbit inside, or close to, the HZ of their host stars. Unfortunately, several photospheric and chromospheric phenomena associated with the presence of active regions at the stellar surface induce RV variations that can mimic or contaminate a planetary signal (e.g. Fischer et al. 2016, and references therein). For M dwarfs in particular, such activityinduced RV variations can reach a few m s^{−1}, even for quiet stars (e.g. Suárez Mascareño et al. 2017). As such, the detection of Earthmass planets imposes a detailed characterisation of stellar activity, and it has motivated the development of precise instrumentation, such as the stateoftheart ESPRESSO spectrograph (Pepe et al. 2021), which combines the large collecting power of the Very Large Telescope (VLT) with a very high level of stability in order to achieve unprecedented RV precision.
One of the first results from ESPRESSO was the independent confirmation of Proxima b (Suárez Mascareño et al. 2020, hereafter SM2020), which also revealed the first hints of a shorterperiod signal that could be caused by a lowmass planet. Here we report on followup observations that confirm the presence of this low amplitude signal, which is probably caused by a planet with a minimum mass of only 0.26 M_{⊕} and an orbital period of 5.12 days. We first describe the observations (Sect. 2) and the methods used to derive (Sect. 3) and analyse (Sect. 4) the RVs. In Sect. 5 we assess the evidence for the planetary nature of the signal, and we discuss our results in Sect. 6.
Stellar properties of Proxima, compiled from the literature.
2 ESPRESSO observations
Proxima (see Table 1 for the stellar parameters) is on the target list of the ESPRESSO Guaranteed Time Observations (GTO) survey, which is monitoring a select group of nearby stars with the goal of discovering lowmass planets in the HZ (Hojjatpanah et al. 2019; Pepe et al. 2021). The first planet discovery from the survey was presented very recently in LilloBox et al. (2021). After an initial campaign aimed at confirming Proxima b, (SM2020) we restarted observations of Proxima in order to confirm a candidate signal close to 5 days. In this new campaign, the observational strategy was tailored to this period, with a typical interval of 1–2 days between observations and two exposures per night when possible.
To add to the 67 observations reported in SM2020, we obtained 52 new ESPRESSO spectra of Proxima, for a total of 117 observations spread over 99 individual nights from 20190210 to 20210506. The measurements were taken in ESPRESSO’s high resolution mode (HR21) with an exposure time of 900 s ^{1}. We used the Fabry Pérot (FP) for simultaneous calibration, which allows the instrumental drift to be monitored with a precision better than 10 cm s^{−1} (Wildi et al. 2010).
After a careful analysis of each spectrum, we decided to exclude three observations that were affected by instrumental issues. On the night of 25 April 2019, there was a cooling water temperature increase, which propagated inside the spectrograph. The drift measured on the red detector was 0.007 pixel (i.e. 3 m s^{−1}). This spectrum also shows a signaltonoise ratio below 1. On the observation of 31 July 2019, the FP exposure was saturated, most likely due to an issue with the neutral density filters. In this situation, the contamination by the FP light on the chargecoupled device may be large, and the drift measurement is likely unreliable. Finally, on the night of 10 April 2021, the reported temperature of the atmospheric dispersion compensator (ADC) was 0 K, revealing an issue with the correction. The barycentric Julian days of these three excluded spectra are 2 458 598.523839, 2 458 695.52992, and 2 459 314.583084.
In June 2019, ESPRESSO underwent an intervention to update the fibre link, which improved the instrument’s efficiency by up to 50% (Pepe et al. 2021). This intervention introduced an RV offset, leading us to consider separate ESPRESSO18 and ESPRESSO19 datasets ^{2}. More recently, operations at Paranal were interrupted due to the COVID19 pandemic and ESPRESSO was taken out of operations between 24 March 2020 and 24 December 2020. This led to a large gap in the observations after the initial campaign. Moreover, a change in one of the calibration lamps after the rampup of the instrument is likely to have introduced another RV offset. Therefore, we consider an independent ESPRESSO21 dataset for data obtained after the rampup. In summary, we have 114 available RVs, divided between ESPRESSO18 (50 points), ESPRESSO19 (15 points), and ESPRESSO21 (49 points) subsets. The full time span of the data is 815 days, which we note is about 2.3 times shorter than the orbital period of Proxima c (Damasso et al. 2020).
3 Radial velocity determination
ESPRESSO data were reduced with the instrument’s data reduction software (DRS), version 2.2.8 (Pepe et al. 2021). The DRS transforms the raw images into wavelengthcalibrated 2D spectra (in orderpixel space) through a series of recipes, which includes bias and dark subtraction, optimal extraction of spectral orders, flatfielding and deblazing, wavelength calibration, and computation of the barycentric correction and instrumental drift.
From the extracted spectra, we used two techniques to determine the RVs, one based on the crosscorrelation function (CCF) and the other on a templatematching (TM) algorithm, described below.
3.1 CCF
The ESPRESSO DRS automatically derives the RV of each extracted spectrum using the crosscorrelation technique (e.g. Baranne et al. 1996; Bouchy et al. 2001). In this case, the mask used to build the CCF was derived from an ESPRESSO spectrum of Proxima itself, with a linesearch algorithm based on the derivative of the spectrum. Each spectral order is crosscorrelated individually, producing one CCF per order, which are then combined into a final CCF per spectrum.
To trace the star’s activity, we extract a number of activity indicators from the ESPRESSO spectra or from the CCF. The DRS calculates the CCF’s full width at half maximum (FWHM), contrast (i.e. the relative depth of the CCF), and a shape indicator called CCF asymmetry (see Pepe et al. 2021). Using actin (Gomes da Silva et al. 2018), we also measure activity indices based on the CaII H&K, HeI, Hα, and NaI lines.
3.2 Template matching
We also consider a TM algorithm for RV extraction where each spectra is compared to a template spectrum built from the available ESPRESSO observations. This technique has been shown to improve RV precision when compared with the CCF method, in particular for M dwarfs (e.g. AngladaEscudé & Butler 2012; Zechmeister et al. 2018; Lafarga et al. 2020). A detailed description of the method may be found in Silva et al. (2021). We provide here a brief overview.
As a first step, a telluric template is created from a synthetic spectrum of Earth’s transmittance, built with the TAPAS webinterface (Bertaux et al. 2014). All regions where the absorption drops by more than 1% from the continuum (determined with a median filter) are flagged as telluric features, to be masked from the observed spectra.
To build a high signaltonoise stellar template, all the observed spectra (as extracted by the DRS and before blaze correction) are combined order by order. The individual observations are RV shifted to a common reference frame, arbitrarily chosen to bethat of the observation with the lowest uncertainty on the CCF RV, using the RVs derived previously with the CCF technique. The flux of each order is interpolated to a common wavelength grid and the mean flux of all observations is then calculated. We calculate the mean, instead of the sum, to avoid possible issues with numerical overflow. Uncertainties in the stellar template, originating from both the photonnoise of the individual spectra and from interpolation, are propagated (see Silva et al. 2021, for details). Importantly, we determine different stellar templates for each set of observations from ESPRESSO18, ESPRESSO19, and ESPRESSO21 to prevent any bias from instrumental systematics.
Finally, and unlike other TM approaches, we use a single RV shift to describe simultaneously the differences between all orders of a given spectrum and the template. Within a Bayesian framework, we estimate the posterior distribution for the (single) RV shift, after marginalising with respect to a linear model for the continuum levels of the spectra and template. For computational tractability, we use a Laplace approximation (see e.g. Bleistein & Handelsman 1986, Chap. 8) to characterise the posterior distributions for the RV shifts.
From the Laplace approximation to the posterior distribution we use the mean as the estimated RV and the standard deviation as the estimated RV uncertainty, for each epoch. These uncertainties include the contributions from the photon noise of each observation, the interpolation used in constructing the template, as well as the ordertoorder RV scatter.
For Proxima, this TM algorithm provides RVs that have lower uncertainties than those derived from the CCFs. Figure 1 shows a comparison of the two RV time series (top panels), highlighting the similar RV dispersion (bottom panel, left) but quite different distribution of RV uncertainties (bottom panel, right). The LombScargle periodograms (Lomb 1976; Scargle 1982) of the two sets of RVs (middle panels), which have been calculated after subtraction of the weighted mean of each subset, also show differences. Most noticeably, the TM RVs present more power close to 5 days and close to the first harmonic of the stellar rotation period (~40 days). We also note the correlation coefficients between the two sets of RVs: r = 0.96 (Pearson’s) and ρ = 0.95 (Spearman’s).
In summary, the two sets of RVs are compatible but the TM RVs are more precise. Throughout the paper, we perform the exact same analyses and describe the results for both the CCF RVs and the TM RVs.
4 Radial velocity analysis
4.1 Simple periodogram analysis
We first consider a very simple prewhitening procedure applied to the CCF RVs in order to determine the main periodicities present in the data. The observed RVs are initially modelled with a Gaussian process (GP) using a standard quasiperiodic (QP) kernel (Rasmussen & Williams 2006; Haywood et al. 2014). For simplicity, the hyperparameters of the kernel are fixed to the values determined in SM2020 from a combined fit (their Table 3). Two RV offsets are adjusted between the three data subsets. The resulting fit is shown in the top panel of Fig. 2 and the periodogram of the residuals is shown in panel a) of the same figure. The highest peak in the periodogram is at 11.19 days with a false alarm probability well below 1%.
A sinusoidal signal is then fit to the residual RVs, with a period equal to that of the periodogram peak. The periodogram of the new residuals from this fit is shown in panel b, with a significant peak now appearing close to 5 days and another peak (also significant) at an alias close to 1.2 days. After one additional sinusoidal fit to the residuals, no significant peaks are found in the periodogram (panel c).
The two sinusoids recovered with this sequence of steps have semiamplitudes equal to 85 and 38 cm s^{−1}, respectively, and the final residuals show an rms close to 50 cm s^{−1}. Apart from the stellar activity signals that are being modelled by the initial GP, the 11day and 5day signals are the two main signals present in the RV data. However, this simple procedure is not sufficient to conclude on the planetary nature of the two signals. Moreover, the GP may be absorbing signals that are not caused by stellar activity and biasing the results. In the following section we build a more robust model considering Keplerian signals and perform a simultaneous analysis.
4.2 Joint RV and FWHM model
For a more robust determination of the model parameters, we analyse the RV and FWHM time series jointly, using a shared model for stellar activity. The FWHM was found to be an excellent proxy for the stellar activity of Proxima, showing clear variations at the stellar rotation period (80–90 days), and closely tracing the photometric behaviour of the star (SM2020). The same FWHM time series will be used for the analysis of RVs derived from the CCFs and from the TM approach.
Our stellar activity model includes a GP with most of the hyperparameters shared between RVs and FWHM, but distinct variances. We tested two covariance functions to model stellar activity variations, the QP kernel mentioned above and the quasiperiodic with cosine (QPC) kernel proposed by Perger et al. (2021). The results described below correspond to the QP kernel. Other assumptions for the activity model provide compatible results and are discussed briefly in Appendix A.
To model planetary RV signals, we assume Keplerian orbits parameterised by the orbital period P, semiamplitude K, eccentricity e, mean anomaly M_{0} at the time of the first observation, and argument of pericentre ω. As a simple (and fast) criterion for dynamical stability, we only accept proposed solutions from the prior that are stable according to the angular momentum deficit (AMD) criterion (Laskar & Petit 2017), using the extension proposed by Petit et al. (2017), which takes meanmotion resonances into account.
The subsets of RVs from ESPRESSO18, ESPRESSO19, and ESPRESSO21 are each assigned an independent RV offset and an additional jitter in the form of uncorrelated (white) noise characterised by its variance, and the same for the FWHM. A quadratic trend is also included in the RVs. The full model has a total of 4 + 5 + 2N_{s} + 2(N_{s} − 1) + 5N_{p} parameters, where N_{s} = 3 is the number of data subsets and N_{p} is the number of Keplerians in the model. The different terms correspond to the average RV and FWHM plus the coefficients of the quadratic trend, the hyperparameters of the GP, the jitters, the instrument offsets, and finally the orbital parameters.
Since we ultimately aim at comparing models with a different number of planets, it is important to carefully choose the prior distributions and that they encode prior knowledge independently of the value of N_{p}. The GP amplitudes for the RVs and the FWHM (η_{1}) were assigned modified log uniform priors (e.g. Gregory 2005) up to the span of each dataset (8.8 m s^{−1} for the RVs and 36.4 m s^{−1} for the FWHM) and with a break point at 1 m s^{−1}. For the evolution timescale, η_{2}, we use a log uniform prior between 60 and 400 days and for η_{3}, the parameter associated with the stellar rotation period, we use a uniform prior between 60 and 100 days. Proxima’s rotation period was estimated to be 89.8 days (independently from RV data) to a precision of 4 days (Klein et al. 2021, see also Suárez Mascareño et al. 2016); thus, this prior seems appropriate, and likely conservative.
The same priors are used for the orbital parameters of all the N_{p} Keplerians in the model. We assign a wide loguniform prior to the orbital periods, between 1 day and the full time span of observations (815 days), and a modified loguniform prior for the semiamplitudes, up to 10 m s^{−1} with a break point at 1 m s^{−1}. The eccentricities were assigned a Kumaraswamy prior (Kumaraswamy 1980), which closely resembles the Beta distribution proposed by Kipping (2013) but is easier to manipulate numerically. The two angular parameters, M_{0} and ω, were assigned uniform priors between 0 and 2π. Other parameters such as offsets and jitters were assigned datadependent but uninformative priors. The full set of priors is listed in Table B.1 and discussed further in Appendix B.
To sample from the posterior distribution, we use the diffusive nested sampling (DNS) algorithm from Brewer et al. (2011), as implemented in kima (Faria et al. 2018). Together with posterior samples, DNS provides an estimate for the marginal likelihood, or evidence, of the model, which we can use for model comparison (e.g. Brewer 2014; Feroz et al. 2011). We obtain at least 50 000 effective samples from the posterior – as measured by the effective sample size (ESS), the number of samples with significant posterior weight – for each model, which is more than enough to accurately characterise it.
Fig. 1 Comparison between RVs derived with the CCF and TM techniques. Top panels: two RV time series, and the middleleft and middleright panels show the generalised LombScargle (GLS) periodograms of CCF and TM RVs, respectively.In the bottom panels we display the histograms of the RVs (left) and of the RV uncertainties (right) for the two time series. 
Fig. 2 Prewhitening procedure applied to the CCF RVs. Top panel: observed RVs together with the GP prediction. The periodogram of the residuals from this fit is shown in panel a, and the periodograms after two successive sinusoidal fits are in panels b and c. The false alarm probability of 1%, calculated with bootstrap randomisation, is shown by the horizontal grey lines. 
4.3 Results from the joint model
The evidence values for models with N_{p} = 0, 1, 2, 3 are shown inTable 2, together with the log Bayes factors between consecutive models. In the analysis of the CCF RVs, the model with N_{p} = 2 is significantly preferred, with a ΔlnZ > 5, which corresponds to decisive evidence in the scale of Kass & Raftery (1995).
From this model, the posterior distribution for the orbital periods is shown in the top panel of Fig. 3, where the count ineach period bin was normalised by the ESS (thus resembling the TIP proposed by Hara et al. 2022). The posterior shows two very clear peaks at 5.12 and 11.19 days. There is residual posterior probability close to 1.2 days, which is a 1day alias of the 5.12day period and an even smaller posterior peak at 43 days (both are not visible at the scale of the figure). The two main posterior peaks correspond to the orbital periods of Proxima b and of a candidate planet that we call Proxima d.
While the period and semiamplitude of the 5day signal are well constrained, the posterior for the eccentricity is quite wide, with a peak at 0.45 (see Fig. C.2). Above this value, it also shows a sharp upper tail due to the AMD stability criterion, which makes higher eccentricities very improbable. The median and 68% credible intervals result in an estimate for the eccentricity of .
In this analysis of the CCF RVs, the 5day signal has a semiamplitude of 59 cm s^{−1} (posterior median), leading to a minimum mass of only 0.36 ± 0.06 M_{⊕}, and a semimajor axis of 0.028 au, assuming it is due to a planetary companion. For both the 5day and the 11day signals, circular orbits are not excluded at the 2σ level, which led us to consider a model where the eccentricities of both planets are fixed to zero. The evidence for this model (with both CCF and TM RVs) is also shown in Table 2, and it is in both cases slightly higher than for the model with free eccentricities, although not significantly. This means that the currently available ESPRESSO data are not sufficient to distinguish between eccentric and circular orbits at the twosigma level.
The analysis of the TM RVs provides qualitatively similar results, although the Bayes factor between the models with N_{p} = 1 and N_{p} = 2, even if still decisive, is not as high (Table 2). The same two clear peaks at 5.12 days and 11.19 days are visible in the posterior forthe orbital periods (shown in the bottom panel of Fig. 3), but the alias at 1.2 days shows a higher posterior probability.
In terms of likelihood alone, the solutions with the two orbital periods of 5.12 days and 11.19 days are well isolated, with relative to any other combination of periods (except for the 1.2day alias). The same is true for the analysis of the CCF RVs. Despite having comparable periodogram power in the analysis from Sect. 4.1 (panel b in Fig. 2), the 5.12day period shows a much larger posterior probability and likelihood when compared to the 1.2day alias, in both RV datasets, leading to the conclusion that the former is the correct periodicity.
In Fig. 4, we show the joint posterior distributions for the orbital period, semiamplitude, eccentricity, and planet minimum mass of the two Keplerians, using the stellar mass in Table 1. In the TM dataset, the semiamplitude of the 5day signal is significantly smaller than for the CCF RVs, at 38 cm s^{−1} (Fig. 4), although the results for the two datasets are compatible at less than 2σ. This ultimately leads to an estimate for the planetary mass that is 30% lower than with the CCF RVs, at 0.24 ± 0.05 M_{⊕} (posterior median). The eccentricity is better constrained by the TM RVs and the solutions for both planets are compatible with circular orbits. We note that, for such a low semiamplitude, the departure from a sinusoidal signal caused by an eccentricity of, say, 0.1, would be of the order of K ⋅ e ≈ 5 cm s^{−1}, which is below the RV precision we can currently probe.
Figure 5 shows the ESPRESSO RVs phasefolded to the orbital periods of the two Keplerians, assuming the maximum a posteriori (MAP) solution from the twoplanet model applied to the TM RVs. We finally adopt this solution and the 68% quantiles of the posterior as the estimates for the parameters of the twoplanet model, listed in Table C.1. We note that the choice for quoting the results from the TM RVs over the CCF RVs is mostly guided by the substantially smaller RV uncertainties on this dataset.
In Fig. C.1, we show the individual fits to the ESPRESSO18, ESPRESSO19, and ESPRESSO21 subsets as well as the combined residuals. The weighted RMS of the RV residuals reaches 27 and 60 cm s^{−1} for the FWHM. We note a smaller scatter in the ESPRESSO19 residuals (in both RV and FWHM), which could be explained by the fewer observations in this subset or the larger time separation between the points. Indeed, the ESPRESSO19 subset does not contain nights with more than one observation, which may suggest some degree of overfitting and an inability of the model to represent intranight variations. In other words, the model may not completely capture all the RV (and FWHM) variations at short periods of 1 or 2 days.
The stellar activity component of the model (the shared GP for RVs and FWHM), is very well constrained both in the CCF and TM datasets. The posteriors for the GP hyperparameters from the analysis of the TM RVs are shown in Fig. 6. In the RVs, the contribution from the GP reaches 2.3 m s^{−1}, well above the amplitudes of either of the two planet signals. In the FWHM, the GP’s contribution is 9.4 m s^{−1}.
The rotation period of the star, parameterised by η_{3}, is constrained at a value of 85.1 days (posterior median), with a 1σ uncertainty below 1 day. This value is close to half of the estimated timescale of evolution, η_{2} = 151 days, suggesting that the active regions at the surface of the star (or at least the activity signal itself) evolve over a timescale of two rotation periods, similarly to what is observed in the Sun. The relatively large uncertainty on η_{2} could be related to our assumption of sharing this parameter between RVs and FWHM, which might not be strictly valid (and indeed the same is true for η_{4}). We discuss briefly this assumption in Appendix A.
Finally, we note that the orbital period prior used in these analyses does not allow for the detection of Proxima c, which has an orbital period of 1900 days (Damasso et al. 2020). If the signal from this planet is present in the ESPRESSO dataset, we would expect it to be absorbed by the quadratic RV trend considered in our models. Nevertheless, for both the CCF and TM RVs, the coefficients of the quadratic trend are found to be compatible with zero. The limited time span, together with the inclusion of RV offsets between the three subsets of data makes the detection of Proxima c using only the currently available ESPRESSO data very challenging, in line with prediction from Damasso & Del Sordo (2020). A complete analysis of the ESPRESSO data together with other RV datasets could shed light on the presence of this planet, but may require a more general activity model that includes longterm activity variations such as the magnetic cycle detected in Proxima (see e.g. Suárez Mascareño et al. 2016; Wargelin et al. 2017).
Fig. 3 Posterior distribution for the orbital periods in the twoplanet model from the analysis of the CCF RVs (top) and TM RVs (bottom). The posteriors are shown normalised by the ESS per histogram bin (so the maximum possible value in the abscissa is 1). The prior is loguniform from 1 day to the time span of the data, which is marked with a dashed line. 
Evidence (ln Z) and Bayes factors (Δ ln Z) for models with a given number of Keplerians, N_{p}, from the analysis of the CCF and TM RVs.
Fig. 4 Joint and marginal posteriors for the orbital periods, semiamplitudes, eccentricities, and minimum planet masses, M_{p} sin i, in Earth masses, of the 11day (left) and 5day (right) Keplerian signals. These results come from the analysis of the TM RVs. The estimates at the top of each panel correspond to the median and 68% quantiles of the posteriors. 
Fig. 5 Templatematching RVs from ESPRESSO phasefolded to the orbital periods of Proxima b (left) and Proxima d (right). TheGP and background components of the model were subtracted assuming the MAP solution. The full fit to the RVs andFWHM is shown in Fig. C.1. 
Fig. 6 Prior and posterior distributions for the hyperparameters of the GP in the twoplanet model. Posterior estimates (median and 68% credible intervals) are show in each panel. It should be noted that for η_{3} and η_{4} the full support of the prior is larger than what is shown. 
5 Assessing the nature of the 5day signal
The RV analysis in the previous section included a GP component to model activityinduced RV variations and used the FWHM of the CCF as an activity indicator. Nevertheless, we cannot discard the possibility that the 5day signal originates from an incomplete model of either stellar activity or instrumental systematics, or both. In this section we study the evolution of the 5day signal as the number of ESPRESSO RVs increased, look at other activity indicators to confirm the planetary nature of the signal and study the chromatic behaviour of the RVs.
5.1 Evolution with number of points
A planetary signal gives rise to strictly periodic RV variations and its stationary nature implies that it should become more significant as the number of observations increases. Using the same model and priors as in Sect. 4, we analysed the datasets built from the first n ESPRESSO observations, with n = 30, 40…, 100, 114. We note that this means the upper limit of the prior for the orbital periods increases with n, but this is of little consequence for this analysis given that the periodicities of interest are well covered already on the shortest dataset.
The results are shown in Fig. 7, which displays the evolution of the Bayes factor between the models with N_{p} = 1 and N_{p} = 2 together with the posterior distributions for the orbital periods and semiamplitudes in the twoplanet model, for increasing values of n. The Bayes factor increases monotonically until it reaches the value obtained for the full dataset. We note the large increase after the initial observing campaign presented in SM2020 and, after about 90 observations, when we adapted the observing strategy to mitigate aliasing by obtaining two RV points per night when possible.
It is clear that the posterior probability of the 5day signal increases with the number of RVs in the dataset. The 11day signal shows the same behaviour, even though its posterior is already well defined with fewer data points, due to the larger amplitude. At the same time, the semiamplitudes of both signals also become better constrained as n increases, as seen in the bottom panel of Fig. 7. These results provide evidence that the 5day signal is stable over time, as expected fora planetary origin.
5.2 Periodicities in activity indicators
Besides the FWHM, a number of other activity indicators are available from the CCF or from individual spectral lines and can be used to identify activity signals. Moreover, all the indicators are obtained simultaneously to the RVs, so that the imprint of the observational sampling, if present, should appear similarly in each time series.
Since any activityrelated signal that might be present in the indicators is likely not periodic, the typical LombScargle periodogram might not be the ideal statistic to measure them. An arguably better model could use a QP GP to try to identify the stellar rotation signal and possibly an additional periodic component. But even this model might not be adequate, especially if it assumes the same priors as for the RVs, since the indicators can be sensitive to activity in different ways than the RVs, or be affected by instrumental effects.
In addition, Proxima is known for its frequent flares (e.g. Davenport et al. 2016), which introduces relatively strong outliers in the measurements of some activitysensitive spectral lines, such as the Ca H&K, Hα or Na lines. These are indeed true outliers, in the sense that the estimated uncertainty of the affected points is not necessarily higher than that ofunaffected observations.
For these reasons, we decided to use a simple surrogate model to study the periodicities present in the activity indicators. We consider up to five sinusoidal signals, with periods constrained between 1 day and the time span of the observations. A Student’s t likelihood is used instead of the typical Gaussian in order to make the analysis more robust to outliers (the degrees of freedom of the distribution is a free parameter). We then plot the normalised posterior distribution for the five periods,which conveys information about the significance of each periodicity in the different time series. In practice, this results in a ‘cleaned’ and robust version of the periodogram that retains only the most significant periods.
The results are shown in Fig. 8 for the FWHM, CCF bisector, CCF contrast, , Na index, and Hα index. Notably, the FWHM,, and CCF contrast show significant periodicities close to the stellar rotation period and its first harmonic and all indicators show more or less significant periodicities at the second harmonic, close to 28 days. No clear peaks are seen around 11 or 5 days, suggesting that the two signals detected in the RVs are not related to activity, at least to the extent it is captured by the indicators. Again, these results add to the conclusion that the two signals are better explained as due to planetary companions.
Fig. 7 Evolution of the Bayes factor between models with one and two Keplerians (top) and of the posterior distribution for the orbital periods (middle) and semiamplitudes (bottom) with an increasing number of ESPRESSO observations. The dashed vertical line corresponds to the number of observations available to SM2020. 
Fig. 8 Posterior distributions for the periods of sinusoids in the analysis of activity indicators. The distributions are normalised and shown on the same vertical scale, which allows for a visual comparison of the significance of the periodicities. For example, the CCF contrast shows more significant periodicities (close to 100 and 40 days) than the Hα index. The stellar rotation period and its first and second harmonics are denoted with filled lines, and the 11day and 5day signals detected in the RVs with dashed red lines. 
5.3 Chromatic RV variations
As shown in SM2020, the wide wavelength coverage of ESPRESSO, together with the collecting power of the VLT, allow for the derivation of socalled chromatic RVs by dividing the spectra into a few wavelength bins. SM2020 divided the ESPRESSO spectra into blue (440–570 nm), green (570–690 nm), and red (730–790 nm) bins, chosen to guarantee a similar RV precision. We select the same spectral orders to calculate chromatic RVs within the TM approach.
The advantage of calculating chromatic RVs with the TM approach comes from the improved RV precision obtained when using this technique (see Fig. 1). The blue, green, and red chromatic TM RVs show median RV uncertainties of 29, 24, and 22 cm s^{−1}, respectively. In contrast, using the CCFs (and the same spectral orders) we achieve median uncertainties of 47, 48, and 48 cm s^{−1} for the three wavelength regions. The latter are comparable to the semiamplitude measured for the 5day signal, making a chromatic analysis of this signal challenging with the CCF RVs but possible with the TM RVs.
We analysed the blue, green, and red TM RV datasets individually, with a similar model to that of Sect. 4. The only difference is that here we assume only circular orbits, to minimise the number of parameters in the model and thus decrease the sensitivity to the larger contribution from photon noise to the error budget. The results are shown in Fig. 9. The top and middle panels show the phase curves of the 11day and 5day signals, respectively, while the bottom panels show the posterior distributions for η_{1} RV and η_{1} FWHM.
In all three wavelength bins the estimated semiamplitudes are similar, varying by 16 and 5 cm s^{−1} for the 11 and 5day signals, respectively (about 13% in both cases). A variation with wavelength would have been expected from a signal that was caused by a stellar spot (e.g. Desort et al. 2007; Figueira et al. 2010). Nevertheless, the posterior probability for 5day signal is lower in the red RVs, resulting in a less significant detection in this wavelength bin. This could be explained by a larger telluric contamination or a more complicated activity signal.
Indeed we find hints that the amplitude of the activity signal in the RVs is larger in the green and red bins than in the blue, at least as captured by the GP model (see bottom panels in Fig. 9). However, the posterior estimates are mostly compatible, so this result should not be overinterpreted. In summary, this analysis of chromatic RVs provides us with additional confidence that the 5day signal is of planetary origin.
Fig. 9 Analysis of chromatic TM RVs. Three columns: results for the blue, green, and red regions, corresponding to the wavelengths between 440–570 nm, 570–690 nm, and 730–790 nm, respectively. Top and middle panels: ESPRESSO RVs phasefolded at the 11 and 5day periods, together with the maximum likelihood solution for each signal. Bottom panels: posteriors for the GP amplitude, η_{1} RV, with the posterior median and 68% quantiles as labels. 
6 Discussion and conclusions
In this work we have analysed a set of 114 ESPRESSO observations of Proxima. The imprint of Proxima b is the clearest signal seen in the ESPRESSO RVs (as in SM2020) and is detected with very high significance in both the CCF and TM datasets. The orbital parameters are consistent in the two analyses, even if the TM RVs suggest a slightly smaller semiamplitude. This smaller amplitude is closer to the one found by SM2020 from the combined analysis of ESPRESSO, HARPS, and UVES data. Using the TM RVs, the minimum mass of this planet is estimated as 1.07± 0.06M_{⊙}, with a relative uncertainty below 6% (see Table C.1) when taking the uncertainty in the stellar mass into account.
A significant 5.12day signal is detected in the ESPRESSO data, which we attribute to a planet candidate, Proxima d, first found in SM2020. We detect the signal with a slightly higher significance on the CCF RVs than on the TM RVs (Table 2). It is unlikely that this signal is directly connected to the stellar variability at the rotation period, which we estimate to be 84.5 days. First, the Keplerian signal becomes more significant as ESPRESSO observations are added, with its amplitude and period becoming increasingly better constrained, in a similar way to those of the 11day signal. Second, none of the activity indicators we analysed show periodicities at either 5 or 11 days, though they do vary with the rotation period and its harmonics. Finally, the 5day signal shows a wavelengthindependent amplitude within the range observed by ESPRESSO. These tests show evidence for the stationary and achromatic nature of the signal and are thus strong arguments in favour of the planetary nature of Proxima d.
The estimated amplitude and eccentricity of Proxima d are smaller on the TM RVs and lead to an estimate for the minimum mass of 0.26 ± 0.05 M_{⊕}. This planet orbits closer to the star (at 0.02885 au) than to the inner edge of the HZ (see Table 1), and it is the lightest planet detected with the RV method so far.
Together with other recent detections (Demangeon et al. 2021; LilloBox et al. 2021), our results showcase the extreme RV precision systematically achieved with ESPRESSO. Even in the presence of stellar activity signals causing RV variations of the order of m s^{−1}, it is now possible to detect and measure precise masses for very lowmass planets that induce RV signals of only a few tens of cm s^{−1}.
Our results bring new interest to the planetary system around Proxima. The star is now the likely host of three lowmass planets. The habitability conditions of Proxima b, which orbits within the HZ of the star, have been extensively studied (e.g. Barnes et al. 2017; Ribas et al. 2016; Turbet et al. 2016; Meadows & Barnes 2018). On the other hand, the candidate Proxima d orbits much closer to the star and outside the HZ range.
The atmospheric properties of these planets may be affected by stellar magnetic activity (e.g. Garraffo et al. 2016; GarciaSage et al. 2017). Recently, Klein et al. (2021) studied the extended magnetosphere of Proxima using circular polarisation spectra measured with HARPSPol and estimated the typical radius of the spherical Alfvén surface to be about 25 R_{*}. They concluded that Proxima b, which orbits at a distance of more than 70 R_{*}, should lie in the superAlfvénic regime, with no direct star–planet magnetic connection. This conclusion was also supported by Kavanagh et al. (2021). The inner planet, planet d, orbits at about 40 R_{*}, still in the superAlfvénic regime.
With the small stellar radius, Proxima d has a transit probability just above 2%. Its equilibrium temperature may reach 360 K, assuming a Bond albedo of 0.3 (e.g. Seager et al. 2010). From the planetary properties and stellar parameters, we can estimate a planetary radius of 0.81 ± 0.08 R_{⊙} using the random forest model from UlmerMoll et al. (2019), leading to a transit depth of about 0.3% (approximately half that of Proxima b; AngladaEscudé et al. 2016). A transit detection would allow a precise measurement of the planetary radius and could place constraints on the planet bulk density and possible atmosphere. However, a transit is unlikely given that Proxima b has not been found to transit (see Jenkins et al. 2019, and references therein) and that transit events at periods below 5 days and depths above 3 mmag have been ruled out (Feliz et al. 2019; Vida et al. 2019).
While the currently available ESPRESSO data do not allow for the detection of Proxima c, a full characterisation of the system may be possible with continued spectroscopic observations or with Gaia astrometry (Damasso & Del Sordo 2020; Damasso et al. 2020). A complete analysis of the available UVES, HARPS, and ESPRESSO datasets (see Damasso et al. 2020; AngladaEscudé et al. 2016, and references therein), perhaps after a homogeneous RV derivation with TM, would further constrain the orbits of all the planets in the system. Given that this combined dataset would span several years, a more comprehensive stellar activity model (maybe including a longterm magnetic cycle) would have to be used, still extracting information from appropriate activity indicators. In summary, further observations of Proxima will help to complete our understanding of this multiplanetary system and may even unveil the presence of additional planets.
Acknowledgements
The authors acknowledge the ESPRESSO project team for its effort and dedication in building the ESPRESSO instrument. This work was supported by FCT – Fundação para a Ciência e Tecnologia – through national funds and by FEDER through COMPETE2020 – Programa Operacional Competitividade e Internacionalização by these grants: UID/FIS/04434/2019; UIDB/04434/2020; UIDP/04434/2020; PTDC/FISAST/32113/2017 and POCI010145FEDER032113; PTDC/FISAST/28953/2017 and POCI010145FEDER028953; PTDC/FISAST/28987/2017 and POCI010145FEDER028987; PTDC/FISAST/30389/2017 and POCI010145FEDER030389. J.P.F. and O.D. are supported in the form of work contracts funded by national funds through FCT with the references DL 57/2016/CP1364/CT0005 and DL 57/2016/CP1364/CT0004. A.M.S. acknowledges support from FCT through the Fellowship 2020.05387.BD and POCH/FSE(EC). A.S.M. acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) under 2018 Juan de la Cierva program IJC2018035229I. A.S.M., J.I.G.H., and R.R. acknowledge financial support from the MICINN project PID2020117493GBI00 and from the Government of the Canary Islands project ProID2020010129. J.I.G.H. also acknowledges financial support from the Spanish MICINN under 2013 Ramón y Cajal program RYC201314875. Y.A. acknowledges the support of the Swiss National Fund under grant 200020_172746. The INAF authors acknowledge financial support of the Italian Ministry of Education, University, and Research with PRIN 201278X4FL and the “Progetti Premiali” funding scheme. J.LB. acknowledges financial support received from “la Caixa” Foundation (ID 100010434) and from the European Union’s Horizon 2020 research and innovation programme under the Marie SkłodowskaCurie grant agreement No 847648, with fellowship code LCF/BQ/PI20/11760023. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. This project has received funding from the European Research Council (ERC) under the following European Union’s Horizon 2020 research and innovation programmes: grant agreement SCORE No 851555 and project FOUR ACES; grant agreement No 724427. D.E. acknowledges financial support from the Swiss National Science Foundation for project 200021_200726. V.A. acknowledges the support from FCT through Investigador FCT contract nr. IF/00650/2015/CP1273/CT0001. N.J.N. acknowledges support form the following projects: UIDB/04434/2020 and UIDP/04434/2020, CERN/FISPAR/0037/2019, PTDC/FISOUT/29048/2017, COMPETE2020: POCI010145FEDER028987 and FCT: PTDC/FISAST/28987/2017. R.A. is a Trottier Postdoctoral Fellow and acknowledges support from the Trottier Family Foundation. This work was supported in part through a grant from FRQNT. This research has made use of the SIMBAD database, the VizieR catalogue operated at CDS, France, and the DACE platform developed in the frame of PlanetS (https://dace.unige.ch).
Appendix A Analysis with different GP covariances
One important choice when modelling the stellar activity signal with a GP is the kernel used to build the covariance matrix. Throughout the paper, we discussed the results obtained with the QP kernel, which is given by (A.1)
in terms of the timelag τ = t_{i} − t_{j}. This covariance is still the most commonly used in the analysis of RV observations (e.g. Grunblatt et al. 2015; Haywood et al. 2014). Recently, Perger et al. (2021) suggested a modification of the QP kernel to ac count for an explicit periodic component at η_{3}∕2, which is typically observed in simulated data. This QPC kernel introduces one additional hyperparameter, η_{5}, which controls the amplitude of the cosine term: (A.2)
We performed the same analysis as in Sect. 4, replacing the QP kernel with the QPC and assuming the same priors (see Ap pendix B). We found virtually identical results to those obtained with the QP kernel, as seen in Fig. A.1, which compares the pos terior distributions for the kernel parameters.
Fig. A.1 Posterior distributions for the hyperparameters of the GP kernels from the analysis of the TM RVs. We quote the posterior median and 68% quantiles for each parameter and kernel. 
Another assumption we made when jointly modelling RVs and FWHM was to share the η_{2}, η_{3}, and η_{4} hyperparameters be tween the two time series. This is motivated by the expectation that the signal in RV and FWHM will share the same correlation structure. However, stellar activity could appear in the two ob servables with different characteristics. We nevertheless expect η_{3} to correspond to the stellar rotation period and thus be the same in RV and FWHM. To test these assumptions, we repeated the analyses considering independentη_{2} and η_{4} parameters for the RVs and the FWHM (with the same priors). For both the CCF and TM RVs, this model revealed no significant differences from the original model with shared parameters.
Fig. A.2 Predictive mean of the GP component of the model versus ro tation phase, assuming a value of η_{3} = 85.1 days. The top and bottom panels show the predictions for the RVs and the FWHM, respectively. The curves are only shown for times in which there are ESPRESSO ob servations, and their colour reflects the time since the first observation. 
A.1 Differential rotation and activity evolution
Wargelin et al. (2017) found evidence for differential rotation (DR) on Proxima by studying Vband observations from the ASAS survey. The stellar rotation period changed between 77 and 90 days over a period of about 10 years, corresponding to a fractional DR estimate of ΔP_{rot}∕⟨P_{rot}⟩≃ 0.16.
Since the ESPRESSO dataset spans almost ten rotation periods and the three subsets each span about one rotation period, it is interesting to check if there are hints of DR showing in these data. While the GP model for stellar activity does not explicitly account for DR, it can in principle model its effects on both the RVs and the FWHM.
We use the maximum likelihood solution from the model with two planets fit to the TM RVs. Figure A.2 shows the predictive mean for the GP component (after subtracting the two Keplerian signals and the quadratic trend), phasefolded on the value of η_{3}, with the colour of the curves corresponding to the time since the first observation. The predictive is only shown for those times in which there are observations (that is, excluding the gap between ESPRESSO19 and ESPRESSO21) because, in the absence of data, the GP tends to the mean of zero and its interpretation is less meaningful.
At a glance, Fig. A.2 highlights the quasiperiodicity of the activity signal, as modelled by the GP. Moreover, as is also visible in Fig. C.1 below, there is a clear time delay between the activity signal on the RVs and on the FWHM, with the maxima (and minima) of the FWHM happening 1015 days, or about 15% of the rotation period, after the RV maxima. This temporal shift has been observed before for the Sun (Collier Cameron et al. 2019)and other stars (e.g. Queloz et al. 2009; Santos et al. 2014) and is related to the effect of active regions on the stellar line profile (see Collier Cameron et al. 2019, and references therein).
Perhaps more interestingly, Fig. A.2 also shows a small variation in the rotation period, as measured by the position of the maxima and minima of the activity signal, over short timescales of two to three rotation periods (c.f. the blue and purple minima seen on the left in the top and bottom panels). Assuming that the same active regions are creating the signal over these timescales, such a variation could be due to DR as the active regions move differentially on the stellar surface.
Appendix B Prior distributions
Here we detail the prior distributions and further justify some of the choices made when defining them. For clarity, we refer to three distinct components of the model: background, Keplerians, and noise.
The background model includes the systemic velocity and the quadratic trend coefficients in the RVs, as wellas the ‘systemic FWHM’, which is a simple mean value for this quantity. Our priors for v_{sys} and f_{sys} are slightly unusual, but simply limit the parameters to the observed ranges and do not influence the results, especially when taking into account that the orbital pe riod of any Keplerian component is limited to the time span of the data (see below). For the RV and FWHM offsets between the (three) subsets of ESPRESSO data, we chose uniform priors within the observed ranges. The Gaussian priors for the quadratic trend coefficients are, again, informed by the data but only to the extent of setting typical scales for these parameters.
Our rationale when setting priors for the orbital parameters was to use the same uninformative distributions for all N_{p} Kep lerians included in the model. In principle, there is prior information about Proxima b that is independent of the ESPRESSO data, and so could be used. In practice, the signal is so clearly detected in the ESPRESSO data that this prior becomes less relevant. On the other hand, there is no prior information about the 5day signal that can be used since its detection relies solely on ESPRESSO data. A loguniform prior limited at the time span of the data is thus a reasonable choice.
The most restrictive prior ends up being the one for K, which we limit to 10 m s^{−1}. This is simply because there is no indication of higheramplitude signals being present in the ESPRESSO data, and a prior that used a simple statistic of the observed RVs (e.g. the RV span) would be virtually the same. More than 99.9% of the posterior probability ends up being below 1.6 m s^{−1}, suggesting that this choice of upper limit does not influence the results. The Kumaraswamy prior for the orbital eccentricities and its parameters are motivated by Kipping (2013) where a Beta distribution was proposed. The Kumaraswamy is virtually identical to the Beta (Kumaraswamy 1980), but its cumulative distribution function (required by DNS) is analytical.
In the noise component, we include the GP and the individual jitters for each subset of data. The hyperparameters of the QP kernel are η_{1−4} (see Eq. A.1) and we consider independent η_{1} parameters for the RVs and the FWHM. These parameters are assigned modified loguniform priors up to the observed range and which extend to zero. The justification for the prior for η_{3} was mentioned before, with the wide uniform prior we use being quite conservative. Since we expect the activity signal to be QP, the prior for η_{2} starts at the lower limit for η_{3}, and extends to a reasonable limit of 400 days, almost five times the rotation period of Proxima. We chose a loguniform prior due to the large span. We also tested a different prior for η_{2} that extended down to zero (an inverse Gamma distribution) but this did not change the results.
Parameters and prior distributions of our model for the RV and FWHM data.
Appendix C Posterior estimates
Table C.1 lists the estimates for the parameters of the twoplanet model resulting from the analysis of the CCF and TM RVs. Each estimate is shown as the (MAP) value of the parameter together with the 68% quantiles of the posterior distribution.In Fig. C.1, we show the fit to the ESPRESSO TM RVs and FWHM, individually for each subset of the data together with the residuals after subtracting the complete model. These panels also show the weighted RMS of the residuals assuming the original RV and FWHM uncertainties (i.e. the jitter parameters have not been added in quadrature).
Finally, Fig. C.2 shows the joint posteriors for the orbital pe riod, semiamplitude, eccentricity, and planet minimum mass of the two Keplerians from the analysis of the CCF RVs. In this dataset, the semiamplitude of the 5day signal is significantly smaller than for the CCF RVs, at 38 cm s^{−1} (Fig. 4), although the results for the two datasets are compatible at less than 2σ. This ultimately leads to an estimate for the planetary mass that is 30% lower than with the CCF RVs, at 0.24± 0.05 M_{⊕}. The eccentricity is better constrained by the TM RVs and the solutions for both planets are compatible with circular orbits.
Fig. C.1 Maximum a posteriori solution for the twoplanet model on the TM RVs. The top panels show the RV observations for ESPRESSO18, ESPRESSO19, and ESPRESSO21, together with the GP component of the model (in pink) as well as the full model (in black). The RV residuals are shown just below, highlighting the full residual rms. The two bottom panels show the same for the FWHM. 
Fig. C.2 Joint and marginal posteriors for the orbital periods, semiamplitudes, eccentricities, and minimum planet masses, M_{p} sin i, in Earth masses, of the 11day (left) and 5day (right) Keplerian signals. These results come from the analysis of the CCF RVs. The estimates at the top of each panel correspond to the median and 68% quantiles of the posteriors. 
Posterior estimates for the parameters of the twoplanet model from the analysis of CCF and TM RVs. For each parameter we show the MAP estimate and the 68% quantiles of the distribution.
References
 AngladaEscudé, G., & Butler, R. P. 2012, ApJS, 200, 15 [Google Scholar]
 AngladaEscudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437 [Google Scholar]
 Artigau, É., Kouach, D., Donati, J.F., et al. 2014, SPIE, 9147, 914715 [NASA ADS] [Google Scholar]
 Baglin, A., Auvergne, M., Barge, P., et al. 2006, ESA SP, 1306, 33 [NASA ADS] [Google Scholar]
 Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barnes, J. R., Jeffers, S. V., AngladaEscudé, G., et al. 2017, MNRAS, 466, 1733 [NASA ADS] [CrossRef] [Google Scholar]
 Benedict, G. F., & McArthur, B. E. 2020, Res. Notes Am. Astron. Soc., 4, 86 [Google Scholar]
 Bertaux, J. L., Lallement, R., Ferron, S., Boonne, C., & Bodichon, R. 2014, A&A, 564, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bessell, M. S. 1991, AJ, 101, 662 [Google Scholar]
 Bleistein, N., & Handelsman, R. A. 1986, Asymptotic Expansions of Integrals (New York: Dover Publications) [Google Scholar]
 Bonfils, X., AstudilloDefru, N., Díaz, R., et al. 2018, A&A, 613, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Borucki, W. J. 2016, Rep. Prog. Phys., 79, 036901 [Google Scholar]
 Bouchy, F., Pepe, F., & Queloz, D. 2001, A&A, 374, 733 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boyajian, T. S., von Braun, K., van Belle, G., et al. 2012, ApJ, 757, 112 [Google Scholar]
 Brewer, B. J. 2014, ArXiv eprints [arXiv:1411.3921] [Google Scholar]
 Brewer, B. J., Pártay, L. B., & Csányi, G. 2011, Stat. Comput., 21, 649 [Google Scholar]
 Broeg, C., Fortier, A., Ehrenreich, D., et al. 2013, EPJ Web Conf., 47, 03005 [Google Scholar]
 Collier Cameron, A., Mortier, A., Phillips, D., et al. 2019, MNRAS, 487, 1082 [Google Scholar]
 Cosentino, R., Lovis, C., Pepe, F., et al. 2012, Proc. SPIE, 8446, 84461v [Google Scholar]
 Damasso, M., & Del Sordo, F. 2020, MNRAS, 494, 1387 [CrossRef] [Google Scholar]
 Damasso, M., Del Sordo, F., AngladaEscudé, G., et al. 2020, Sci. Adv., 6, eaax7467 [Google Scholar]
 Davenport, J. R. A., Kipping, D. M., Sasselov, D., Matthews, J. M., & Cameron, C. 2016, ApJ, 829, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Delfosse, X., Forveille, T., Ségransan, D., et al. 2000, A&A, 364, 217 [NASA ADS] [Google Scholar]
 Demangeon, O. D. S., Osorio, M. R. Z., Alibert, Y., et al. 2021, A&A, 653, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Desort, M., Lagrange, A.M., Galland, F., Udry, S., & Mayor, M. 2007, A&A, 473, 983 [CrossRef] [EDP Sciences] [Google Scholar]
 Díaz, R. F., Delfosse, X., Hobson, M. J., et al. 2019, A&A, 625, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ehrenreich, D., Lovis, C., Allart, R., et al. 2020, Nature, 580, 597 [Google Scholar]
 Faria, J. P., Santos, N. C., Figueira, P., & Brewer, B. J. 2018, JOSS, 3, 487 [NASA ADS] [CrossRef] [Google Scholar]
 Feliz, D. L., Blank, D. L., Collins, K. A., et al. 2019, AJ, 157, 226 [NASA ADS] [CrossRef] [Google Scholar]
 Feroz, F., Balan, S. T., & Hobson, M. P. 2011, MNRAS, 415, 3462 [Google Scholar]
 Figueira, P., Marmier, M., Bonfils, X., et al. 2010, A&A, 513, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fischer, D. A., AngladaEscude, G., Arriagada, P., et al. 2016, PASP, 128, 066001 [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2016, A&A, 595, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 GarciaSage, K., Glocer, A., Drake, J. J., Gronoff, G., & Cohen, O. 2017, ApJ, 844, L13 [NASA ADS] [CrossRef] [Google Scholar]
 Garraffo, C., Drake, J. J., & Cohen, O. 2016, ApJ, 833, L4 [NASA ADS] [CrossRef] [Google Scholar]
 Gomes da Silva, J., Figueira, P., Santos, N. C., & Faria, J. P. 2018, JOSS, 3, 667 [NASA ADS] [CrossRef] [Google Scholar]
 Gratton, R., Zurlo, A., Le Coroller, H., et al. 2020, A&A, 638, A120 [EDP Sciences] [Google Scholar]
 Gregory, P. C. 2005, ApJ, 631, 1198 [NASA ADS] [CrossRef] [Google Scholar]
 Grunblatt, S. K., Howard, A. W., & Haywood, R. D. 2015, ApJ, 808, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Hara, N. C., Unger, N., Delisle, J.B., Díaz, R., & Ségransan, D. 2022, A&A, in press https://doi.org/10.1051/00046361/202140543 [Google Scholar]
 Haywood, R. D., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517 [Google Scholar]
 Hojjatpanah, S., Figueira, P., Santos, N. C., et al. 2019, A&A, 629, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jao, W.C., Henry, T. J., Subasavage, J. P., et al. 2014, AJ, 147, 21 [Google Scholar]
 Jeffers, S. V., Dreizler, S., Barnes, J. R., et al. 2020, Science, 368, 1477 [Google Scholar]
 Jenkins, J. S., Harrington, J., Challener, R. C., et al. 2019, MNRAS, 487, 268 [Google Scholar]
 Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773 [Google Scholar]
 Kavanagh, R. D., Vidotto, A. A., Klein, B., et al. 2021, MNRAS, 504, 1511 [Google Scholar]
 Kervella, P., Arenou, F., & Schneider, J. 2020, A&A, 635, L14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kipping, D. M. 2013, MNRAS, 434, L51 [Google Scholar]
 Klein, B., Donati, J.F., Hébrard, É. M., et al. 2021, MNRAS, 500, 1844 [Google Scholar]
 Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Kumaraswamy, P. 1980, J. Hydrol., 46, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Lafarga, M., Ribas, I., Lovis, C., et al. 2020, A&A, 636, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J., & Petit, A. C. 2017, A&A, 605, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LilloBox, J., Figueira, P., Leleu, A., et al. 2020, A&A, 642, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LilloBox, J., Faria, J. P., Mascareño, A. S., et al. 2021, A&A, 654, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lomb, N. R. 1976, Astrophys. Space Sci., 39, 447 [Google Scholar]
 Mahadevan, S., Ramsey, L., Bender, C., et al. 2012, Proc. SPIE, 8446, 84461S [Google Scholar]
 Mann, A. W., Feiden, G. A., Gaidos, E., Boyajian, T., & von Braun, K. 2015, ApJ, 804, 64 [Google Scholar]
 Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20 [NASA ADS] [Google Scholar]
 Meadows, V. S., & Barnes, R. K. 2018, in Handbook of Exoplanets, eds. H. J. Deeg & J. A. Belmonte (Cham: Springer International Publishing), 2771 [CrossRef] [Google Scholar]
 Pavlenko, Y., Suárez Mascareño, A., Rebolo, R., et al. 2017, A&A, 606, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perger, M., AngladaEscudé, G., Ribas, I., et al. 2021, A&A, 645, A58 [EDP Sciences] [Google Scholar]
 Petit, A. C., Laskar, J., & Boué, G. 2017, A&A, 607, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Queloz, D., Bouchy, F., Moutou, C., et al. 2009, A&A, 506, 303 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Quirrenbach, A., Amado, P. J., Mandel, H., et al. 2010, SPIE, 7735, 773513 [NASA ADS] [Google Scholar]
 Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian Processes for Machine Learning (Cambridge: MIT Press) [Google Scholar]
 Ribas, I., Bolmont, E., Selsis, F., et al. 2016, A&A, 596, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ribas, I., Tuomi, M., Reiners, A., et al. 2018, Nature, 563, 365 [Google Scholar]
 Ricker, G. R., Latham, D. W., Vanderspek, R. K., et al. 2010, AAS Meeting Abs., #215, 450.06 [Google Scholar]
 Santos, N. C., Mortier, A., Faria, J. P., et al. 2014, A&A, 566, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Scargle, J. D. 1982, ApJ, 263, 835 [Google Scholar]
 Seager, S., Dotson, R., & Lunar and Planetary, 2010, Exoplanets, The University of Arizona Space Science Series (Tucson: University of Arizona Press) [Google Scholar]
 Silva, A. M., Faria, J. P., Santos, N. C., et al. 2021, A&A, submitted [Google Scholar]
 Suárez Mascareño, A., Rebolo, R., & González Hernández, J. I. 2016, A&A, 595, A12 [Google Scholar]
 Suárez Mascareño, A., Rebolo, R., González Hernández, J. I., & Esposito, M. 2017, MNRAS, 468, 4772 [Google Scholar]
 Suárez Mascareño, A., Faria, J. P., Figueira, P., et al. 2020, A&A, 639, A77 [Google Scholar]
 Tuomi, M., Jones, H. R. A., Butler, R. P., et al. 2019, AAS, submitted [astroph] [arXiv:1906.04644] [Google Scholar]
 Turbet, M., Leconte, J., Selsis, F., et al. 2016, A&A, 596, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 UlmerMoll, S., Santos, N. C., Figueira, P., Brinchmann, J., & Faria, J. P. 2019, A&A, 630, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vida, K., Oláh, K., Kővári, Z., et al. 2019, ApJ, 884, 160 [Google Scholar]
 Wargelin, B. J., Saar, S. H., Pojmański, G., Drake, J. J., & Kashyap, V. L. 2017, MNRAS, 464, 3281 [NASA ADS] [CrossRef] [Google Scholar]
 Wildi, F., Pepe, F., Chazelas, B., Lo Curto, G., & Lovis, C. 2010, SPIE, 7735, 77354X [Google Scholar]
 Wildi, F., Blind, N., Reshetov, V., et al. 2017, SPIE, 10400, 1040018 [NASA ADS] [Google Scholar]
 Zechmeister, M., Reiners, A., Amado, P. J., et al. 2018, A&A, 609, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
Evidence (ln Z) and Bayes factors (Δ ln Z) for models with a given number of Keplerians, N_{p}, from the analysis of the CCF and TM RVs.
Posterior estimates for the parameters of the twoplanet model from the analysis of CCF and TM RVs. For each parameter we show the MAP estimate and the 68% quantiles of the distribution.
All Figures
Fig. 1 Comparison between RVs derived with the CCF and TM techniques. Top panels: two RV time series, and the middleleft and middleright panels show the generalised LombScargle (GLS) periodograms of CCF and TM RVs, respectively.In the bottom panels we display the histograms of the RVs (left) and of the RV uncertainties (right) for the two time series. 

In the text 
Fig. 2 Prewhitening procedure applied to the CCF RVs. Top panel: observed RVs together with the GP prediction. The periodogram of the residuals from this fit is shown in panel a, and the periodograms after two successive sinusoidal fits are in panels b and c. The false alarm probability of 1%, calculated with bootstrap randomisation, is shown by the horizontal grey lines. 

In the text 
Fig. 3 Posterior distribution for the orbital periods in the twoplanet model from the analysis of the CCF RVs (top) and TM RVs (bottom). The posteriors are shown normalised by the ESS per histogram bin (so the maximum possible value in the abscissa is 1). The prior is loguniform from 1 day to the time span of the data, which is marked with a dashed line. 

In the text 
Fig. 4 Joint and marginal posteriors for the orbital periods, semiamplitudes, eccentricities, and minimum planet masses, M_{p} sin i, in Earth masses, of the 11day (left) and 5day (right) Keplerian signals. These results come from the analysis of the TM RVs. The estimates at the top of each panel correspond to the median and 68% quantiles of the posteriors. 

In the text 
Fig. 5 Templatematching RVs from ESPRESSO phasefolded to the orbital periods of Proxima b (left) and Proxima d (right). TheGP and background components of the model were subtracted assuming the MAP solution. The full fit to the RVs andFWHM is shown in Fig. C.1. 

In the text 
Fig. 6 Prior and posterior distributions for the hyperparameters of the GP in the twoplanet model. Posterior estimates (median and 68% credible intervals) are show in each panel. It should be noted that for η_{3} and η_{4} the full support of the prior is larger than what is shown. 

In the text 
Fig. 7 Evolution of the Bayes factor between models with one and two Keplerians (top) and of the posterior distribution for the orbital periods (middle) and semiamplitudes (bottom) with an increasing number of ESPRESSO observations. The dashed vertical line corresponds to the number of observations available to SM2020. 

In the text 
Fig. 8 Posterior distributions for the periods of sinusoids in the analysis of activity indicators. The distributions are normalised and shown on the same vertical scale, which allows for a visual comparison of the significance of the periodicities. For example, the CCF contrast shows more significant periodicities (close to 100 and 40 days) than the Hα index. The stellar rotation period and its first and second harmonics are denoted with filled lines, and the 11day and 5day signals detected in the RVs with dashed red lines. 

In the text 
Fig. 9 Analysis of chromatic TM RVs. Three columns: results for the blue, green, and red regions, corresponding to the wavelengths between 440–570 nm, 570–690 nm, and 730–790 nm, respectively. Top and middle panels: ESPRESSO RVs phasefolded at the 11 and 5day periods, together with the maximum likelihood solution for each signal. Bottom panels: posteriors for the GP amplitude, η_{1} RV, with the posterior median and 68% quantiles as labels. 

In the text 
Fig. A.1 Posterior distributions for the hyperparameters of the GP kernels from the analysis of the TM RVs. We quote the posterior median and 68% quantiles for each parameter and kernel. 

In the text 
Fig. A.2 Predictive mean of the GP component of the model versus ro tation phase, assuming a value of η_{3} = 85.1 days. The top and bottom panels show the predictions for the RVs and the FWHM, respectively. The curves are only shown for times in which there are ESPRESSO ob servations, and their colour reflects the time since the first observation. 

In the text 
Fig. C.1 Maximum a posteriori solution for the twoplanet model on the TM RVs. The top panels show the RV observations for ESPRESSO18, ESPRESSO19, and ESPRESSO21, together with the GP component of the model (in pink) as well as the full model (in black). The RV residuals are shown just below, highlighting the full residual rms. The two bottom panels show the same for the FWHM. 

In the text 
Fig. C.2 Joint and marginal posteriors for the orbital periods, semiamplitudes, eccentricities, and minimum planet masses, M_{p} sin i, in Earth masses, of the 11day (left) and 5day (right) Keplerian signals. These results come from the analysis of the CCF RVs. The estimates at the top of each panel correspond to the median and 68% quantiles of the posteriors. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.