Measuring masses of transiting Earth-like planets: first insights from the PLATO RV data challenge

Introduction

The transit method is one of the most efficient and prolific techniques for detecting exoplanets today. Previous space missions such as Kepler (launched in 2009, Koch et al. 2010; Borucki et al. 2010) and TESS (launched in 2018, Ricker et al. 2015) have considerably improved our understanding of planetary demographics, particularly for short-period planets (Petigura et al. 2013; Batalha 2014; Burke et al. 2015; Bryson et al. 2020; Royer & Steffen 2025; Cui et al. 2026).

The transiting configuration naturally biases detections towards short-period planets. However, thanks to their large fields of view monitoring thousands of stars simultaneously, both missions have produced numerous and prominent detections. This success is largely due to the fact that the signature of a transiting planet is highly orthogonal to most stellar photometric variability signals, such as granulation and rotational modulation. However, these two missions have shortcomings. Most stars in the Kepler field are faint, which makes the mass measurements difficult. Whereas TESS targets are brighter, the observation strategy focuses on short period planets.

PLATO (PLAnetary Transits and Oscillations of stars) is an upcoming space-based photometric mission designed to observe bright stars, such that follow-up observations are possible, and its observation strategy is designed to enable the detection of Earth analogs (Rauer et al. 2014, 2024, 2025). Indeed, the survey strategy is broadly similar to that of Kepler, with a long-duration pointing towards a southern field observed for about at least two years. This baseline enables the detection of at least two transits of Earth-like planets, or potentially single-transit (mono-transit) events. The nominal mission is expected to continue for an additional two years, either by switching to a northern field or by extending the southern observations, resulting in a total time baseline comparable to Kepler’s four-year monitoring. Furthermore, its photometric precision is expected to be close of Kepler, with an anticipated noise level of about 20–30 ppm on hourly timescales for bright targets (Jannsen et al. 2024). Unlike Kepler, however, PLATO will primarily target brighter stars (Cabrera et al. 2026; Montalto et al. 2026), a key advantage for enabling precise ground-based radial-velocity follow-up and mass characterization of the detected planets.

The radial velocity (RV) follow-up is an integral part of the mission and aims at measuring the mass of the exoplanets detected by PLATO, including Earth analogs. Combined with the radius estimate, this will enable a mean density measurement, and yield constraints on the interior structure. Breaking the degeneracy between several interior compositions requires a high accuracy on the density, ideally $\sim10\%$ for terrestrial planets (Plotnykov & Valencia 2024), which translates to stringent requirements on the measurement of RV semi-amplitude, as seen below.

The PLATO mission requirements on mass measurements

The PLATO mission requirements R-SCI-Lo-15² reads: PLATO shall observe stars that are bright enough to allow for a determination through radial velocity measurements of the mass of a terrestrial planet orbiting a G0V star with an accuracy of 10% or better. “Terrestrial planet” is defined as having $R<2\,R_{\oplus}$ and $M<10\,M_{\oplus}$ – a goal that remains challenging for current RV capabilities.

In principle, the recent ultra-stable spectrographs are capable of reaching such a precision with $\mathrm{m\,s^{-1}}$ stability. These include ESPRESSO (Pepe et al. 2014), NEID (Schwab et al. 2018), EXPRES (Jurgenson et al. 2016), HARPS-N (Cosentino et al. 2012), MAROON-X (Seifahrt et al. 2018) and KPF (Gibson et al. 2016), building upon the legacy of instruments such as HARPS (Mayor et al. 2003). For instance, ESPRESSO has demonstrated an intra-night precision of the order of 10 cm/s (Pepe et al. 2021; Figueira et al. 2025). Unfortunately, at these precision levels, stellar activity dominates the RV signal, significantly complicating the detection of Earth-like planets around Sun-like stars (see e.g. Meunier & Lagrange 2020a; Crass et al. 2021; Hara & Ford 2023).

The challenge in RVs is greater than in photometry because the planetary signal is not a so-called gated signal (i.e., a component that appears and disappears over time), but instead is continuously superposed with all other contributions. This has important implications: for instance, machine-learning algorithms cannot be trained on subsets of the data that are “free” of the planetary signal, but also it means that any incomplete or inaccurate model will always interact with the planetary signals in unpredictable ways. In addition, the contaminating signals are more numerous, often larger in amplitude, and more problematic in the time and frequency domains (Meunier et al. 2015; Fischer et al. 2016; Cretignier et al. 2021; Crass et al. 2021; Hara & Ford 2023).

Some attempts have been made to evaluate the RV limitations based on simulations (Meunier & Lagrange 2019, 2020a,b), but the results are subject to caution. On the one hand, the simulations might not be realistic enough, on the other hand, the data analysis methods used in these studies might be too simplistic. From being discouraging, they motivate further intensive investigations based on real data sets and gathering more independent experts: ultimately the organisation of an RV challenge.

In this paper, we aim to assess whether the PLATO mass-precision requirements can be fulfilled with current methods and the extent to which they can be met across the $(K,P)$ parameter space. We have therefore designed an RV challenge involving several competitors and collaborators within the PLATO community, that used some state-of-the-art methodologies currently available. This challenge aims to mimic as closely as possible a PLATO follow-up. One of its most important aspects is that it makes use of real data, where planets are injected at the spectral level.

In Sect. 2, we introduce the main concept of the challenge. We begin by reviewing previous challenges (Sect. 2.1), followed by a description of the design of the present challenge, including the injected planetary systems and the datasets used (Sect. 2.2). The participating teams and their methods are then presented in Sect. 2.3. The results are described in Sect. 3, including the definition of the scoring metrics (Sect. 3.1), the performance on the different datasets (Sect. 3.2), and the overall ranking of the methods (Sect. 3.3). We then discuss the broader implications for the RV field (Sect. 3.4) before concluding in Sect. 4.

A PLATO RV challenge

Previous RV challenges

Few RV challenges have been conducted over the past decade. The main advantage of such community exercises is that they allow different methods to be tested on identical datasets. This enables a meaningful benchmarking of the field, provides a snapshot of the current state of the art, and allows the performance of the various approaches to be assessed within a consistent framework.

The first RV challenge was organized by Dumusque et al. (2017) using simulations as well as real HARPS data sets in which planetary signals were injected directly into the $RV(t)$ time series. Although realistic, the simulations were not capturing all the intricacies of stellar variability. Indeed, a linear fit of activity indicators worked extremely well to remove stellar activity signals in simulated systems, but performed worse in practice. Furthermore, this challenge was unfortunately limited to high-level data products from the data reduction software (DRS), namely the extracted RV time series and common activity proxies. In recent years, however, several studies have suggested that significant improvements can be achieved by working directly in the velocity domain [129] (Collier Cameron et al. 2021a; Zhao & Tinney 2020; Zhao et al. 2022b; Cretignier et al. 2022; Yu et al. 2024; Klein et al. 2025), in the wavelength domain (Cretignier et al. 2021; Shahaf & Zackay 2023; Zhao et al. 2024b) or by combining multiple RV time series derived from different spectral line selections (Lienhard et al. 2022; Dravins & Ludwig 2024; Al Moulla et al. 2024) or obtained with different extraction methods (Zechmeister et al. 2018; Silva et al. 2022).

More recently, Zhao et al. (2022a) organized another challenge based on four EXPRES data sets, comparing different methodologies. In this case, data products at multiple processing levels were made available to participants, offering a broader and more up-to-date perspective on the field. However, the conclusions were less definitive. First, the metric used to compare methods is the root mean square of the corrected RVs, which might favor methods that are overfitting (Hara & Ford 2023). Second, no participant consistently outperformed the others across all four data sets, which complicates any attempt at establishing a robust ranking metric.

Although it is well known that the true planetary content of real data sets can never be fully known since low undetected signals may always exists ³, the absence of injected signals in any of the data sets represents a limitation of the analysis. Injecting realistic planetary signals at the low-level data stage is technically challenging — and propagating those signals consistently through to high-level products is even more complex. However, we are aware that such an attempt is in progress (Zhao et al. in prep.).

Design of the PLATO challenge

We designed the PLATO RV challenge to closely mimic the conditions expected for the RV follow-up of PLATO targets, implementing some recommendations of Hara & Ford (2023) on the design of RV challenges. To achieve this, several high-level design choices were made.

First, the challenge exclusively relies on real data sets. While simulated data allow complete control over the ground truth, they inevitably lack certain complexities present in real observations, such as instrumental systematics and subtle correlated noise sources. These effects remain part of the real problem and should not be underestimated.

Second, planetary signals corresponding to transiting planets are injected at the low-level data stage on the spectra, rather than directly into derived RV time series. This ensures that injected signals propagate consistently through the full reduction and analysis chain.

Third, to place the RV challenge within PLATO’s core objective of detecting Earth-mass planets in the habitable zone (HZ) of Sun-like stars, we ensured that at least one transiting exoplanet is always located in the Sun’s HZ. This guarantees that the target star would be considered a high-priority candidate for follow-up within the Ground-based Observation Program (GOP) of PLATO.

Finally, additional planetary signals from non-transiting planets are injected. This enables the methods to be evaluated under realistic conditions, where both known (transiting) and unknown (non-transiting) signals coexist. As a reminder, in the Solar System no more than three planets can be seen transiting from any given vantage point in the Universe (Wells et al. 2018). This is important as the number of measurements needed to estimate the mass at a given precision of transiting exoplanets increases if non-transiting companions are present (He et al. 2021). Also, the detection of non-transiting planets is a more stringent test of RV methods. While it seems possible for a data analysis method to retrieve by chance a small semi-amplitude, close to the one of an injected planet, the odds of retrieving by chance the correct period of the injected signal are very low.

Injected planetary system

The parameter space of multi-Keplerian systems is extremely large. In this challenge, we therefore restricted ourselves to the production of three multi-planetary data sets and limited the analysis to purely circular orbits. Without loss of generality, the term Keplerian will be used hereafter, even if it should be understood as a sinusoid, rather than a strict Keplerian. All the participants were aware of the null eccentricity hypothesis and therefore removed it from their Keplerian fit adjustment. This is certainly one of the least realistic aspects of the challenge, since in this configuration only the semi-amplitudes of the transiting signals need to be adjusted. However, as shown in Shen & Turner (2008) and Baluev (2015), it is not significantly more difficult to detect eccentric planets than circular ones up to $e \approx0.7$. This is also the limit from which computational schemes might struggle to estimate eccentricities (Hara et al. 2019). There is no clear evidence that such eccentricities are common for low mass planets.

Transit measurements are known to yield far more precise constraints on phase $T_0$ and orbital period $P$ than RV observations. In practice, this precision would differ from system to system, adding an additional layer of complexity to the challenge analysis — particularly given the limited number of data sets we generated. Therefore, for all the transiting planets, we provided the exact orbital periods and transit phases such that all the participants set their priors on them to an arbitrary large precision. The performance of the teams will be solely evaluated based on their ability to retrieve the semi-amplitudes of the transiting planets, and find the transiting signals.

Based on Doppler shifts injection tests performed in Appendix B, we restricted the values of the $K$ semi-amplitudes below $\sim25$ m/s for the injected planets. Above this threshold, the Doppler wavelength shift interpolation used to inject the planets in the spectra may lose their accurate validity to a degree that is critical for the RV challenge.

The three injected multi-planetary systems are represented in the $(K, P)$ parameter space in Figure 1 and are hereafter referred to as systems A, B, and C. The Keplerian parameters are indicated in Table 1. All terrestrial planets range between 2 to 7 $M_{\oplus}$, that are already challenging planets to characterize today.

Syst.	$P$ [days]	$K$ [m/s]	$m_p$ [$M_{\oplus}$]	$T_c$ [BJD]	Type	$K_{\mathrm{raw}}$ [m/s] (%)
A	26.47	0.53	2.49	57152.2818	T	0.11 (21%)
A	53.90	0.50	2.97	57163.1609	T	0.15 (30%)
A	105.5	1.88	13.9	57166.2604	-	0.11 (6%)
A	373.0	0.30	3.4	57142.0194	T	0.12 (40%)
A	627.3	1.82	24.4	57186.9142	-	0.27 (15%)
A	921.5	7.48	114	57000.0000	-	0.92 (12%)
B	41.21	0.38	2.06	54971.5052	T	0.05 (13%)
B	52.10	0.71	4.17	54971.1748	T	0.11 (15%)
B	90.74	0.50	3.51	54981.4550	T	0.20 (40%)
B	128.0	0.50	3.97	54987.1590	-	0.20 (40%)
B	186.8	0.34	3.05	54964.6487	-	0.15 (44%)
B	473.5	0.56	6.78	55120.2901	M	0.29 (52%)
C	63.79	0.22	1.38	59640.1579	T	0.13 (59%)
C	369.5	0.19	2.18	59646.0919	T	0.35 (184%)
C	615.9	17.1	228	59641.5763	-	0.80 (5%)
C	3286	20.6	480	59656.0921	-	2.86 (14%)

For all the systems, we always assumed the solar properties for the host star, including the HZ limits even if we extrapolated the results to cooler type stars later (see Sect. 3.5).

System A

System A is a hierarchical six-planet system composed of three transiting terrestrial planets, two non-transiting Neptune-sized planets, and one non-transiting Saturn-mass planet. The system architecture is inspired by the HIP 41378 system, which hosts both transiting and non-transiting planets spanning a wide range of orbital periods, and thus provides a realistic benchmark for complex planetary configurations (e.g. Vanderburg et al. 2016; Santern et al. 2019; Grouffal et al. 2025). To adapt this architecture to our study, we rescaled the orbital periods by a factor of 1.7 for all planets, except for the fourth planet.

As a result, the orbital period of the innermost planet is close to the solar rotation period around 27 days, providing a stringent test for disentangling planetary signals from stellar activity due to increased signal crosstalk in both photometric and RV time series.

The third planet is non-transiting, similarly to HIP 41378 g, and contributes only through dynamical and RV signatures. The outermost planet is also non-transiting and induces a long-period dominant signal in the RV data.

The fourth planet was adjusted to transit and to lie within the HZ of a solar-type (G2V) star, corresponding to orbital periods between 338 and 586 days. To satisfy this constraint while preserving the system hierarchy, a specific scaling factor of 1.36 was applied to its orbital period. This planet is designed to mimic a high-priority PLATO candidate suitable for follow-up observations.

Planetary radii were scaled by factors between 2 and 2.5 in order to control the resulting RV semi-amplitudes. A relation $m_p \propto R_p^{3.7}$ was assumed to convert radii to masses following Zeng et al. 2016. By doing so, we ensured that all injected signals remain below 1 m/s, consistent with the low-amplitude regime expected for PLATO targets.

2.2.1.2 System B System B is a compact six-planet system composed exclusively of rocky planets, inspired by the Kepler-11 system (Lissauer et al. 2011). Such systems are representative of tightly packed inner planetary systems commonly observed by transit surveys, making them a relevant benchmark for testing detection performance in crowded configurations.

The system includes four transiting planets (the three innermost and the outermost) and two non-transiting planets. The outermost planet is located within the HZ. In this configuration, the habitable-zone planet produces a mono-transit event. Only the transit epoch was provided, while the orbital period was left unconstrained⁴, mimicking realistic observational limitations.

The planetary radii were scaled by factors of 1.6, 2.1, 2.4, 2.5, 1.9, and 2.1, while all orbital periods were increased by a factor of 4. This makes all the RV semi-amplitudes below 70 cm/s, placing the system in a regime where stellar activity dominates the RV signal, thereby significantly complicating planet detection and characterization.

2.2.1.3 System C System C is a four-planet system composed of two transiting terrestrial planets and two long-period non-transiting giant planets.

The system architecture was inspired by the HD 110067 system (Luque et al. 2023), with modified orbital periods and two planets removed to simplify the configuration (planets d, e, f, g only). Their orbital periods were rescaled by factors of 7, 12, 15, and 60. Planetary radii were kept unchanged, except for the innermost planet, whose radius was reduced by a factor of 2.

The two inner planets are terrestrial and transiting. The second planet is located within the HZ to mimic a high-priority PLATO target for follow-up observations. The two outer planets are long-period, non-transiting giant planets.

Despite the reduced number of injected planets compared to Systems A and B, the terrestrial signals in System C are among the most challenging to detect, with RV semi-amplitudes as low as 20 cm/s.

The outermost giant planet has an orbital period close to the 11-year solar magnetic cycle (Schwabe cycle), introducing additional difficulty in disentangling long-period planetary signals from stellar magnetic variability.

HARPS-N solar data set

To create the RV challenge, we used the HARPS-N solar telescope data set (Dumusque et al. 2015, 2021). HARPS-N has been observing the Sun daily since 2015 with 5-minute integration times continuously during the day. We used the data release corresponding to Dumusque et al. (2021), although a more recent release became available this year (Dumusque et al. 2026). The baseline of the observations span eight years between July 2015 (BJD = 2,457,220) up to October 2023 (BJD = 2,460,238).

The use of HARPS-N solar data offers several advantages compared to archival stellar RV databases. A key benefit is that, for the Sun, the planetary signals of the Solar System can be removed, thereby providing a well-defined ground truth. In addition, the exceptional temporal sampling of the solar observations enables flexible control of the cadence, allowing it to be reshaped in a realistic manner to closely match the expected sampling of the PLATO RV follow-up.

However, several important caveats must be considered when extrapolating solar results to stellar observations.

First, although difficult to quantify, it should be noted that some participants may have previously used or worked with the HARPS-N solar dataset. This potentially affects the strictly “blind” nature of the challenge⁵, and may therefore introduce a slight imbalance between participants. This limitation is arguably an inherent trade-off when relying on a dataset of such importance, which has been widely used in the literature over the past years.

The datasets used in the challenge have an average visible bandpass⁶ signal-to-noise ratio of $SNR = 450 \pm40$. The impact of lower S/N on the results was not explored, although such an analysis would be valuable and could be performed through the injection of photon noise.

Third, because the Earth orbits the Sun, Sun-as-a-star data sets are less sensitive to some instrumental systematics that are fixed in the laboratory rest frame. When observing the Sun, the barycentric Earth radial velocity (BERV) span is limited to approximately $\sim2$ km/s (due to the Earth’s orbital eccentricity), whereas for typical stellar observations the BERV span can reach dozens of km/s. Given that the typical photospheric line width of a slow rotating star is $\sim6$ km/s for an HARPS-like instrumental resolution, the instrumental systematics in solar observations remain effectively “frozen” with respect to the line profile, whereas in stellar observations they are swept across the spectrum over time. This difference allows instrumental effects couple to the spectral lines (see e.g. Silva et al. 2025).

Last but not least, the solar telescope follows a different optical path than standard stellar observations. On the one hand, this removes some sources of systematic related to telescope guiding (see e.g. the recently unperfect fibre scrambling identified on HARPS-N in Anna John et al. 2025). On the other hand, the disk-resolved nature of the Sun introduces additional complications, particularly differential extinction across the solar disk. Although this effect has been calibrated at the RV level through an analytical correction term $RV_{\mathrm{corr}}(t)$ in Collier Cameron et al. (2019), its correction at the spectral level has never been fully demonstrated or implemented. As a result, a simple Doppler shift by the measured $RV(t)$ is not strictly sufficient to correct the spectra.

Aware of this issue, we limited this effect doing the following selection of the spectra. First of all, we only selected spectra having a good quality flag. Bad data are usually identified as outliers in the daily magnitude vs. airmass relation which is expected to be linear (Al Moulla et al. 2023) and as such diagnose the periods in time when clouds occulted the Sun. We restricted the effect of the differential extinction by only selecting solar observation close to the meridian with an airmass smaller than $Z = 1.3$. At such airmass, the differential extinction effect is lower than 10 cm/s according to the $RV_{\mathrm{corr}}(t)$ vector. Note that this choice also recreates a realistic season length sampling since the Sun in winter is not going below an airmass $Z = 1.3$.

Each day, we randomly selected three consecutive observations that satisfy the previous requirements to create an effective 15 minutes exposure time which mimics the expected exposure time for PLATO targets and strongly mitigates the contribution from the 5 minutes p-modes pulsations (Chaplin et al. 2019). Stacking the three consecutive spectra produces an effective visual bandpass SNR $\sim450$ or 15 cm/s of RV photonic precision. As an important note, the Sunas-a-star on HARPS-N$^{6}$ is equivalent to a magnitude $m_{\mathrm{V}} = 5.4$ on the 3.6m TNG telescope or to a magnitude $m_{\mathrm{V}} = 7.2$ on the 8.2m VLT telescope. Such magnitude represents the bright tail of the PLATO sample (Nascimbeni et al. 2026).

In order to avoid obvious reverse engineering (such as doing the difference of two RV time series), three independent realisations of spectra selection was done for each of the planetary systems described in Sect. 2.2.1. Each data set contains a different number of observations. We will call those three data sets hereafter A1050, B1300 and C800 related to the respective planetary system and number of data points each of those data set have (see Figure 3 for the exact number).

Note that the data set A, B and C represents a huge time investment with respectively 26, 33, and 21 equivalent nights of telescope time which is realistically the amount of time needed to be able measuring such HZ signals. The season length $S$ of each data set is also slightly different even though longer than half-a-year, with $S = 225$, $S = 236$ and $S = 207$ days. This is roughly two months shorter than the season length expected from PLATO field from Paranal ($S \sim280$ days).

We displayed in Figure 2 the Generalized Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009) in semi-amplitude of the three RV data sets without the planets injected and also indicate the ratio between the planetary semi-amplitude and the native semi-amplitude power at this period. Similarly, we fitted the planetary systems on the RVs without the planet injected in order to quantify the initial amount of power and degeneracy between the injected systems and the data sets. We reported the semi-amplitudes $K_{\mathrm{raw}}$ in Table 1, larger is the number, higher is the cross-talk between the native power spectrum and planetary signals.

The RV time-series of the three data sets with the planets injected are shown in Figure 3. We also displayed as a comparison below the MgII index⁷ converted in plage filling factor using the time-series published in Cretignier et al. (2024a). We also displayed the solar TSI (see Sect. 2.2.2.5) that is a mixture of faculae and sunspot filling factor.

Processing of the data

Injection of the planetary signals at the spectrum level may be challenging since to be exact, the planetary signal should be injected on spectra free of instrumental systematics. Indeed, an incorrect processing that would strongly underestimate the EPRV problem would be to simply shift the spectra by the Keplerian solution, injecting the planetary signals on the laboratory rest-frame systematics.

We used YARARA (Cretignier et al. (2021)) that is a post-processing pipeline analysing spectra time series to correct instrumental systematics at the spectrum level. As a short recall, YARARA starts by a preliminary standardization of the S1D order-merged spectra by fitting their continuum with the RASSINE code § (Cretignier et al. 2020b). During this stage, all the spectra are also interpolated linearly on a common wavelength grid in order to analyse and correct each wavelength bins time series independently. The code then creates a master spectrum, and fit on the residual spectra time series multi-linear bases from proxies extracted by the pipeline itself.

The YARARA correction of HARPS-N includes the corrections of cosmics rays, telluric lines, ghosts (Dumusque et al. 2015), interference pattern (e.g Thompson et al. 2020), instrumental point spread function (PSF) distortions (see Stalport et al. (2023) or Klein et al. (2024)) such as those introduced by the change of the solar $v\sin i$ along the year (Collier Cameron et al. 2019) and stellar activity. Note that any correction can be canceled and added back on the spectra easily. This is mandatory since the stellar activity and change of the PSF should be on the spectra prior to their RV shift by the Keplerian solution.

For each data set, we first ran YARARA on the spectra time series, remove the stellar activity and PSF correction, shift the spectra according to the Keplerian solutions, before adding back all the contaminations. We then ran a second time the pipeline on this new data set and offer to the participant two flavours of the data: a) the $YV0$ flavour that is uncorrected of anything and b) the $YVA$ flavour corrected of all the instrumental systematics except the stellar activity. Participants were free to choose between the two versions of the data they wanted to use. The data ($\sim15$ Gb) are freely accessible on the Google Drive⁹ and Zenodo¹⁰.

HARPS-N solar S-index

We provide the S-index time series extracted from the HARPS-N spectra, which is a standard diagnostic typically delivered alongside RVs for visible-band spectrographs. The S-index is directly related to the filling factor of active regions (Cretignier et al. 2024a), although the line response remains relatively weak for solar-like activity, requiring careful spectral extraction.

Indeed, the Sun-as-a-star observations obtained with HARPS-N are known to be affected by significant contamination from optical ghosts crossing the cores of the Ca II H&K lines, from which this activity proxy is derived (Dumusque et al. 2021).

For the Sun, this contamination would in principle be of limited concern, as instrumental systematics are expected to remain largely stable over time. However, regular warm-up cycles of the HARPS-N detector modify the relative intensities of these ghosts reflection with the main physical orders, introducing not only an absolute contamination but also a time-variable component.

This contamination can be mitigated by estimating the flux between physical spectral orders on the CCD, as described in Dumusque et al. (2021), or by using the YARARA framework, which includes dedicated corrections for such effects (Cretignier et al. 2021). We opted for this latter for simplicity.

In addition, the contamination is further reduced by applying the emission-profile extraction method presented in Cretignier et al. (2024b). In contrast to the classical Mount Wilson S-index $S_{MW}$, which relies on a triangular bandpass integration, this approach directly fits on the spectra the emission profile of plage regions to the data (see Oranje 1983; Cretignier et al. 2024a). This method provides improved uncertainties and is less sensitive to outliers compared to standard bandpass integration. In this framework, the plage emission profile used for the extraction can be interpreted as an optimally weighted integration of the signal.

Note that, due to the spectrograph spectral efficiency response of HARPS-N, the SNR at 550 nm is equivalent to twice the value of the one in the continuum region surrounding the CaII H&K lines at 420 nm.

Solar TSI

We provided the total solar irradiance (TSI) as obtained from the LISIRD¹¹ database. We did not inject the different transits in the TSI because this later is not necessarily representative of the precision expected from the PLATO lightcurves, adding an extra layer of complexity in the challenge and its possible interpretation that we wanted to avoid. Instead of that, as explained previously, we provided the exact phase and period of the transiting systems and mainly provided the TSI as a proxy that participants could use to fit their stellar activity model.

ARVE solar photospheric components

In order to benchmark the RV challenge and assess performance on a less demanding dataset, we generated a simplified red-noise model free of instrumental systematics and active regions: it only contains white noise, oscillations, granulation and supergranulation. For this purpose, we used the ARVE code (Al Moulla 2025), which simulates realistic time series from their power spectrum, as fitted on the HARPS-N solar dataset. The individual RV contributions are shown in Figure 4 for a three weeks baseline. The model predicts a minimum total RV jitter of 86 cm/s, composed of 15 cm/s from photon noise, 5 cm/s from p-modes, 24 cm/s from granulation, and 62 cm/s from supergranulation. These contributions combine to produce the total RV budget of 86 cm/s similar to the budget found in Al Moulla et al. (2023) and Meunier & Lagrange (2020b). We note that we did not include potential long-term magnetic cycle modulation of these components. Recent work suggests that supergranulation properties may vary with magnetic activity (O’Sullivan et al. 2025), but this effect is not accounted for in the present model. Using ARVE, we simulated the same 8-year RV time series with 15-minute integration times as the three data sets. This theoretical red-noise time series will be used (see next Sect. 2.3.1) for benchmarking purposes and represents an optimistic case scenario.

Participants

In total, eight participants participated in the challenge, seven provided a full set of results and one participant provided a partial result analysing only one (B1300) of the three systems. We compared the results of the candidates with a fictive processing represented by a benchmark participant. The output expected from the participants were the $K_i$ samplers from their guess Keplerian solution. We present here below the methodologies of all the participants to deliver those products.

Benchmark (exact Keplerian solution + ARVE)

This fictive participant is not using the YV0 or YVA data sets, but the ARVE RV budget described in Sect. 2.2.2 that includes the same planetary signals on it. A Keplerian signal was added for each transiting planet with the period and phase specified by the challenge. Non-transiting planets were added if visible as a significant peak in the residuals periodograms (obtained after fitting the transiting planets amplitudes), which was the case for all the signals due to the lack of power for the ARVE time-series at long ($P > 3$ days) periods. All the eccentricities were assumed null and uniform prior was used on all the remaining parameters. Note that we also investigated a benchmark with free eccentricities that provided results very similar and we therefore stuck to the null eccentricity model to a fair comparison with the participants.

Under the assumption that no participant is able to deal with supergranulation (which is a question none of them have tried to answer so far), this participant represents the best performance that can be achieved in theory since the Keplerian model is close to the exact one and there is no long-term stellar activity and instrumental systematics. It will allow us to check if the amount of information in the time-series was sufficient or if solar photospheric components are sufficient to hinder the success of the challenge. All the MCMC analysis were done from the DACE¹² implemented tool. Transit phases and periods were kept fixed for the transiting planets during the MCMC exploration.

Linear regression of the S-index - (S. Aigrain)

The participant analysed both the YV0 and YVA data sets, assuming that all planets were on circular orbits, as the participants were informed of this a priori. The analysis consisted of two steps: 1) detecting the non-transiting planets, and 2) jointly fitting for the planetary and stellar signals.

In both steps, the stellar RV variations were assumed to be proportional to the variations of the S-index. In step 1, the RVs were first modelled as a linear combination of the S-index time-series, and sinusoidal terms corresponding to the known ephemerides of the transiting planets (ignoring mono-transits). The coefficients of the resulting, linear basis model were fitted using scikit-learn.LinearRegression¹³ (version 1.2.2), and the best fit model subtracted from the RVs. The GLS periodogram of the residuals was then computed using astropy.timeseries.LombScargle¹⁴ (version 6.1.7), on a regularly spaced grid of frequencies from $\nu= 0.0002$ to 0.9 cycles/day ($P = 1.1$–5000 days), with nyquist_factor = 10. If the highest peak in the periodogram had a False Alarm Probability FAP $\le10^{-5}$, the corresponding period was taken to correspond to an additional planet. A sine and cosine term at this new period were added to the linear basis model, and the process was repeated until no further signals were found below the FAP threshold.

The fitted coefficients of these additional terms were then converted to RV semi-amplitudes and $T_0$ values for the additional planets. For system 3, whose RV variations are dominated by non-transiting planets rather than activity signals or transiting planets, the initial fit was omitted and the initial GLS periodogram was computed directly on the YV0/YVA RVs.

Step 2 consisted of a Bayesian modelling of the planetary and stellar signals using Hamiltonian Monte Carlo (HMC) to explore the joint posterior over all the parameters. The free parameters of the fit were thus: the systemic RV$_0$ offset, the semi-amplitudes $K_i$ of all the planets, the periods $P_i$ and times of inferior conjunction $T_{0,i}$ of the additional planets identified at step 1 only, the coefficient $S$ of the S-index term, and a jitter term $\sigma_{\mathrm{jit}}$ added in quadrature to the RV uncertainties. This participant used jaxoplanet¹⁵ (version 0.1.0) to evaluate the RV signatures of the planets.

The HMC sampling was implemented in numpyro¹⁶ (version 0.18.0) using the No U-Turn Sampler (NUTS). For the semi-amplitudes, half-normal priors with a mean of zero and a standard deviation of 50 m/s were used. For the periods, the prior was a truncated log-normal with a standard deviation of 0.05 (in log), and bounds mid-way between the logarithms of the initial period guesses. This choice of prior allowed the model to explore a fairly wide range of periods for each planet without allowing the planets to "switch" order. A uniform prior between 0 and $2\pi$ was used for the orbital phase of the additional planets, using numpyro-ext.distx.Angle16 (version 0.0.4) to sample this circular variable correctly. Wide normal priors (with a standard deviation of 10 m/s) were used over $RV_0$ and $S$, and a log-uniform prior between 0.1 and 3 m/s over $\sigma_{\mathrm{jit}}$. After a 1000-step burn-in phase, 2 chains of 10000 steps each were run, which was sufficient to ensure that the Gelman-Rubin statistic $R \leq1.01$ and the Effective Sample Size (ESS) was more than half of the number of samples for all parameters and all 3 systems. The resulting samples were provided as output along with various diagnostic plots (best-fit model and residuals, chains and corner plots).

Linear regression of the S-index - (M. Damasso)

The team analysed the YVA dataset only, assuming that all planets have circular orbits, which is a prior information provided to the participants. For all the targets, each test model included a linear regression with the $S_{\mathrm{MW}}$ index and the planetary Doppler signals, after verifying that a significant correlation exists between the RV and $S_{\mathrm{MW}}$ datasets. They also included a constant uncorrelated jitter term, which was added in quadrature to the nominal RV uncertainties.

A frequency content analysis of the RV data was initially performed with GLS to identify potential signals due to non-transiting planets in the original datasets and pre-whitened residuals, having a FAP<0.01% determined through a bootstrap analysis. These signals defined the maximum number of planets to be fitted for each simulated system through a Monte Carlo (MC) simulation. The MC modelling was performed by exploring the parameter space with the publicly available nested sampler and the Bayesian inference tool MultiNest v3.10 (e.g., Feroz et al. 2019), through the pyMultiNest wrapper (Buchner et al. 2014). included 500 live points, a sampling efficiency of 0.3, and a Bayesian evidence tolerance of 0.5. The package radvel (v1.4.9 Fulton et al. 2018) was used to model the RV Doppler planetary signals.

Bayesian model comparison was performed to by calculating the difference $\Delta\ln Z$ between the values of the Bayesian evidence $\ln Z$ given by MultiNest for each model, assuming that the a-priori probability is the same for each model, and following the scale in Feroz et al. (2011) to interpret the statistical significance $\Delta\ln Z$. Values of $\Delta\ln Z$ were calculated between pairs of models with a number of planets $n$ and $n+1$ The team adopted uninformative priors for all the free parameters, except for the orbital period and epoch of inferior conjunction of the transiting planets, for which they adopted Gaussian prior distributions using the information provided by the organisers. In particular, we explored the semi-amplitudes of the planetary Doppler signals within intervals with lower limit zero and upper limits from 5 to 50 ms$^{-1}$, depending on the $\chi$-square best-fit values calculated with GLS for each potential planetary signal.

The team analysed the residual signal in RVs and found a significant stellar activity component that could be fit with a quasi-periodic GP. However, inclusion of the GP itself with the previously fit Keplerian systems did not significantly affected their results and decision was made ultimately to adopt in this work models with no GP term included in order to provide results from a framework not too sophisticated that can be regarded as a benchmark against which to compare those obtained using more complex models.

Multi GPs - (L. Malavolta, D. Nadiello, G. Mantovan)

The participant analysed the YVA data set only and performed a joint modelling of photometric and spectroscopic data within a Bayesian framework using PyORBIT17 (Malavolta 2016; Malavolta et al. 2016, 2018). The RV signal of each planet was modelled using a Keplerian function. For known transiting planets, the eccentricity was fixed to zero and the argument of pericenter to 90°. Arbitrary tight Gaussian priors were imposed on the orbital period and the time of inferior conjunction around the provided values.

For planets with unknown periods, a range from 10 to $10^2$ days was explored in uniform log$_2$ space, with eccentricity allowed to vary up to 0.5. When the time of inferior conjunction was not available, the mean longitude was fitted instead, assuming a longitude of the ascending node equal to zero. The analysis included up to four additional non-transiting planets alongside the known transiting planet. For planets with known periods, the period and RV semi-amplitude were explored in linear space; otherwise, they were explored in log$_2$ space.

Stellar activity was mitigated using GP regression through the S+LEAF 2 package (Delisle et al. 2020, 2022). An Exponential-Sine Periodic (ESP) kernel (including up to four harmonics) was adopted as a computationally efficient approximation of the squared-exponential periodic kernel. For the spectroscopic dataset, a multivariate GP regression was implemented using full width at half maximum (FWHM), bisector-span (BIS), equivalent-width (EW), S-index, and H$\alpha$ as activity indicators to constrain the underlying process.

For the photometric data, three scenarios were considered: exclusion of the photometry; GP regression using the same kernel and shared hyperparameters as the spectroscopic dataset but with an independent covariance matrix; and full inclusion within the multivariate GP framework, treated analogously to the spectroscopic activity indicators.

The first derivative of the GP was included only for the RV, BIS, and EW time series (see e.g. Mantovan et al. 2024). To avoid degeneracies in the GP coefficients, the convective coefficient associated with the RV was constrained to be positive, while uninformative priors were adopted for all other parameters. Each dataset included independent offset and jitter terms. Additional tests considered both unidimensional GP models and configurations using only BIS and S-index as activity indicators.

The parameter space was explored using the affine-invariant MCMC ensemble sampler emcee (Foreman-Mackey et al. 2013), initialized from the results of the differential evolution optimizer pyDE. The choice of sampler was motivated by the need for controlled execution time. The optimisation phase was run for 50,000 generations with pyDE, followed by 100,000 MCMC steps with emcee, adopting a thinning factor of 100 and a burn-in of 80,000 steps. The number of chains was set dynamically to four times the dimensionality of the model.

Model selection was performed by comparing Bayesian Information Criterion (BIC) values derived from the RV-only likelihood, with the penalty term computed using only the parameters involved in the RV modelling. This approach was specifically implemented for the RV data challenge.

CCFs SCALPELS {#sec:2.3.5} - (A. Collier-Cameron)

This participant applied the SCALPELS algorithm on the CCFs. Since some instrumental systematics can also induce variations in the CCF shape, both the YV0 and YVA realisations of the CCF time series are analysed.

The autocorrelation function (ACF) of a spectrum (or of a CCF) is invariant under translation. The SCALPELS algorithm (Collier Cameron et al. 2021a) exploits this property to identify orthonormal modes of variation in the CCF shape that are insensitive to planetary Doppler shifts. The ACF, denoted $A$, is computed for each CCF in the time series using a circular shift and cross-multiplication, with each row normalised by its mean. A singular value decomposition (SVD) of the resulting time series is then performed:

The rows of the right-singular matrix $V^{T}$ represent orthonormal modes of variation in the ACF shape, while the columns of the left-singular matrix $U$ provide the corresponding orthonormal time-domain coefficients. The diagonal matrix $S$ contains eigenvalues ordered by decreasing contribution to the variance of $A$. In practice, these eigenvalues decrease rapidly, such that most of the signal can be captured by a low-rank reconstruction using only the first few principal components, typically four.

The solar RV time series contains a small number of outliers, generally arising from anomalous distortions of the CCF that are poorly described by a low-rank model. These are identified and removed using an iterative scheme based on low-rank reconstruction and outlier rejection.

The shape-driven component of the RV signal is obtained by projecting the RV time series onto a basis defined by the leading columns of the matrix $U$. These basis vectors are reordered according to their contribution to the Bayesian Information Criterion (BIC) of the RV fit. After this reordering, the first four columns of $U$ provide a satisfactory reduced-rank representation of the shape-driven RV signal without significant overfitting.

The search for planetary signals is carried out using the KIMA package (Faria et al. 2018). This framework simultaneously fits Keplerian orbital models together with linear combinations of the selected columns of $U$, treated as activity indicators, using the diffusive nested sampling algorithm DNEST4 (Brewer & Donovan 2015). Gaussian priors are adopted for the ephemerides of known transiting planets. When only the time of conjunction from a single transit is available, an uninformative prior is used for the orbital period, allowing the algorithm to identify any RV signal consistent with the observed transit time.

The trans-dimensional nature of the sampler allows the number of non-transiting planets to vary up to a maximum value $N_{\mathrm{p,max}}$. Bayesian odds ratios between models with different numbers of planets are derived from their relative occurrence in the posterior samples.

CCFs LINSTAR {#sec:2.3.6} - (J. Blanco-Pozo, M. Perger, I. Ribas)

The participants applied the LINSTAR (LINear Regression for STellar Activity Removal) method, which is an adaptation of the CANSTAR method developed by Blanco-Pozo et al. (2026) for this specific challenge. The CANSTAR activity correction relies on training a convolutional attention network on StarSim (Herrero et al. 2016) generated data; however, the lack of sufficiently realistic synthetic models of the Sun at the time of the RV challenge led to an adapta-

tion of the method into a purely data-driven approach. The method operated through a two-step procedure:

(i) Separation of line shifts from line-shape distortions in the CCF is achieved through the following decomposition:

where $I$ is the unity function, $T$ is the time-averaged template CCF, $T' = \mathrm{d}T/\mathrm{d}v$ is its first derivative (capturing the line shift $l(t)$), and $E$ is the residual error. The functions $G_n$ form a basis describing line-shape distortions. A Gram–Schmidt orthogonalisation is applied to ensure that $I(v)$, $T(v)$, and $G_n$ are strictly orthogonal to $T'$, enabling a clear separation between Doppler shifts ($l$) and pure line-shape distortions ([$c_0, c_1, g_2, \ldots, g_N$]).

(ii) A joint fit of Keplerian orbits and a linear combination of the selected line-shape distortion coefficients ([$c_0, c_1, g_2, \ldots, g_N$]) using the dynamic nested sampler dynasty (Speagle 2020).

For FGK stars, where $T$ is approximately Gaussian—as in the present study—a Hermite basis is adopted. This basis is constructed from the $n$-th derivatives of a Gaussian function fitted to $T$, starting from $n = 2$, since the $n = 0$ and $n = 1$ components correspond to $T$ and $T'$. As shown in Blanco-Pozo et al. (2026), these coefficients naturally map onto classical indicators ($c_1$ to contrast, $l$ to RV, $g_2$ to FWHM, $g_3$ to BIS), while extending to higher-order shape variability ($n \ge4$), which captures both symmetric and asymmetric distortions.

Following a SCALPELS-like strategy (Collier Cameron et al. 2021b), the optimal number of distortion coefficients is determined iteratively by minimizing the RV root-mean-square (RMS) while controlling model complexity through the Bayesian Information Criterion (BIC), thereby avoiding overfitting. Giant planets dominating the time-series variability are pre-whitened via a preliminary least-squares fit prior to this selection.

Once the optimal number of distortion coefficients is established, Keplerian signals are iteratively fitted to the data. For transiting planets, informative priors are applied on periods and phases. For non-transiting planets, uninformative priors are adopted for the orbital phase and semi-amplitude, while period priors are restricted to broad intervals around the dominant peaks in the periodogram to prevent overlap between signals in multi-planet fits.

The dynamic nested sampling is initialized with 500 live points, using a random-walk proposal to explore the parameter space. The sampler computes the log-Bayesian evidence ($\Delta\ln Z$) by integrating the likelihood over the prior volume, progressively focusing on regions of higher likelihood. The log-likelihood values of the sampled points are stored at each iteration to compute importance weights, which are then used to derive the posterior distributions of the model parameters.

Model complexity is increased iteratively, evaluating both $\Delta\ln Z$ and the residuals after subtracting the best-fitting model. The addition of further components is halted when no significant signals remain in the residual GLS periodogram (FAP $> 1\%$), or when the increase in model evidence no longer justifies the added complexity ($\Delta\ln Z \le5$; Trotta 2008).

SHELL CNN {#sec:2.3.7} - (Y. Zhao)

The shell-CNN is a convolutional neural network (CNN)-based method designed to model stellar activity signals at the level, thereby improving the detection of Earth-like planets (Zhao et al. 2024a). The approach relies on the distinction that planetary signals induce pure Doppler shifts in stellar spectra, while stellar activity produces variations in spectral line profiles, such as changes in asymmetry and depth.

For each dataset, the 1D solar spectra are projected into the spectral shell domain ($f_0$, $\frac{d f_0}{d v}$). A solar line mask derived from Cretignier et al. (2020a) is used to select strong, unblended spectral lines that are free from telluric contamination.

Next, for each dataset, we trained a convolutional neural network from scratch to model the stellar activity induced RVs, full width at half maximum (FWHM), and bisector span (BIS) using a 10-fold cross-validation scheme. Specifically, the spectral shell data were randomly divided into 10 groups. The model was then trained over 10 iterations, where in each iteration one group was designated as the test set while the remaining nine groups were used for training. The final fitting results were obtained by concatenating the predictions on the test sets from all iterations.

The trained neural network is used to estimate and subtract the stellar activity–induced RV component from the total RVs, which are measured using a classical CCF applied to the solar spectra. The planetary signal is assumed to follow a circular orbit. A periodogram of the resulting RV residuals is then computed, and a numerically estimated FAP is derived, adopting a detection threshold of 0.1%. A signal is considered detected if the corresponding peak exceeds this threshold. In such cases, a sinusoidal model is fitted to recover the amplitude and phase of the injected planet. Posterior distributions of the best-fit Keplerian parameters are subsequently obtained using the UltraNest package (Buchner 2021).

For datasets A1050 and C800, the RV time series extracted from the solar spectra exhibit strong planetary signals. These dominant signals are first modelled with a single-planet fit and removed prior to training the neural network on the residuals, in order to isolate the stellar activity component.

YARARA V2 - (M. Cretignier)

YARARA V2, presented in Cretignier et al. (2023), is an extension of YARARA V1 (Cretignier et al. 2021). The first version primarily consists of the corrections described in Sect. 2.2.3, applied directly in the wavelength domain at the spectral level. Once the spectra are corrected, line-by-line (LBL) RVs are extracted following Dumusque (2018). Time-domain proxies are then derived by projections onto compact spaces.

The first compact space, called the shell, corresponds to a change in the variables of the spectrum $f(\lambda)$ from the wavelength $\lambda$ to the flux gradient, $f\left(\frac{d f}{d \lambda}\right)$. This transformation, introduced in Cretignier et al. (2022), provides two main advantages. First, it preserves information on the normalised flux, which is an important diagnostic of stellar activity as an indirect proxy for photospheric depth (Cretignier et al. 2020a; Al Moulla et al. 2022, 2024). Second, in this representation, a pure Doppler shift has a well-defined form, referred to as the Doppler shell, $DS(f, d f/d\lambda)$. This transformation is derived from the first-order approximation of the Doppler shift. A principal component analysis (PCA) is performed in this space to obtain components $PC(f, d f/d\lambda)$. The number of components is selected using a modified $p$-leave-out cross-validation scheme (see Appendix of Cretignier et al. (2022)). A Gram–Schmidt orthogonalisation is then applied between these components and the Doppler shell before fitting them to the data. The resulting coefficients form time series that have

been shown to correlate with RV variations induced by activity or instrumental systematics.

The second compact space is constructed from the LBL RVs. Under the assumption that a planetary signal corresponds to a pure Doppler shift affecting all spectral lines identically, the variance across LBL RVs is free of planetary contributions (Cretignier et al. 2023; Ould-Elhkim et al. 2023). Consequently, a PCA performed on LBL RVs, which operates on mean-centered data and captures differential variations between spectral lines, should not contain planetary signals. This again defines a time-domain basis that can be fitted to the data.

YARARA V2 consists of fitting these two time-domain bases to the LBL RVs. Transiting planetary signals are explicitly included in the model to prevent absorption of their amplitudes and to enable the search for additional signals in the residuals. A sinusoidal signal search is performed using an iterative GLS periodogram with a 1% FAP threshold, limited to a maximum of six signals and orbital period range $P$ between 3 to 5000 days. Once the Keplerian model and the amplitudes of the basis vectors are fixed, the residual RVs are analysed using an MCMC framework, assuming circular orbits and fitting only for the amplitudes of the transiting signals with the DACE platform as for the benchmark in Sect. 2.3.1.

1D GP - (E.C. Gonçalves)

The participant only analysed the data set B1300 by lack of time. They used the KIMA code (Faria et al. 2018) to model the provided RVs with a quasi-periodic kernel. The model included 4 Keplerians for the transiting planets and a quadratic polynomial for the mean trend. The Keplerians had variable period $P$, $K$ semi-amplitude and $T_c$ time, the prior for the periods was gaussian with the mean being the given period value with a std of 0.1 (the mono-transiting case was set with a uniform prior from 100 days to 5 times the time span of the time series), the $K$ semi-amplitude had a uniform prior from 0 to 5 times the variation of the RVs and the prior for the $T_c$ were gaussian with the mean being the given value and an std of 0.001. (eccentricities were set to 0 and argument of periastron to pi/2).

A standard quasi-periodic kernel (Rasmussen & Williams 2006) was used for the GP with the following priors for the hyperparameters $\eta_1$ set to default by KIMA (uniform from 0 to the variation of the data max-min), $\eta_2$ set with a uniform prior between 15 and 600 days, $\eta_3$ set with a uniform prior between 15 and 50, $\eta_4$ was left as KIMA’s default which is a uniform prior between 0.2 and 5.

The model was run using diffusive nested sampling (Brewer & Foreman-Mackey 2016) for 1 000 000 steps with 16 threads.

Results of the challenge

In this section, we present the analysis of the three data sets processed by the participants. We start by defining a score metric (Sect. 3.1) that can be used in the context of the PLATO mission in order to compare participants (Sect. 3.2 and Sect. 3.3). Later, the results of all the participants will be put together in order to get broader conclusions for the RV field (Sect. 3.4).

Defining a score metric

The output asked to the participants were the posterior samples of the Keplerian fit on the three data sets. We will define a false detection as any fitted signal with a period farther than 10% of the injected orbital periods $P_{\mathrm{inj}}$. On the contrary, a missed signal will be an injected planetary signal for which no signal were closer than those 10% in period.

We then developed a metric using the semi-amplitude samples $K_i$ to evaluate the methods aiming to closely match the PLATO requirements written in Sect. 1.1. We note that there is an ambiguity in the requirement, it is unclear whether the 10% figure pertains to the precision or accuracy of the measurement. On the one hand, the accuracy cannot be checked in practice since the ground truth is never known. On the other hand, precise measurements can be biased, such that the confidence in the estimated value is overestimated. The present challenge provides a good benchmark to link the accuracy, that can be computed here since the ground truth is known, to the precision, computable also for real PLATO follow-up data.

In this work, we will define the accuracy by $\Delta K/K$ as the relative difference between the injected ground truth $K_{\mathrm{inj}}$ and the median of the posterior on the samples $K_i$:

Note that our accuracy definition is signed, but we also use an unsigned one. It will be clear when signed versus unsigned accuracy was used.

On the other hand, we defined the precision $\sigma_K$ as the sigma width of the posterior, measured by the median absolute deviation (MAD) of the samples, normalised by the ground truth injected value:

where the factor 1.48 is just a conversion factor from MAD to an usual standard deviation for a Gaussian distribution.

While it is trivial that neither accuracy or precision alone are sufficient to assess the score of a participant, it is certain that a method for which both the accuracy and precision are good should be considered as successful in general. We therefore defined the maximum error $\varepsilon_K$, of the challenge as the worst value between the two:

Note that with those definitions, the usual "number of sigma" $N_\sigma$ is just given by the ratio of the accuracy over the precision:

If in practice this maximum error $\varepsilon_K$ definition is an attempt to obtain a one-dimensional metric to rank the participants, as described later, this ruler is not perfect and we therefore need to always present the pair $(\Delta K/K,\sigma_K)$ or their alternative representation $(\varepsilon_K,N_\sigma)$ to get a global view of a method.

Results of the participants

In Figure 5, we display the $K_i$ posteriors distributions for the data set A1050 for all the correctly identified signals for each participant. The question of false positives will be discussed later. For later comparison (see Sect. 3.4), we also created the average posterior from all the participants. Each participant is represented by a unique color and marker shape displayed on the left and preserved along the paper.

Until the end of the paper, we defined two arbitrary levels as relevant for the PLATO mission. The first level is 10% level that matches the PLATO requirements, the second level is 25% maximum error that is the score level needs to achieve accurate planetary atmosphere modelling ((Seager & Deming 2010; Crass et al. 2021)). We display in the Figure 5, the results satisfying 10% by a green box, and the ones satisfying 25% by a yellow box.

As visible in the figure, for all the signals larger than 50 cm/s in amplitudes the 10% maximum error $\varepsilon_K$ is achieved by most participants. However, such a score is not achieved for the HZ transit planets that reach at best the 25% maximum error. On a positive note, the transiting planet close to the rotational period of the Sun is recovered at 10% by most participants.

In order to better represent the accuracy and precision, we extracted the value for each participant and represented all of them in Figure 6. The same was done for the data set B1300 and C800, respectively in Figure 7 and Figure 8. The 10% and 25% level are now represented by colored boxes. The score reached by the benchmark is also shown as a black box for comparison.

Such a visualization allow a rapid diagnostic of the easy and challenging planets. For all the data sets, a maximum error of 10% was never achieved for any of the HZ planets and only the HZ planets of the system A1050 reach the 25% maximum error. This confirms the results obtained by Meunier et al. (2019a) based on simulations that achieving such a precision level on the mass, even in the case of a transiting planet, is extremely challenging today. Our result even goes one step further, first by using this time real data, secondly by offering a broader range of methods tested. We first discuss below the results obtain on each system individually ordered by number of observations.

Results on C800

This system was certainly the most difficult with semi-amplitudes as small as 20 cm/s. As visible in Figure 8, none of the participants reached the 25% level for the transiting planets. On the other hand, most participants correctly recovered the amplitude of the giants, despite the existing cross-talk between the 11-year gaseous planet and the solar magnetic cycle. The phase of the 2 $M_{\mathrm{Jup}}$ signal was slightly anti-correlated with the solar magnetic cycle ($R=-0.61$) indicating some risk of interaction between the signals that was correctly accounted by the participants.

The performance on the transiting HZ planet is relatively poor, with a best score of only $\sim30\%$. However, even the benchmark analysis does not exceed $\sim45\%$ for this signal¹⁹, indicating an intrinsic difficulty in modelling the data given the available information.

This low performance is likely driven by the proximity of the signal to a one-year period, which results in poor phase coverage. This effect is exacerbated by the observational window of C800, which spans only about 7 months, leaving approximately $\sim40\%$ of the orbital phase unobserved.

A more detailed analysis shows that one of the two quadratures is entirely missing. These phases are particularly informative for constraining the semi-amplitude $K$, and their absence significantly degrades the ability to recover the signal accurately.

The only bad results are obtained for the 64-day transiting signal for which the average participant score is around 40% while the benchmark reaches 20% with a global overestimation of the mass. This likely shows that with realistic data, even 800 measurements

(21 nights) are not enough to recover such small amplitudes signals even when phases and periods are known.

As for C800 (see previous § 3.2.1), the worst result is obtained for the 373-day transiting HZ signal with an average score around 25%, while the benchmark grazes the 10% limit. Opposite to the signal in C800, this signal seems better constrained by both the benchmark and the participants. While the presence of 250 new observations (5 extra nights of telescope time) in A1050 compared to C800 should certainly help, it is doubtful that such an improvement can only be obtained by more data and part of the explanation can also be found in the phase coverage of the signal. Indeed, opposite to the signal in C800, the 8 months season length is now placing the seasonal gap no more on the quadrature of the signal but on the rising section of the phase. The presence of both quadratures certainly helping in the determination of the $K$ semi-amplitude.

Results on A1050

The results obtained on this data set are generally speaking the best as displayed in Figure 6. The levels of 25% (and even 10%) is reached several times by the participants. This is notably the case for the 26-day period planet close to the solar rotation period. While it should be reminded that most of our HARPS-N dataset is observing a quiet Sun from among the quietest cycles history (Jiang et al. 2015; Pesnell 2016), this seems to confirm than, nonetheless, 50 cm/s can be relatively well constrained in RVs with an HARPS-like stability precision instrument.

Despite the good average score value, this 373-day signal is also

the one for which the score of the participants is most fluctuant, with score spanning between 15% to 40%.

Results on B1300

This data set is the least successful. All the participants performed significantly worse than the benchmark, in particular for the mono-transit HZ signal at 474 days, the inner most planet at 41 days, and the almost non-detection by all of the 187 days signal. This system definitely contained more complexity than expected.

This result was not anticipated since this data set is the only one for which all the signals are correctly covered in phase, contains the largest number of observations ($N = 1300$ or 33 nights of telescope time) and has the longest season (8 months), so understanding why the signals were not recovered is a key point.

It should be recalled that this HZ planet was the unique mono-transit test performed, meaning that no tight prior on the period was used by the participants, even if this signal was also the most massive transiting HZ of the challenge ($7M_{\oplus}$) and correctly identified by most participants. As an interesting observation, the partial submission participant has attributed the mono-transit to the wrong planets (the 128-d signal) which can be explained by the lower boundary of the period prior spanning down to 100 days.

Furthermore, in this system, almost all the participants used a wrong planetary system model by missing the injected signal at 184 days of 34 cm/s and included a false detection at 630 days. It is very likely that such a wrong planetary model combined by the broad period prior is at the origin of the poor performance achieved on this data set.

By analysing the system architecture, we can notice a few peculiarities. The first is that three of the signals are close to the 1-year harmonics at 183, 122 and 91 days respectively. While not enough to explain the discrepancy (in particular since no eccentric signals were involved in this challenge), we can imagine that some 1-year systematics interact in some unpredictable way with the planetary signals. Such cross-talk were already mentioned for the previous data set except than this cross-talk is now occurring on a non-transiting planet.

We can question if the reduction by YARARA producing the YVA data set is not affecting those 1-year harmonics signals. This observation is partially motivated by the fact that the only participant that detected the 187-day signal used YV0 and not YVA, but we note that other participants using YV0 also missed this signal though. Furthermore, this explanation would not be coherent with the fact that no significant difference was observed between YV0 and YVA for the system A and C that were both containing 1-year planets.

Global comments

There are several culprits that may explain why no method is today yet able to reach the required precision. For some of the transiting signals, the difficulty may arise due to their proximity with 1-year that is an expected very difficult case scenario due to aliasing, poor phase coverage of the signal with season gaps and systematics all producing strong interaction around the 1-year period. Such an issue mainly arise for the HZ of solar-type stars, which can be mitigated when observing cooler K-dwarf for which the HZ moves inwards (see Sect. 3.5).

Part of the long-period difficulty may be due to imperfect stellar activity correction, even if we tested a rotational period transiting planet scenario. Similarly, supergranulation that is expected to produce part of the jitter observed in the HARPS-N solar data set (Moulla et al. 2023) may not leak power on such long timescale except if the properties of the supergranulation may also evolve on long timescale (O’Sullivan et al. 2025).

While other exotic stellar activity phenomena may be invoked (such as meridional flows, see e.g. Meunier & Lagrange (2020a)), the most likely culprit is the instrument stability itself, as well as granulation and supergranulation. It should be noted that while short-term stability of HARPS-N is known to be better than 50 cm/s (Dumusque et al. 2026), the long-term stability of instruments is very difficult to assess. The result of this challenge may indicate that the long-term stability of HARPS-like instruments is not sufficient to achieve the PLATO requirement for such small signals and more stable spectrograph (such as ESPRESSO) have certainly to be considered or further improvements should be achieved in this domain.

Given the result of the participants, one could infer the level of white noise which would yield the same $1\sigma$ uncertainties on the orbital elements (see Table A1). To compute this equivalent white noise, we use the analytical formula giving covariance of a least square fit. It has been tested that for the systems at hand, with over 1000 points, this analytical formula is consistent with a full MCMC computation of uncertainties. For system $A$ and $B$ and the time series generated with RV, we find that the corresponding level of white noise is 80 cm/s, in accordance with the simulation, and the YARARA YV2 method yields the same equivalent white noise value, while maintaining a comparable accuracy. In both cases, the amplitude is consistant with the combined effect of p-modes, granulation and super-granulation, the latter one being dominating. To go further, it is therefore extremely important to focus on these signals.

Planets

Injected

Exact Kep. (ARVE)

S-INDEX LINEAR (YVA)

Multi GPs (YVA)

SHELL CNN (YVA)

CCF SCALPELS (YVA)

CCF SCALPELS (YV0)

CCF LINSTAR (YVA)

S-INDEX LINEAR (YVA)

S-INDEX LINEAR (YV0)

YARARA YV2 (YVA)

1D GP (YVA)

Participant average

A26d

53

$55 \pm3$

$54 \pm4$

$57 \pm4$

$68 \pm2$

$54 \pm3$

$58 \pm6$

$58 \pm4$

$56 \pm4$

$54 \pm5$

$57 \pm3$

$56 \pm6$

A54d

50

$53 \pm3$

$52 \pm4$

$51 \pm4$

$52 \pm2$

$54 \pm4$

$52 \pm6$

$54 \pm4$

$54 \pm4$

$54 \pm5$

$52 \pm3$

$52 \pm5$

A106d

188

$195 \pm3$

$196 \pm5$

$194 \pm4$

$204 \pm2$

$193 \pm5$

$195 \pm4$

$196 \pm4$

$197 \pm5$

$197 \pm3$

$196 \pm5$

A373d

30

$30 \pm3$

$33 \pm4$

$38 \pm4$

$38 \pm2$

$42 \pm3$

$39 \pm5$

$40 \pm4$

$33 \pm4$

$34 \pm5$

$35 \pm3$

$36 \pm7$

A627d

182

$183 \pm4$

$193 \pm4$

$191 \pm4$

$223 \pm2$

$196 \pm6$

$224 \pm9$

$195 \pm5$

$189 \pm5$

$187 \pm5$

$197 \pm4$

$193 \pm8$

A922d

749

$757 \pm4$

$775 \pm4$

$770 \pm4$

$823 \pm2$

$768 \pm4$

$767 \pm8$

$772 \pm7$

$769 \pm5$

$767 \pm5$

$765 \pm4$

$770 \pm7$

FP1

$41$ (2185d)

$29$ (222d)

$250$ (3465d)

$52$ (2626d)

B41d

38

$44 \pm3$

$48 \pm5$

$45 \pm4$

$45 \pm2$

$48 \pm4$

$49 \pm4$

$46 \pm5$

$48 \pm5$

$49 \pm5$

$44 \pm3$

$49 \pm11$

$46 \pm4$

B52d

71

$72 \pm3$

$75 \pm4$

$77 \pm4$

$74 \pm2$

$74 \pm4$

$74 \pm4$

$73 \pm5$

$75 \pm5$

$75 \pm5$

$75 \pm3$

$74 \pm8$

$74 \pm4$

B91d

50

$57 \pm4$

$51 \pm5$

$50 \pm4$

$52 \pm2$

$54 \pm5$

$50 \pm5$

$51 \pm5$

$50 \pm5$

$51 \pm5$

$50 \pm3$

$53 \pm10$

$51 \pm5$

B128d

50

$51 \pm4$

$51 \pm4$

$57 \pm5$

$54 \pm2$

$52 \pm3$

$60 \pm5$

$51 \pm4$

$51 \pm5$

$52 \pm5$

$53 \pm3$

$56 \pm66$

$51 \pm6$

B187d

34

$42 \pm4$

$44 \pm5$

$44 \pm5$

B474d

56

$57 \pm4$

$80 \pm5$

$80 \pm2$

$80 \pm7$

$80 \pm5$

$80 \pm5$

$72 \pm5$

$83 \pm4$

$80 \pm5$

FP1

$44$ (620d)

$25$ (14d)

$42$ (627d)

$42$ (612d)

$82$ (642d)

$45$ (621d)

$45$ (624d)

$58$ (626d)

$30$ (630d)

$43$ (1073)

FP2

$35$ (586d)

$28$ (1381d)

$212$ (2528d)

$24$ (166d)

FP3

$79$ (1498d)

$301$ (3768d)

FP4

$77$ (3860d)

C64d

22

$23 \pm4$

$28 \pm7$

$32 \pm5$

$39 \pm2$

$33 \pm6$

$33 \pm6$

$25 \pm6$

$33 \pm6$

$31 \pm6$

$32 \pm4$

$32 \pm7$

C370d

19

$29 \pm6$

$20 \pm10$

$26 \pm8$

$20 \pm8$

$32 \pm6$

$29 \pm5$

$20 \pm8$

$34 \pm5$

$22 \pm6$

$25 \pm10$

C616d

1714

$1718 \pm4$

$1754 \pm11$

$1755 \pm5$

$1919 \pm2$

$1742 \pm8$

$1750 \pm7$

$1744 \pm14$

$1755 \pm6$

$1759 \pm5$

$1733 \pm5$

$1752 \pm13$

C3286d

2063

$2070 \pm6$

$2114 \pm26$

$2213 \pm34$

$2139 \pm20$

$2142 \pm16$

$2115 \pm27$

$2134 \pm13$

$2120 \pm12$

$2081 \pm6$

$2107 \pm86$

FP1

$27$ (264d)

$44$ (190d)

$43$ (1090d)

$64$ (1000d)

FP2

$72$ (733d)

However, a fundamental information to bear in mind is also provided by the benchmark that also struggle sometime to reach a good score as for the 370-day signal on C800 or 187-day signal in B1300. This seems to indicate that the amount of information in the data was simply insufficient to reach a better level than 25%. Recalling that such a benchmark was obtained using the ARVE solar red noise time-series without rotational modulation or instrumental stability consideration, this indicates that supergranulation alone when combined with suboptimal sampling can forbid such precision level.

We display in Figure 9 the maximum error $\varepsilon_K$ obtained by all the participants for each injected planet. The planets (columns) are ordered by the maximum error $\varepsilon_K$ of the benchmark. The participants (rows) are ordered according to their average $\varepsilon_K$. By carefully scrutinizing this figure, it becomes clear that the maximum error of the participants is in fact closely matching the benchmark even if slightly worse than it where the only exception is the dataset B1300.

Ranking of the methods

There exists no good and unbiased ranking metric to order the participants. First because the ranking may depend on the type of signals, secondly because the ranking metric can be biased by systematic bias in the methodologies used by the participants (see below). We started by categorizing the injected signals into distinct classes. Four classes were defined to evaluate participant performance under different conditions: (1) terrestrial signals (2) transiting signals, (3) habitable-zone transiting signals, and (4) all the signals.

As an example of it, we computed for each class the average maximum error $\varepsilon_K$ of all participants²⁰ and present the results in Figure 10. Participants are ranked according to their maximum error evaluated over the "terrestrial signals" category (top panel), and this ordering is preserved across the remaining panels for consistency and ease of comparison. While the HZ transiting planets constitutes the primary focus of this challenge, it can not be used due to the low statistic number ($N = 3$). The second best category of reference should have been the transiting planets, but due to the partial submission of one group, this metric is slightly biased for them and we therefore stuck to the terrestrial category for the reference. This choice is without consequence since all the transiting planets are terrestrial and as visible the ordering remains essentially the same using this group or the transiting one.

As expected, the best performance is achieved by the benchmark participant described in Sect. 2.3.1. With a few exceptions, the ranking of the participants is largely insensitive to the signal category considered. The average maximum error $\varepsilon_K$ is slightly higher when all signals are included, since this subset contains the large gaseous planets (see Figure 1), which are, by construction, the easiest signals to characterize.

No participant reaches a 10% maximum error for the transiting signals, and three methods fall below 25%: (1) YARARA V2, (2) Multi-GPs (Malavolta), and (3) Linear S-index (Damasso). It is worth noting that even the ARVE benchmark does not reach the 10% level, indicating that such performance was simply not achievable for the injected signals given the number of observations and the adopted realistic noise level.

However, drawing conclusions solely from this global metric would be misleading. In Appendix C, we also present the rankings based separately on precision $\sigma_K$ (Figure C1) and accuracy $\Delta K/K$ (Figure C2). The precision metric is particularly useful for diagnosing methods that underestimate their uncertainties. Such cases can be identified when a method reports a precision better than the benchmark, which effectively represents a near best-case scenario.

Two methods exhibit a better precision than the benchmark analysis (1) Shell CNN and (2) YARARA V2. This is suspicious, and likely indicates that the uncertainties are underestimated. A plausible explanation is that these methods do not simultaneously fit the Keplerian model and the red-noise component within a joint MCMC framework. Such rigid models might overlook the degeneracies between model components, and lead to underestimated uncertainties.

Systematic underestimation of uncertainties can bias the maximum error $\varepsilon_K$ ranking. Indeed, a method that overestimates its precision in a systematic way can achieve a better maximum error simply because the metric becomes dominated by its accuracy while the other participants continue to be evaluated on both metrics. To enable a fair comparison, it is therefore essential to complement the maximum error analysis with an evaluation of the number-of-sigma compatibility $N_\sigma$ for each participant. This analysis is presented in Figure 11.

As visible by the comparison with the previous Figure 10, while YARARA V2 was among the best scoring methods, it is also among the methods that systematically overestimate its precision with one of the largest $N_\sigma$ deviation around $N_\sigma= 1.5$ for terrestrial planets and even $N_\sigma\simeq3.0$ for HZ transiting ones.

In conclusion, we can not rely on a single one-dimensional ruler to decide of the success of a method and there always subsist a trade-off between accuracy and precision (or their alternative representation). We therefore display the average of each method in the 2D $(\Delta K/K, \sigma_K)$ and $(\varepsilon_K, N_\sigma)$ in the Figure 12. When analysing the methods that have both the smallest maximum error and are the most closely compatible to 1-sigma two methods can be highlighted that are the linear S-index detrending from Sect. 2.3.3 and the multi-GP from Sect. 2.3.4. Very interestingly, we see that the "best" participants themselves may depend if the metric used are the $(\varepsilon_K, N_\sigma)$ from the top row or $(\Delta K/K, \sigma_K)$ from the bottom row.

In definitive, we do not clearly attest the existence of an outstanding

method to the day and all methods seems to compete perfectly on such a data set. This is not extremely surprising given that some method may reveal their plain potential when working on more active stars, or on different SNR datasets.

Assessing the current RV precision in the RV field

Because this challenge involved a sufficiently large number of participants, employing a wide variety of methodologies — among the most innovative and representative currently available — we attempted to assess the present-day precision of the RV method.

The relatively small number of injected planetary signals is, of course, a limitation and should be kept in mind when interpreting the results, as should the caveats raised earlier regarding the extrapolation of results from the solar data set (see Sect. 2.2.2).

It is becoming increasingly standard practice not to rely on a single pipeline or analysis when publishing planetary signals, but instead to run several independent pipelines in parallel and compare their outcomes. This healthy practice is now routine for JWST, where independent pipelines are executed in parallel during data processing (e.g Kirk et al. 2024), and we advocate that a similar standard should progressively become the norm in RV studies.

However, combining results from different reductions or analyses remains a matter of debate. One cannot simply average the outputs of multiple pipelines, since they typically start from (or partially share) the same input data. This makes posterior propagation and model averaging non-trivial.

We therefore adopted a simpler approach: for each injected signal, we selected the best score achieved by any participant. Naturally, this choice introduces a bias and slightly overestimates the final maximum error $\varepsilon_k$. As a counterbalance, however, we treated all false positive detections as valid detections, even if reported by only a single participant. By doing so, a statistic of 27 false detections is obtained by all the participants combined. The results of this analysis are shown in Figure 13, and reveal an interesting trend.

Two distinct regions clearly emerge: one where detection scores are high and false positives are essentially absent, and another where scores are low and false detections²¹ become likely. Although the number of injected planets was insufficient to densely populate the $(K,P)$ parameter space, a rough empirical boundary $K_{\mathrm{crit}}(P)$ can nevertheless be identified as follows for an HARPS-like instrument:

The number of false positives emerging from this challenge is a major concern, although they primarily occur — as expected — at low amplitudes and long orbital periods. This result may be concerning since a significant number of planets have been published in this region of the parameter space over the past years, Among the systems falling into the suspicious region are for instance: HD10700 f, HD20794 f, HD134060 c, HD95735 c, HD34445 g, GJ433 c and GJ15A c.

Although we did not perform a case-by-case reassessment of those published detections, for instance to verify whether they were supported by independent methods, it is important to emphasize that our challenge was conducted under highly favorable conditions. The data sets typically contains an average of $N \sim1000$ observations and were based on one of the most stable spectrographs currently available, operating under conditions and SNR that outperform most real stellar RV data sets obtained to date.

Given these considerations, our results suggest that caution is warranted when interpreting low-amplitude, long-period RV detections. In particular, we recommend that planets previously labeled as “confirmed” based solely on RV measurements, especially in this region of parameter space, should be reconsidered back as “candidates” unless supported by independent lines of evidence.

This further strengthens the case for the use of multiple independent pipelines and reduction strategies when assessing RV detections. Cross-validation between methods should become standard practice. However, our results also demonstrate that such sanity checks are not, by themselves, a guarantee against false positives.

In particular, among the false detections identified in the challenge, we find a $\sim600$-day Super-Earth candidate ($\sim5M_{\oplus}$) in data set B1300 that was independently detected by all participants. This illustrates that methodological diversity alone does not fully protect against systematic biases or activity-induced signals that can convincingly mimic planetary signatures.

K-dwarf stars: the RV sweet spot?

In this section, we present a simple scaling analysis of the signals from the RV challenge to assess the impact of spectral type on the PLATO requirements.

Such an exercise should, however, be interpreted with caution. Changing the spectral type affects not only the stellar mass and the location of the HZ, but also the stellar activity budget and the spectral signatures. As a result, it is not straightforward to predict whether the various methods used by the participants would perform better or worse for different spectral types. Moving towards cooler spectral types nevertheless provides several advantages.

First, cooler stars naturally yield larger RV semi-amplitudes for a given planetary mass. In the HZ, the semi-amplitude scales as

reflecting both the stellar-mass dependence and the inward shift of the HZ.

Second, the inward shift of the HZ leads to shorter orbital periods, which helps mitigate observational window effects linked to Earth’s orbit. This improves phase coverage and therefore enhances both detectability and characterization, as signals can be sampled over multiple orbital cycles.

Finally, the expected RV variability induced by convective flows is reduced. The amplitudes of p-modes, granulation (and supergranulation) are predicted to decrease by about 30% (Meunier et al. 2019b). In addition, cooler stars tend to rotate more slowly and exhibit lower spot contrast. As a consequence, stellar activity signals are expected to be more strongly dominated by the inhibition of convective blueshift (ICB) described in Meunier et al. (2010), rather than by flux variations induced by active-region contrast (Meunier & Lagrange 2013). The ICB amplitude being itself smaller since the convective flows are slower for cooler type-stars (Liebing et al. 2021; Al Moulla et al. 2024).

We illustrate the expected gain in Figure 14. We adopt the same KP detection limit as the red curve shown in Figure 12 and defined in Eq. (7). This limit is then rescaled with stellar mass for two additional spectral types: K1 ($M_{\star} = 0.85\,M_{\odot}$) and K5 ($M_{\star} = 0.70\,M_{\odot}$). The HZ is correspondingly shifted inwards following Kopparapu et al. (2013, 2014). The resulting HZ period ranges are [272, 800] days for a G2 star, [136, 472] days for a K1 star, and [65, 237] days for a K5 star.

We apply the same scaling to all planetary systems in the challenge by keeping the planetary masses fixed while rescaling their orbital periods such that HZ transiting planets receive the same incident flux as the Earth–Sun system. This analysis reveal that while for a G2 type star the requirement are achieved for $\sim8\,M_{\oplus}$ planets for the innermost part of the HZ, the PLATO 10% requirements may be fulfill down to $\sim3\,M_{\oplus}$ for a K1 and $\sim2\,M_{\oplus}$ for a K5 star.

We note that similar scaling arguments have previously been used to motivate the search for HZ planets around M dwarfs, although several important limitations were not anticipated. In particular, M dwarf spectra in the visible are often heavily blended and obtained at relatively low signal-to-noise ratio, which has driven the development of infrared spectrographs. However, the technical challenges associated with infrared observations were initially underestimated (Plavchan et al. 2015; Fischer et al. 2016), and their instrumental stability has only reach very recently the level achieved in the visible (Artigau et al. 2024; Bouchy et al. 2025; Varas et al. 2025).

Furthermore, M dwarfs typically host stronger magnetic fields, leading to enhanced stellar activity signals, in particular through Zeeman broadening (Reiners et al. 2013; Kochukhov 2021). Finally, the question of habitability itself has been rapidly called into question: the habitable zone lies very close to the host star, and young M dwarfs undergo intense flaring phases that may be sufficient to erode atmospheres or prevent the emergence of life.

Conclusions

In this paper, we presented an RV data challenge organized within the PLATO Stellar Variability Mitigation Working Group to compare different analysis methods and establish a robust analysis strategy in the context of the future PLATO follow-up, in particular regarding the precision and accuracy achievable on RV semi-amplitudes.

The solar HARPS-N dataset was used to generate three synthetic datasets designed to closely mimic realistic follow-up observations. Three consecutive 5-minute exposures were combined to produce an effective integration time of 15 minutes, yielding datasets with $S/N_{550} \sim540$. This signal-to-noise ratio corresponds to observations of a star with magnitude $m_V \approx5.4$ using a 3.6m telescope such as HARPS, and $m_V \sim7.2$ with ESPRESSO on the VLT.

Planetary signals were injected directly into the spectra, after which CCFs and RV time series were derived. Particular care was taken to ensure that the injected signals affected only the stellar component of the spectra, and not systematics tied to the Earth laboratory rest frame.

Three multi-planet systems were simulated, each containing between four and six planets, including one transiting habitable-zone planet as well as additional non-transiting companions. To limit the dimensionality of the problem, all planetary orbits were assumed to be circular, such that only the semi-amplitudes needed to be fitted for the known transiting planets. While this assumption simplifies the analysis, the impact of eccentricity on signal recovery could be explored in future studies.

Four metrics were used to evaluate participant performance based on the recovered K semi-amplitudes: (1) accuracy $\Delta K/K$, (2) precision $\sigma_K$, (3) maximum error $\epsilon_K$, and (4) significance in units of sigma $N_\sigma$. No overall ranking was established, as performance was found to depend strongly on both the chosen metric and the specific signals. Moreover, the limited number of injected signals prevents any statistically robust ranking.

For terrestrial-mass signals, no participant achieved a precision better than 10%, and even the benchmark time series reached only ∼ 15%. This indicates that supergranulation alone is sufficient to prevent reaching the PLATO requirement at this level. Most participants exhibit maximum error $\varepsilon_k$ in the range of 20–25%.

By analysing both false positives and best-performing detections, a region of the $(K, P)$ parameter space was identified where planetary detections should be treated with caution for instruments with HARPS-like stability. Our results highlight the importance of comparing independent analysis methods when claiming low-amplitude planetary signals. In this regime, signals that are not consistently recovered by different methods should currently be considered as planet candidates rather than confirmed RV detections until corroborated by additional evidence.

Additional challenges arise from observational sampling, particularly for signals near a one-year period, which complicates the detection of HZ planets around solar-type (G2) stars. In contrast, for a similar magnitude, K dwarfs benefit from shorter HZ periods. Simple scaling relations suggest that the 10% precision requirement could be met for HZ planets of ∼ 8–10, $M_{\oplus}$ around G2 stars, and down to ∼ 2–5, $M_{\oplus}$ around K dwarfs. These results further strengthen the case for K dwarfs as the most favourable targets for the first robust detections of habitable-zone planets with the RV technique. However, we noticed that this difficulty was enhanced in the RV challenge due to the differential extinction of the Sun producing season length of the three data sets around $S \sim220$ days, while the visibility of the PLATO south field of view from Paranal is $S \sim280$ days.

This challenge provides three new benchmark data sets publicly available for the RV community to test and compare analysis methods. It should be noted, however, that future participants will be aware of the injected systems, whereas this challenge was conducted under near-blind conditions.

Future extensions of this work include investigating the impact of lower SNR, eccentric signals and incorporating data from additional instruments such as ESPRESSO. The priority seems, to us, to improve the correction of granulation and supergranulation. The latter appears to be one of the dominant limitations for future high-precision RV measurements. A dedicated data challenge focusing specifically on this component would be a natural next step.

Some improvements on the planetary mass determination can be expected for the PLATO mission in particular if TTVs are used jointly to the Keplerian fit in RVs (Leleu et al. 2024), even if the improvement is yet unclear in presence of non-transiting planets. Best chance of improvements are expected if the supergranulation is better understood in the coming decade since this component alone is sufficient to make the 10% difficult to achieve, in line with the results of Meunier & Lagrange (2020b). Such advancements are expected either thanks to PoET, the solar telescope of ESPRESSO (Santos et al. 2025) or intensive future RV survey such as the second Earth project 2ES (Stürmer et al. 2024) or the Terra Hunting Experiment (Thompson et al. 2016; Hall et al. 2018) both entering in operation this year and having supergranulation as part of their open scientific question.

In addition, the long-term stability of high-precision spectrographs should remain a priority, as instrumental interventions and ageing may introduce systematics that become critical on the multi-year timescales required to detect habitable-zone planets.

Acknowledgements

MC and SA acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 865624). JBP, IR and MP were funded by Spanish grants PID2021-125627OB-C31, PID2023-150491NB-I00, and PID2024-158486OB-C31, and by “ERDF A way of making Europe”, by the programme Unidad de Excelencia María de Maeztu CEX2020-001058-M financed by MCIN/AEI/10.13039/501100011033 and by the MaX-CSIC Excellence Award MaX4-SOMMA-ICE, by the Generalitat de Catalunya/CERCA programme, and by the European Research Council (ERC) under the European Union’s Horizon Europe programme (ERC Advanced Grant SPOTLESS; no. 101140786). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. A.M. acknowledges funding from a UKRI Future Leader Fellowship, grant number MR/X033244/1 and a UK Science and Technology Facilities Council (STFC) small grant ST/Y002334/1. STH acknowledges support from UKRI-STFC grants ST/X001571/1 and UKRI974. ACC acknowledges support from STFC consolidated grant number ST/V000861/1, UKRI/ERC Synergy Grant EP/Z000181/1 (REVEAL), and UKSA grant number ST/X002217/1. NCS is funded by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This work was supported by Fundação para a Ciência e a Tecnologia (FCT) through national funds under the research grant UID/04434/2025 (DOI 10.54499/UID/04434/2025). The French researchers acknowledge support on PLATO by the Centre national d’études spatiales (CNES). TGW acknowledges support from the University of Warwick and UKSA. D.N. acknowledges support from PLATO ASI-INAF agreements no. 2022-28-HH.0. KA acknowledges support from the Swiss National Science Foundation (SNSF) under the Postdoc Mobility grant P500PT_230225. NCS, KA, and EG acknowledge funding by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This work was supported by Fundação para a Ciência e a Tecnologia (FCT) through national funds under the research grant UID/04434/2025 (DOI 10.54499/UID/04434/2025). C.A.W. would like to acknowledge support from the UK Science and Technology Facilities Council (STFC, grant number ST/X00094X/1).

This research has made use of NASA’s Astrophysics Data System (ADS) bibliographic services. We acknowledge the community efforts devoted to the development of the following open-source packages that were used in this work: numpy (numpy.org), matplotlib (matplotlib.org), and astropy (astropy.org). This publication makes use of The Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva (CH) dedicated to extrasolar planets data visualisation, exchange and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. The DACE platform is available at https://dace.unige.ch.

Data availability

The data underlying the RV challenge are publicly available in the Google Drive repository: https://drive.google.com/drive/u/1/folders/1Jbz-i2ID8I1Ir90w3AVSDzygctGuFnqo. Data are also publicly available on Zenodo (10.5281/zenodo.20151946).

Tables of participants output

Evaluation of the spectra RV shift injection on the CCFs moments

A major update of this RV challenge compared to previous efforts is the injection of planetary signals directly into the spectra. To implement this, we first evaluated the impact of the RV-shift injection algorithm in order to validate its use for the challenge.

In principle, the cleanest way to inject a planetary signal is to apply a Doppler shift to the wavelength solution prior to computing the CCFs using the standard extraction algorithm. However, this approach was not compatible with our framework, as the RV challenge relies on YARARA to process and clean the spectral time series.

YARARA operates on a spectral time-series matrix in which all spectra share a common wavelength grid. This uniform grid is a core assumption of the data format required by the code. Consequently, any injected RV shift must be followed by a reinterpolation step back onto the native wavelength grid.

While reinterpolation is sometimes regarded with caution in high-resolution spectroscopy, there is no strong evidence that it significantly degrades the data quality when applied carefully. In the worst case, interpolation may introduce a slight increase in the noise level, but this effect is typically negligible and does not materially affect the scientific interpretation.

In our implementation, the shifted spectra were reinterpolated onto the original grid using cubic splines. Naturally, the validity of such an interpolation depends on the amplitude of the Doppler shift. What constitutes a “small” shift must therefore be evaluated in the specific context of the challenge.

To quantify the impact of the method, we selected a representative spectrum from the data set and created an exact copy. Photon noise was then injected to reproduce the final RV precision of the solar data set used in the challenge, namely $\sigma_{\mathrm{RV}} \simeq15$ cm/s. We applied Doppler shifts ranging from 0 to 100 m/s, followed by reinterpolation, exactly as performed in the construction of the RV challenge.

The results are shown in Figure B1. We did not observe any significant deviation between the injected and recovered RV values over the full tested range up to 100 m/s. However, a weak linear trend is visible in the full width at half maximum (FWHM) of the CCF, indicating a small leakage of the Doppler shift into the line-profile shape.

We quantified this leakage by fitting a linear trend to the CCF moments. For the RV itself, we obtained:

This implies that for an injected shift of 100 m/s, the recovered RV is underestimated by 17 cm/s. This is a relative difference of 0.172% on the injected amplitude, far smaller than the 10% objective of this challenge. Such a bias is extremely small which demonstrate that the method selected here will not introduce any significant bias for the participants and confirm the validity of the method. However, given the precision requirements of the challenge, we therefore limited the Keplerian semi-amplitudes to 20 m/s, ensuring that the maximum bias remains below 3.4 cm/s.

Furthermore, restricting the injected amplitudes to this range also prevents significant contamination of the CCF moments, since over such RV ranges the leakage remains below the intrinsic jitter precision of the data set for all considered CCF diagnostics.

Ranking of the participants with precision and accuracy

While an attempt was made to develop a single one-dimensional ruler to rank the participants, for instance by the use of the maximum error $\epsilon_k$ metric that combines both the accuracy and precision we realise that all the score metric are imperfect to some degree and still require a 2D parameters space to analyse the results.

We display the ranking of all the participants based on each individual metric, respectively the accuracy in Figure C2 and the precision in Figure C1.

This paper has been typeset from a \TeX/\LaTeX{} file prepared by the author.