# The SPIRou Legacy Survey: Detection of a nearby world orbiting in the habitable zone of Gl 725 B achieved by correcting strong telluric contamination in near-infrared radial velocities with wapiti

M. Ould-Elhkim, C. Moutou, J.-F. Donati, P. Cortés-Zuleta, X. Delfosse, É. Artigau, C. Cadieux, P. Charpentier, A. Carmona, I. Boisse, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

M. Ould-Elhkim, C. Moutou, J.-F. Donati, P. Cortés-Zuleta, X. Delfosse, et al. The SPIRou Legacy Survey: Detection of a nearby world orbiting in the habitable zone of Gl 725 B achieved by correcting strong telluric contamination in near-infrared radial velocities with wapiti. Astronomy & Astrophysics - A&A, 2026, 705, pp.A234. 10.1051/0004-6361/202555469. hal-05490007

HAL Id: hal-05490007

https://hal.science/hal-05490007v1

Submitted on 26 Feb 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# The SPIRou Legacy Survey

## Detection of a nearby world orbiting in the habitable zone of Gl 725 B achieved by correcting strong telluric contamination in near-infrared radial velocities with wapiti

M. Ould-Elhkim<fnref n="1" />,*, C. Moutou<fnref n="1" />, D., J.-F. Donati<fnref n="1" />, P. Cortés-Zuleta<fnref n="2" />, X. Delfosse<fnref n="3" />, É. Artigau<fnref n="4" />,<fnref n="5" />, C. Cadieux<fnref n="4" />, P. Charpentier<fnref n="1" />, A. Carmona<fnref n="3" />, I. Boisse<fnref n="6" />, C. Reylé<fnref n="7" />, E. Gaidos<fnref n="8" />, R. Cloutier<fnref n="9" />,<fnref n="10" />, G. Hébrard<fnref n="11" />, L. Arnold<fnref n="12" />, J.-D. do Nascimento, Jr.<fnref n="13" />, N. J. Cook<fnref n="4" />, and R. Doyon<fnref n="4" />

Received 9 May 2025 / Accepted 17 November 2025

## ABSTRACT

**Context.** M dwarfs are prime targets in the search for exoplanets because of their prevalence and because low-mass planets can be better detected with radial velocity (RV) methods. In particular, the near-infrared (NIR) spectral domain offers an increased RV sensitivity and potentially reduced stellar activity signals. Precise NIR RV measurements are strongly affected by telluric absorption lines from the Earth's atmosphere, however.

**Aims.** We searched for planets orbiting Gl 725 B, a nearby late-M dwarf located at 3.5 pc, using high-precision SPIRou RV observations. We also assessed the effect of telluric contamination on these measurements and evaluated the performance of the weighted principal component analysis reconstruction (wapiti) method, which is a weighted principal component analysis (wPCA) approach designed to mitigate these systematics and to improve the sensitivity of planet detections.

**Methods.** Using synthetic and observational SPIRou data, we simulated the effect of telluric lines on RV data under varying barycentric Earth radial velocity (BERV) conditions. We then applied the wapiti method for identifying and correcting telluric-induced systematics in line-by-line RVs. The method was tested through an injection-recovery test on simulated data and was subsequently applied to real SPIRou observations of Gl 725 B.

**Results.** wapiti successfully corrects telluric contamination in simulated and real datasets. This enhances the detectability and accuracy of planetary signals. In the corrected Gl 725 B dataset, we identified a two-planet system composed of a candidate inner planet (Gl 725 Bb), with periods of 4.765 ± 0.004 days and a semi-amplitude of 1.4 ± 0.3 m s⁻<fnref n="1" />, and a confirmed planet Gl 725 Bc, with a period of 37.90 ± 0.17 days and a semi-amplitude of 1.7 ± 0.3 m s⁻<fnref n="1" />. Their minimum mass is 1.5 ± 0.4 M⊕ and of 3.5 ± 0.7 M⊕, respectively, and the outer planet is located in the habitable zone of its host star. Using a multi-dimensional Gaussian process framework to model and correct for stellar activity, we also recovered a stellar rotation period of 105.1 ± 3.3 days.

**Key words.** methods: data analysis – techniques: radial velocities – techniques: spectroscopic – stars: low-mass – planetary systems – stars: individual: Gl 725 B

## 1. Introduction

High-precision radial velocity (RV) is a powerful technique for detecting and characterizing exoplanets by measuring their

masses and orbital parameters. The accuracy of RV measurements is affected by photon noise, instrumental stability, telluric contamination, and stellar activity (Plavchan et al. 2015; Fischer et al. 2016). Stellar surface inhomogeneities such as spots and plages can introduce RV signals of a few m s⁻<fnref n="1" />

<footnote n="1">¹ Université de Toulouse, UPS-OMP, IRAP, CNRS, 14 avenue E. Belin, Toulouse 31400, France</footnote>

<footnote n="2">² SUPA School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews KY16 9SS, UK</footnote>

<footnote n="3">³ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="4">⁴ Trottier Institute for Research on Exoplanets, Université de Montréal, Département de Physique, C.P. 6128 Succ. Centre-ville, Montréal, QC H3C 3J7, Canada</footnote>

<footnote n="5">⁵ Observatoire du Mont-Mégantic, Université de Montréal, Département de Physique, C.P. 6128 Succ. Centre-ville, Montréal, QC H3C 3J7, Canada</footnote>

<footnote n="6">⁶ Aix Marseille Univ, CNRS, CNES, Institut Origines, LAM, Marseille, France</footnote>

<footnote n="7">⁷ Université Marie et Louis Pasteur, CNRS, Observatoire des Sciences de l'Univers THETA Franche-Comté Bourgogne, Institut UTINAM (UMR 6213), 25000 Besançon, France</footnote>

<footnote n="8">⁸ Department of Earth Sciences, University of Hawaii at Manoa, Honolulu, HI 96822, USA</footnote>

<footnote n="9">⁹ Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="10">¹⁰ Department of Physics & Astronomy, McMaster University, 1280 Main Street West, Hamilton, ON L8S 4L8, Canada</footnote>

<footnote n="11">¹¹ Institut d'astrophysique de Paris, UMR7095 CNRS, Université Pierre & Marie Curie, 98bis boulevard Arago, 75014 Paris, France</footnote>

<footnote n="12">¹² Canada France Hawaii Telescope (CFHT) Corporation, UAR2208 CNRS-INSU, 65-1238 Mamalahoa Hwy, Kamuela 96743 HI, USA</footnote>

<footnote n="13">¹³ Universidade Federal do Rio Grande do Norte (UFRN), Departamento de Física, 59078-970, Natal, RN, Brazil</footnote>

<footnote>* Corresponding author: merwan.ould-elhkim@lam.fr</footnote>
(Dumusque et al. 2011; Dumusque 2018; Cegla et al. 2019). These activity-induced variations depend on the wavelength and are often reduced in the red, optical, and NIR (Marchwinski et al. 2015; Robertson et al. 2020; Cale et al. 2021), at least for the most active M dwarfs (Cortés-Zuleta et al. 2023). This makes multiwavelength RV studies important for distinguishing between stellar and planetary signals (Huélamo et al. 2008; Mahmud et al. 2011; Carmona et al. 2023; Cortés-Zuleta et al. 2023).

The interest in M dwarfs has grown with the realization that they host more rocky planets than solar-type stars (Bonfils et al. 2013; Dressing & Charbonneau 2015; Gaidos et al. 2016; Hsu et al. 2020; Sabotta et al. 2021; Mignon et al. 2025). Their lower masses enhance RV signals from small planets (Anglada-Escudé et al. 2016; Faria et al. 2022; Suárez Mascareño et al. 2023), and they are both numerous and close to us (Henry et al. 2006). Moreover, the flux contrast between active regions and the photosphere is typically lower than in Sun-like stars, which might reduce the RV jitter induced by stellar activity (Berdyugina 2005; Crockett et al. 2012; Carmona et al. 2023; Cortés-Zuleta et al. 2023), especially when this is combined with variations in line depths between spots and photospheres (Larue et al. 2025).

These favorable characteristics have motivated the development of a new generation of NIR spectrographs, such as GIANO (Carleo et al. 2016), HPF (Metcalf et al. 2019), SPIRou (Donati et al. 2020), CARMENES-NIR (Bauer et al. 2020), IRD (Hirano et al. 2020), iSHELL (Rayner et al. 2022), NIRPS (Thibault et al. 2022), and CRIRES+ (Dorn et al. 2023). It remains difficult to achieve high RV precision in the NIR, however, because of telluric absorption and emission lines from the Earth's atmosphere (Bean et al. 2010). These lines are mainly due to water vapor, methane, carbon dioxide, and oxygen, and they are more prominent in the NIR and do not share the stellar Doppler shift.

Telluric lines have traditionally been corrected for using observations of hot featureless stars (Artigau et al. 2014), and water line residuals as low as 2% were achieved in the 900–1350 nm range (Sameshima et al. 2018). More recent approaches relied on synthetic telluric models such as TAPAS (Bertaux et al. 2014), with best-case residuals at 1–2% (Ulmer-Moll et al. 2019).

The residual telluric contamination can still significantly affect the RV precision, especially for targets with poor coverage of the barycentric Earth radial velocity (BERV; i.e., a BERV span ≤10 km s⁻<fnref n="1" /> over a year). The BERV corresponds to the projection of the Earth's velocity vector onto the line of sight to the star, and its amplitude varies depending on the stellar declination and meridian transit. Poor BERV coverage limits the modulation of telluric lines in observations and makes it harder to distinguish them from true stellar signals. Simulations performed by Wang et al. (2022); Latouf et al. (2022) using synthetic spectra have shown that with optimized mitigation, the RV scatter from telluric residuals can be reduced to below 1.5 m s⁻<fnref n="1" /> in the NIR.

In this context, the case of Gl 725 B is particularly relevant. This nearby late-type M dwarf is a target of the SPIRou Legacy Survey Planet Search (SLS-PS<fnref n="1" />) program, which aims to detect and characterize planetary systems around ~50 M dwarfs using high-precision NIR RV measurements (Moutou et al. 2023). The source Gl 725 B has a narrow BERV coverage of less than 5 km s⁻<fnref n="1" />, however, which limits our ability to separate the effects of telluric contamination over time.

We simulate the effect of telluric absorption on SPIRou RV data specifically for Gl 725 B using the SPIRou Data Reduction Software APERO (Cook et al. 2022) and a line-by-line (LBL) RV extraction method (Dumusque 2018; Artigau et al. 2022).

<table><caption>Table 1. Stellar parameters of Gl 725 B.</caption><thead><tr><th>Parameter</th><th>Gl 725 B</th><th>Ref.</th></tr></thead><tbody><tr><td>T<sub>eff</sub> (K)</td><td>3379 ± 31</td><td>1</td></tr><tr><td>[M/H]</td><td>-0.28 ± 0.10</td><td>1</td></tr><tr><td>[α/Fe]</td><td>0.14 ± 0.04</td><td>1</td></tr><tr><td>log g</td><td>4.82 ± 0.06</td><td>1</td></tr><tr><td>Radius (R<sub>☉</sub>)</td><td>0.280 ± 0.005</td><td>1</td></tr><tr><td>Mass (M<sub>☉</sub>)</td><td>0.25 ± 0.02</td><td>1</td></tr><tr><td>log(L/L<sub>☉</sub>)</td><td>-2.038 ± 0.003</td><td>1</td></tr><tr><td>P<sub>rot</sub> (d)</td><td>135 ± 15</td><td>2</td></tr><tr><td>P<sub>rot</sub> (d)</td><td>105.1 ± 3.3</td><td>This work</td></tr><tr><td>Age (Gyr)</td><td>8.7 ± 2.1</td><td>This work</td></tr></tbody></table>

Notes. 1: Cristofari et al. (2022); 2: Donati et al. (2023b).

We then assess the capability of the *wapiti* (Weighted principal component analysis) reconstruction of the Gl 725 B system to mitigate telluric systematics in these simulations. Finally, we apply *wapiti* to real SPIRou observations of Gl 725 B, which leads to the detection of a two-planet system, including a candidate that is located within the habitable zone.

The structure of this paper is as follows. Section 2 introduces the M dwarf Gl 725 B and presents its SPIRou RV and TESS photometric data. Section 3 describes our simulation setup and the effect of limited BERV coverage. In Sect. 4 we perform an injection-recovery test to evaluate the performance of *wapiti*. Section 5 presents the analysis of SPIRou observations of Gl 725 B and reports the detection of a two-planet system. We summarize our conclusions in Sect. 6.

## 2. Observations of Gl 725 B

The M5.0V dwarf star Gl 725 B (Struve 2398 B, HD 173740) (Kirkpatrick et al. 1991) lies at a distance of 3.5231 ± 0.0003 pc from Earth, as reported by Gaia Collaboration (2020). A summary of its physical properties is available in Table 1. Additionally, the longitudinal magnetic field ($B_t$) of this star using circularly polarized spectropolarimetry was analyzed by Donati et al. (2023b) and led to the identification and characterization of a rotational period for Gl 725 B, estimated to be 135 ± 15 days.

The source Gl 725 B is part of a binary system. Its companion Gl 725 A is classified as an M3V dwarf star (Kirkpatrick et al. 1991), and its analysis with SPIRou and SOPHIE data led to the discovery of a super-Earth was reported by Cortés-Zuleta et al. (2025). These authors estimated the separation between these two stars as 63 ± 1 au, with an eccentricity of 0.29 ± 0.01 and an orbital period of 871 ± 108 years. A companion with a mass of 15 $M_J$ and an orbital period of 5.5 years was reported around Gl 725 B by Baize (1976), but this claim was not confirmed by the more recent analysis of Kervella et al. (2019). A potential planet with a period of 2.7 d and a 1.2 m s⁻<fnref n="1" /> semi-amplitude was also previously reported around Gl 725 B from high-cadence observations with HARPS-N in the optical, but was considered to be comparable to the noise floor (Berdiñas et al. 2016).

### 2.1. SPIRou data

The high-resolution NIR spectropolarimeter SPIRou is installed at the Cassegrain focus of the CFHT in 2018. It covers a spectral range from 0.98 to 2.35 µm with a resolving power of 70 000. The instrument was specifically designed for the detection and

<footnote n="1">¹ https://spirou-legacy.irap.omp.eu/doku.php/</footnote>
characterization of exoplanets around low-mass stars in the NIR domain, where M dwarfs are brightest (Donati et al. 2020).

Observations of Gl 725 B with SPIRou were conducted between February 14, 2019, and April 19, 2022, resulting in a total of 208 measurements. All spectra were processed with version 0.7.288 of the SPIRou data reduction software APERO (Cook et al. 2022). APERO includes dedicated pipelines for correcting telluric absorption and sky emission lines. It relies on a combination of TAPAS atmospheric models (Bertaux et al. 2014) and on a library of hot-star observations acquired under varying conditions (Artigau et al. 2014). These telluric-corrected spectra form the basis for the subsequent RV extraction.

The RVs were obtained using the Line-by-Line (LBL) method<fnref n="2" />, version 0.65.003. This technique, initially developed by Dumusque (2018) for optical data and extended to the NIR by Artigau et al. (2022), extracts RVs by tracking the displacement of individual spectral lines. This offers an improved sensitivity for M dwarfs compared to classical cross-correlation techniques (Artigau et al. 2022). It also generates time series for each spectral line individually, providing powerful diagnostics to isolate and characterize systematic effects.

Instrumental systematics were mitigated in two steps. First, a Fabry-Pérot (FP) etalon, used as a simultaneous reference, allowed us to track short-term instrumental drifts during the night; the measured drift was subtracted from the raw RVs to correct for intra-night velocity variations. Second, a nightly zero point (NZP) correction was applied to account for long-term instrumental drifts, such as those arising from temperature variations or aging of optical fibers. These drifts were monitored using observations of long-term stable stars, and the corresponding NZP correction was subtracted from the nightly RVs.

Three outliers were identified and removed using a median absolute deviation (MAD) criterion. After these corrections, the dataset achieved a median signal-to-noise ratio per 2.28 km s⁻¹ pixel bin at the center of the H band of 161.9, with a median airmass of 1.355.

The individual per-line RVs obtained by the LBL method were then analyzed using a weighted principal component analysis (wPCA; Delchambre 2015). This technique was recently applied to SPIRou (Ould-Elhkim et al. 2023) and HARPS (Cretignier et al. 2023) data and identifies and removes systematics such as residual telluric contamination. In the next section, we use simulated datasets based on Gl 725 B to assess the effect of residual tellurics on the RV precision and to demonstrate the efficiency of the wapiti correction method.

## 2.2. TESS data

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) has observed nearly the complete sky in its search for exoplanets around bright stars. Its favorable position in the sky allowed it to observe the binary system Gl 725 in Sectors 14–26, 40, 41, 47–57, 59 and 60, covering a period between July 2018 and January 2023. The data from Sector 14 were not considered in our analysis because the position of the star was too close to the edges of the detector. We obtained the presearch data-conditioning (PDC) flux processed by the TESS Science Processing Operations Center (SPOC). The aperture used to perform the flux extraction was not consistent across the sectors, however, meaning that in some sectors, the aperture does not include the two system components. Because the binary

components lie so closely to each other and their magnitudes are so similar, it is not possible to independently extract their fluxes. For the sake of consistency, we performed a custom aperture photometry extraction with an aperture that included the fluxes of Gl 725 A and B.

The resulting flux time series underwent a cleaning process that removed 3σ outliers and a detrending using the wōtan<fnref n="3" /> Python package (Hippke et al. 2019), which includes a method based on a time-windowed sliding filter with an iterative robust location estimator. The edges of the TESS time series are usually affected by strong systematics, and we therefore did not consider them in the detrending model. The standard deviation of the detrended light curves is 130 ppm.

# 3. Simulating the effect of telluric absorption lines on LBL data

## 3.1. Outline of the simulation method

The limited BERV coverage of Gl 725 B (BERVMAX = 4.6 km s⁻¹) means that telluric contamination leaves systematic residuals in its RV time series, even after correction. To evaluate and quantify this effect and to assess the efficiency of correction methods, we performed dedicated simulations centered on Gl 725 B.

Our approach required an accurate stellar template with minimum telluric residuals. We therefore used the M5.0V dwarf Gl 905 (HH And, Ross 248; Alonso-Floriano et al. 2015), which exhibits a large BERV excursion (BERVMAX = 23 km s⁻¹). Its extensive BERV range and high data quality make it suitable for constructing a clean stellar template $T_★$ with negligible residual telluric contamination.

At each epoch, we replaced the observed Gl 725 B spectrum with the Gl 905 template and added photon noise for that epoch. To simulate the effect of telluric contamination, we created a 2D (wavelength times order) telluric template $T_⊕$ from the output of APERO by taking the median (in the Earth frame) of the telluric spectra derived from Gl 905 observations (Cook et al. 2022).

To simulate realistic telluric contamination, we multiplied the noisy stellar template with $T_⊕$, modulated by a parameter $\alpha$ that encodes imperfect telluric correction,

$$T_{\text{contaminated}} = T(t) \times T_{\oplus}^{\alpha(t)}. \qquad (1)$$

The $\alpha$ parameter was drawn at each epoch from a normal distribution $N(\mu_\alpha, \sigma_\alpha)$, with $\mu_\alpha = 3\%$ and $\sigma_\alpha = 1\%$, to reflect typical residual levels after telluric correction (Ulmer-Moll et al. 2019) and negative values of $\alpha$ mimic overcorrection.

We considered two cases: (i) the actual BERV values of Gl 725 B (narrow coverage), and (ii) scaling the BERV values to match those of Gl 905 (wide coverage). This allowed us to isolate the specific effect of the BERV excursion on telluric-induced RV signals.

## 3.2. Simulation results

We generated the RV time series using the LBL algorithm. To ensure realism, we did not reuse the original mask and template from Gl 905, but recomputed them from the telluric-contaminated spectrum $T(t)$.

For each BERV scenario (narrow and wide), we ran three independent simulations. As the differences between simulations

<footnote n="2">² https://lbl.exoplanets.ca/</footnote>

<footnote n="3">³ https://github.com/hippke/wotan</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Simulation results when telluric residuals were injected at a 3 ± 1% level. Top: resulting RV time series from the simulations in time (left) and BERV (right) spaces. Bottom: their respective GLS periodograms. The RVs computed with a wide and narrow BERV coverage are indicated in red and green, respectively. The significance levels of FAP = 10⁻³ and FAP = 10⁻⁵ are indicated by the solid and dashed horizontal lines, respectively. The gray line in the periodogram corresponds to the window function.</figcaption></figure>

were minor, we show results from only one simulation per case in Fig. 1, along with RVs as a function of time and BERV. We computed GLS periodograms (Zechmeister & Kürster 2009) with 10 000 frequency steps between 1.1 and 1000 days, using the *astrology* package (The Astropy Collaboration 2018). We determined the false-alarm probability (FAP) of the highest peak following Baluev (2008).

The wide and narrow BERV cases show periodic structures at telluric harmonics, with RMS values of 1.23 and 1.53 m s⁻¹, respectively. This agrees with earlier simulations by Latouf et al. (2022). As expected, even low residual telluric contamination at the percent level induces spurious RV signals. Importantly, this remains true even with wide BERV coverage, which highlights the need for advanced correction methods. Although the RMS is slightly lower in the narrow BERV case, this results from our choice to inject an identical percentage of telluric contamination in both scenarios. This might be misleading because it assumes a relatively optimistic situation for the narrow BERV case: that ~3% level residuals remain after correction. This is less realistic than for the wide BERV case, where the wider coverage makes low residuals like this more achievable. We also note that telluric contamination affects other LBL observables, such as the differential line width (dLW) that we discuss in Appendix C.

# 4. Correcting the effect of telluric absorption lines on LBL data

## 4.1. The wapiti method

As discussed in the previous section, Gl 725 B presents an interesting challenge through its limited BERV coverage, which increases the effect of residual telluric contamination on the RV measurements. To address this, we tested whether the wapiti method (Ould-Elhkim et al. 2025), a data-driven approach based

on wPCA, can correct for these systematics in the specific context of Gl 725 B. A user guide for wapiti is available online<fnref n="4" />.

In summary, the wapiti method consists of performing a wPCA on the per-line RVs and selecting the most relevant principal vectors $\vec{V}_i$ to serve as proxy time series to model and correct systematic effects in the RV data. This selection procedure is essential, rather than simply choosing an arbitrary number of leading components, because wPCA captures overall variability in the data and not specifically the systematics. As a result, the first principal component might explain the strongest variance, but not necessarily the dominant systematic signal we aim to correct for. A step-by-step description of this selection process is provided in Appendix A.

The components selected for the narrow BERV coverage were the first, second and tenth components, while for the wide BERV coverage, they were the first, third, and fourth components. One caveat of data-driven methods such as wapiti is that the selected components are not always straightforward to interpret. Ideally, we would understand the effect of telluric absorption lines without relying on wPCA-derived proxies for correction. Nonetheless, in Appendix B, we provide an interpretation of the first component, which is common to both cases.

## 4.2. Injection-recovery test

To evaluate the performance of wapiti for signal detection in RV data affected by tellurics, we carried out an injection-recovery test on the same simulated time series as constructed for Gl 725 B. We injected planetary-like signals into the per-line RVs using the form

$$K \sin\left(\frac{2\pi t}{P} + \phi\right), \qquad (2)$$

<footnote n="4">⁴ https://github.com/HkmMerwan/wapiti</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Left: detection maps showing the percentage of recovered injected signals as a function of period and semi-amplitude for the narrow and wide BERV coverages. Middle: same detection tests after applying the wapiti correction. Right: difference between the recovered and injected K values ($\Delta K = K_{out} - K_{in}$) at 1 m s⁻¹, with the shaded area indicating one standard deviation. The color bar indicates the detection rate in percent.</figcaption></figure>

where K is the semi-amplitude, P the orbital period, and φ the phase. We explored a grid of 100 periods logarithmically spaced between 2 and 800 days and 100 semi-amplitudes linearly spaced between 0.5 and 5 m s⁻¹. For each (P, K) pair, we injected the signal ten times using different phases linearly sampled between 0 and 2π.

Detection was performed using a Bayesian periodogram (Delisle et al. 2018), which compares a null model μ₀ and a planet model μₚ for each trial period,

$$RV(t) = \gamma_{SPIRou} + \sum_{i=1}^{N} a_i \hat{V}_i, \qquad (3)$$

$$RV(t) = \gamma_{SPIRou} + \sum_{i=1}^{N} a_i \hat{V}_i + A \sin(2\pi t/P) + B \sin(2\pi t/P). \quad (4)$$

We computed the Bayesian information criterion (BIC) for both models and derived the logarithm of the Bayes factor as

$$\log BF(P) = \frac{BIC_0 - BIC(P)}{2}.$$

We used a grid of 1000 points ranging from 1.5 d to 1000 d with a logarithmic spacing, and a threshold of log BF > 5 was used to claim a detection (Delisle et al. 2018).

We performed the test for both BERV scenarios (narrow and wide), with and without the wapiti correction. A signal was considered recovered when the maximum log BF exceeded 5 (Delisle et al. 2018) and was within 5% of the injected period. The resulting detection maps are shown in Fig. 2. The detection rates clearly improve after correction, particularly in the low semi-amplitude regime. The signals near one-year aliases and

at long periods remain more difficult to recover in the narrow BERV scenario, however.

To assess whether the semi-amplitude was also correctly retrieved, we tested recovery at Kₘ = 1 m s⁻¹ and compared it to the recovered value Kₐut. Fig. 2 shows that wapiti achieved an accurate amplitude recovery within the uncertainties.

In summary, the wapiti method appears to be effective in correcting systematic effects caused by telluric absorption lines. It improves the detectability of planetary signals and allows for an accurate characterization of their semi-amplitude.

# 5. Analysis of Gl 725 B

## 5.1. Activity indicators

We investigated the stellar activity of Gl 725 B using three SPIRou-derived indicators: the dLW, the longitudinal magnetic field (Bₜ), and the differential effective temperature (dET; Artigau et al. 2024). The dET indicator is a recently developed proxy for stellar activity that was derived from the LBL analysis of spectral line depth modulations in SPIRou spectra.

The corresponding time series and Bayesian periodograms for Bₜ and dET are shown in Fig. 3, and both reveal strong signals near 117 and 107 days, respectively. These periods are close to the stellar rotation period of 135 ± 15 days reported by Donati et al. (2023b). The dLW periodogram is heavily contaminated by signals of telluric origins due to the limited BERV coverage of the target (see Appendix C), which hampers the reliability of this indicator for this source.

## 5.2. Correction for stellar activity

We used the multidimensional Gaussian process (GP) framework proposed by Rajpaul et al. (2015), which was inspired by
<figure><img src="image_5.png" /><figcaption>Fig. 3. Time series (top) and Bayesian periodograms (bottom) for the longitudinal magnetic field $B_{\ell}$ (left) and the differential effective temperature activity indicator dET (right). The dashed red line shows the significance threshold of 5. The vertical dashed lines mark the periods of the two signals in RV data at 37.9 and 4.77 days.</figcaption></figure>

the work of Aigrain et al. (2012). It enables the simultaneous modeling of stellar activity and Keplerian signals in the RV time series by incorporating information from stellar activity indicators.

In this framework, the stellar activity affecting the RVs and activity indicators is described by a common latent function $GP(t)$ and its time derivative $\dot{GP}(t)$. The activity-induced RV variations are then expressed as linear combinations of $GP(t)$ and $\dot{GP}(t)$. For this analysis, we employed the s+leaf 2 GP framework<fnref n="5" />, which is designed for an efficient joint modeling of multiple time series with a computational cost that scales linearly with the dataset size (Delisle et al. 2022). We adopted a Matérn 3/2 exponential periodic (MEP) kernel, which approximates the squared-exponential periodic (SEP) kernel that was widely used to constrain stellar activity signals (Aigrain et al. 2012; Haywood et al. 2014),

$$k(\tau) = \sigma \exp \left( -\frac{\tau^2}{2\rho^2} - \frac{\sin^2\left(\frac{\pi}{P}\tau\right)}{2\eta^2} \right). \quad (5)$$

We tested two multidimensional GP models: one model including the RV and dET time series (RV + dET), and another model with the RV, dET, and $B_{\ell}$ time series (RV + dET + $B_{\ell}$). The models are described as follows:

$$RV = M_0 + \alpha_1 GP(t) + \beta_1 \dot{GP}(t),$$

$$dET = \text{offset}_{dET} + \alpha_2 GP(t) + \beta_2 \dot{GP}(t),$$

$$B_{\ell} = \text{offset}_{B_{\ell}} + \alpha_3 GP(t) + \beta_3 \dot{GP}(t),$$

$$M_0 = \gamma + \text{trend} \cdot t + \sum_{i=1}^{N} a_i \hat{V}_i + \sum_{i=1}^{N_p} \left[ A_i \sin\left(\frac{2\pi t}{P_i}\right) + B_i \cos\left(\frac{2\pi t}{P_i}\right) \right], \quad (6)$$

where $M_0$ is the base model incorporating a constant offset $\gamma$, a trend term, systematic components $\hat{V}_i$, and previously identified planetary signals at periods $P_i$ (if any).

To assess the importance of the derivative term $\dot{GP}(t)$, we explored all combinations of setting the $\beta$ coefficients to zero or allowing them to vary, using likelihood maximization to fit the models. The results (Table 2) show that in the RV + dET model, fixing both $\beta_1$ and $\beta_2$ to zero is favored. Physically, $GP(t)$ is interpreted as related to the area of the visible stellar disk covered by active regions, while its time derivative corresponds to the temporal evolution of these regions. Therefore, the lack of preference for including the derivative in the RV and dET model suggests that convective blueshift suppression is the dominant activity-induced effect in these time series (Dumusque et al. 2014). In contrast, for the RV + dET + $B_{\ell}$ model, the favored configuration also sets $\beta_1$ and $\beta_2$ to zero, but allows $\beta_3$ to vary freely. This indicates that for the $B_{\ell}$ time series, the evolution of active regions on the stellar surface plays a more significant role in the observed activity signal.

### 5.3. Search for signals in the RVs

To search for signals in the RV time series, we employed an iterative approach based on Bayesian periodograms (Delisle et al. 2018; Ould-Elhkim et al. 2025). In this method, significant signals are sequentially added to the model, and the log BF periodogram is recomputed at each step. The GP hyperparameters were initially taken from the best fit of the multidimensional model described by Eq. (6). During the period search, all GP hyperparameters were kept fixed except for the amplitude $\alpha$, while we computed the log BF over a dense grid of 10 000 trial periods between 1.5 and 1000 days by adding a circular planetary signal at each fixed period.

When a significant signal at period $P$ was detected during iteration $N_{p+1}$, a new model $M_P$ was tested by adding a sinusoidal component at that period. Because this search assumes

<footnote n="5">⁵ https://dace.unige.ch/pythonAPI/?tutorialId=703</footnote>
<table><caption>**Table 2.** Comparison of the model performances for Gl 725 B using different combinations of time series.</caption><thead><tr><th>Time series</th><th>Model</th><th>log L<sub>max</sub></th><th>χ<sup>2</sup><sub>red</sub></th><th>BIC</th><th>log BF</th><th>RMS (m/s)</th></tr></thead><tbody><tr><td rowspan="4">RV only</td><td>Offset + wapiti + trend</td><td>-552.04</td><td>1.050</td><td>1141.33</td><td>0</td><td>3.6</td></tr><tr><td>Offset + Wapiti + Trend + 1p (circular)</td><td>-540.16</td><td>1.068</td><td>1133.55</td><td>3.9</td><td>3.4</td></tr><tr><td><b>Offset + wapiti + trend + 2p (circular)</b></td><td><b>-530.63</b></td><td><b>1.090</b></td><td><b>1130.47</b></td><td><b>5.4</b></td><td><b>3.3</b></td></tr><tr><td>Offset + wapiti + trend + 2p (eccentric)</td><td>-529.67</td><td>1.113</td><td>1149.83</td><td>-4.3</td><td>3.2</td></tr><tr><td rowspan="4">RV + dET</td><td>Offset + wapiti + trend</td><td>-1023.37</td><td>1.231</td><td>2100.96</td><td>-36.6</td><td>3.8</td></tr><tr><td>Offset + wapiti + trend + GP (β<sub>1</sub> = β<sub>2</sub> = 0)</td><td>-971.71</td><td>0.768</td><td>2027.76</td><td>0</td><td>3.5</td></tr><tr><td>Offset + wapiti + trend + GP (β<sub>1</sub> = β<sub>2</sub> = 0) + 1p (circular)</td><td>-958.56</td><td>0.770</td><td>2019.52</td><td>4.1</td><td>3.3</td></tr><tr><td><b>Offset + wapiti + trend + GP (β<sub>1</sub> = β<sub>2</sub> = 0) + 2p (circular)</b></td><td><b>-947.44</b></td><td><b>0.751</b></td><td><b>2015.35</b></td><td><b>6.2</b></td><td><b>3.1</b></td></tr><tr><td rowspan="4">RV + dET + B<sub>ℓ</sub></td><td>Offset + wapiti + trend</td><td>-1680.02</td><td>1.754</td><td>3430.79</td><td>-34.8</td><td>3.8</td></tr><tr><td>Offset + wapiti + trend + GP (β<sub>1</sub> = β<sub>2</sub> = 0)</td><td>-1622.75</td><td>1.229</td><td>3361.27</td><td>0</td><td>3.5</td></tr><tr><td>Offset + wapiti + trend + GP (β<sub>1</sub> = β<sub>2</sub> = 0) + 1p (circular)</td><td>-1609.92</td><td>1.286</td><td>3354.89</td><td>3.2</td><td>3.3</td></tr><tr><td><b>Offset + wapiti + trend + GP (β<sub>1</sub> = β<sub>2</sub> = 0) + 2p (circular)</b></td><td><b>-1598.76</b></td><td><b>1.324</b></td><td><b>3351.86</b></td><td><b>4.7</b></td><td><b>3.1</b></td></tr><tr><td></td><td>Offset + wapiti + trend + GP (β<sub>1</sub> = β<sub>2</sub> = 0) + 2p (eccentric)</td><td>-1597.59</td><td>1.356</td><td>3381.69</td><td>-10.2</td><td>3.1</td></tr></tbody></table>

**Notes.** Maximum log-likelihood, reduced chi-squared, BIC, log BF, and RMS of the residuals for the best model applied to each combination of time series. The model in bold indicates the best-fit model to the combination of time series.

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Iterative Bayesian periodograms for the Gl 725 B RV time series under different models. From left to right: (1) Raw RVs after trend correction, (2) after removal of systematics via wapiti, (3) with the 38-day planet, and (4) residuals after modeling the 38-day and 4.77-day signals. The dashed red line marks the detection threshold at log BF = 5.</figcaption></figure>

fixed periods and GP hyperparameters, it can slightly overestimate log BF. After each detection, we therefore refit the full model to allow the orbital period and all GP hyperparameters to vary freely in order to verify the statistical robustness of the signal.

The wapiti method applied on Gl 725 B dataset selected the second, ninth, first, and seventh components (see Appendix A for a detailed explanation of this component selection), which successfully removed systematic signals at periods of 180 days, 1 year, and 120 days. The iterative Bayesian periodograms are shown in Fig. 4. After applying the wapiti correction, we detected a 38-day signal, and when this was added to the model, another 4.77 d signal was detected in the residuals. Finally, after adding this second signal, no further signals were detected. We found no trace of the 2.7 d signal reported by Berdiñas et al. (2016), nor evidence of the 5.5-year companion of Baize (1976). Based on their parameters, if such a companion existed, it would produce a semi-amplitude of about 600 m s⁻¹. This amplitude is clearly absent from our SPIRou RVs (Fig. 5).

Consequently, we fit a one-planet (1p) model including the 38-day signal and a two-planet (2p) model including the 38-day and 4.77-day signals, considering both circular and Keplerian orbits, by maximizing their likelihood. The results are presented in Table 2. The RV + dET model yielded the best result, with log BF of 4.1 and 6.2 for the one-planet and two-planet cases,

respectively, compared to the model without planetary signals (see Table 2).

These log BF values were computed from the maximum likelihood values alone and might misrepresent the true evidence (Perrakis et al. 2014), however. A more rigorous assessment of the signal significance, for example, using the false-inclusion probability (FIP; Hara et al. 2022) periodogram, which makes use of the full posterior distribution, is presented in Appendix D.2. In this posterior-based analysis, only the 38-day signal is robustly recovered, while the 4.77-day signal does not reach the detection threshold.

Additionally, we further verified the robustness of these periodicities with respect to the choice of the NZP correction. This sensitivity analysis, presented in Appendix D.3, confirmed that the 38-day signal was consistently recovered in the NZP realizations, whereas the short-period signal is more sensitive to the specific NZP choice and was therefore considered as a candidate rather than a firm detection.

## 5.4. Characterization of the system

The absence of the two signals in the stellar activity indicators, combined with their robustness and statistical significance even when stellar activity was modeled in the RV data, supports a planetary interpretation for the 38-day signal and suggests that
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Two-planet GP model for Gl 725 B. **Top:** raw RVs with trend. **Middle:** phase-folded RVs for the 38-day and 4.77-day signals with binned points (black). **Bottom:** modeled RV and dET time series (with GP) and corresponding residuals.</figcaption></figure>

the 4.77-day signal is a plausible planetary candidate. To characterize the detected signals, we used the adaptive MCMC sampler samsam<fnref n="6" /> (Delisle et al. 2018). The sampler was configured to perform a total of 1000 000 iterations, and we used the final 750 000 iterations for the statistical analysis. We ran two cases, one case with circular orbits, and the one with Keplerian orbits.

The results are presented in Table 3 and show that the eccentricities of the confirmed planet and the candidate signal cannot be well constrained with the current dataset. This is expected given the low amplitude of the signals, and it is consistent with

Table 2, where the Keplerian model was not favored over the circular model. Nevertheless, all other orbital parameters remained consistent between the circular and eccentric fits. The circular model therefore likely provides a more reliable representation of the system. The corner plot for this final circular fit is shown in Appendix D.1.

Our multidimensional GP analysis yields a stellar rotation period of $P_{rot} = 105.1 \pm 3.3$ d, which is more precise than the previously reported value of 135 ± 15 d based on the $B_{\ell}$ time series (Donati et al. 2023b). The difference in value between these two estimates (but they are still compatible within $2\sigma$) might arise from our earlier finding that for the $B_{\ell}$ time series,

<footnote n="6">⁶ https://gitlab.unige.ch/Jean-Baptiste.Delisle/samsam</footnote>
<table><caption>**Table 3.** Orbital parameters of Gl 725 B, showing the best-fit model including one confirmed planet (c) and one candidate (b) with circular and Keplerian fits.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="2">2p</th><th rowspan="2">Priors</th></tr><tr><th>Circular</th><th>Keplerian</th></tr></thead><tbody><tr><td colspan="4"><strong>Instrumental parameters</strong></td></tr><tr><td>γ<sub>S PIRou</sub> (m s<sup>-1</sup>)</td><td>1143.7 ± 5.7</td><td>1143.9 ± 5.3</td><td>U(-∞, ∞)</td></tr><tr><td>trend (m s<sup>-1</sup> d<sup>-1</sup>)</td><td>-0.016 ± 0.001</td><td>-0.016 ± 0.001</td><td>U(-∞, ∞)</td></tr><tr><td>σ<sub>S PIRou</sub> (m s<sup>-1</sup>)</td><td>2.9 ± 0.2</td><td>2.9 ± 0.2</td><td>LU(10<sup>-4</sup>, 10<sup>4</sup>)</td></tr><tr><td colspan="4"><strong>wapiti systematics</strong></td></tr><tr><td>a<sub>1</sub> (m/s)</td><td>-34.1 ± 3.4</td><td>-34.2 ± 3.4</td><td>U(-∞, ∞)</td></tr><tr><td>a<sub>2</sub> (m/s)</td><td>-25.5 ± 3.3</td><td>-25.7 ± 3.3</td><td>U(-∞, ∞)</td></tr><tr><td>a<sub>3</sub> (m/s)</td><td>-19.7 ± 3.4</td><td>-19.6 ± 3.4</td><td>U(-∞, ∞)</td></tr><tr><td>a<sub>4</sub> (m/s)</td><td>15.8 ± 3.5</td><td>16.1 ± 3.5</td><td>U(-∞, ∞)</td></tr><tr><td colspan="4"><strong>dDET parameters</strong></td></tr><tr><td>dDET (K)</td><td>0.5 ± 8.9</td><td>0.2 ± 7.8</td><td>U(-∞, ∞)</td></tr><tr><td>σ<sub>dDET</sub> (K)</td><td>1.6 ± 0.1</td><td>1.6 ± 0.1</td><td>LU(10<sup>-4</sup>, 10<sup>4</sup>)</td></tr><tr><td colspan="4"><strong>Stellar activity hyperparameters</strong></td></tr><tr><td>P<sub>rot</sub> (days)</td><td>105.1 ± 3.3</td><td>105.1 ± 3.4</td><td>LU(0.1, 10<sup>5</sup>)</td></tr><tr><td>ρ (days)</td><td>2160 ± 1697</td><td>2072 ± 1644</td><td>LU(0.1, 10<sup>5</sup>)</td></tr><tr><td>η</td><td>1.3 ± 0.5</td><td>1.3 ± 0.5</td><td>LU(0.05, 20)</td></tr><tr><td>α<sub>1</sub> (m/s)</td><td>-5.3 ± 2.8</td><td>-5.0 ± 2.6</td><td>N(0, 10)</td></tr><tr><td>α<sub>2</sub> (K)</td><td>10.1 ± 4.9</td><td>9.6 ± 4.7</td><td>N(0, 10)</td></tr><tr><td colspan="4"><strong>Orbital parameters</strong></td></tr><tr><td>P<sub>b</sub> (days)</td><td>4.765 ± 0.004</td><td>4.765 ± 0.004</td><td>LU(0.1, 10<sup>5</sup>)</td></tr><tr><td>K<sub>b</sub> (m s<sup>-1</sup>)</td><td>1.4 ± 0.3</td><td>1.4 ± 0.4</td><td>LU(0.1, 10<sup>4</sup>)</td></tr><tr><td>λ<sub>b</sub>(0) (deg)</td><td>-0.3 ± 3.8</td><td>-0.3 ± 3.9</td><td>U(0, 2π)</td></tr><tr><td>e<sub>b</sub> cos ω<sub>b</sub></td><td>—</td><td>0.00 ± 0.12</td><td>modBeta(0.867, 3.03)</td></tr><tr><td>e<sub>b</sub> sin ω<sub>b</sub></td><td>—</td><td>-0.03 ± 0.13</td><td>modBeta(0.867, 3.03)</td></tr><tr><td>P<sub>c</sub> (days)</td><td>37.90 ± 0.17</td><td>37.91 ± 0.16</td><td>LU(0.1, 10<sup>5</sup>)</td></tr><tr><td>K<sub>c</sub> (m s<sup>-1</sup>)</td><td>1.7 ± 0.3</td><td>1.7 ± 0.4</td><td>LU(0.1, 10<sup>4</sup>)</td></tr><tr><td>λ<sub>c</sub>(0) (deg)</td><td>5.2 ± 2.6</td><td>5.4 ± 2.5</td><td>U(0, 2π)</td></tr><tr><td>e<sub>c</sub> cos ω<sub>c</sub></td><td>—</td><td>0.05 ± 0.14</td><td>modBeta(0.867, 3.03)</td></tr><tr><td>e<sub>c</sub> sin ω<sub>c</sub></td><td>—</td><td>-0.03 ± 0.14</td><td>modBeta(0.867, 3.03)</td></tr><tr><td colspan="4"><strong>Derived parameters</strong></td></tr><tr><td>M<sub>p,b</sub> sin i (M<sub>⊕</sub>)</td><td>1.5 ± 0.4</td><td>1.4 ± 0.4</td><td>—</td></tr><tr><td>a<sub>b</sub> (AU)</td><td>0.035 ± 0.001</td><td>0.035 ± 0.001</td><td>—</td></tr><tr><td>e<sub>b</sub></td><td>—</td><td>0.13 ± 0.11</td><td>Beta(0.867, 3.03)</td></tr><tr><td>ω<sub>b</sub></td><td>—</td><td>-0.4 ± 1.8</td><td>—</td></tr><tr><td>T<sub>eq,b</sub> (K)</td><td>430 ± 8</td><td>430 ± 8</td><td>—</td></tr><tr><td>M<sub>p,c</sub> sin i (M<sub>⊕</sub>)</td><td>3.4 ± 0.7</td><td>3.5 ± 0.7</td><td>—</td></tr><tr><td>a<sub>c</sub> (AU)</td><td>0.139 ± 0.004</td><td>0.139 ± 0.004</td><td>—</td></tr><tr><td>e<sub>c</sub></td><td>—</td><td>0.14 ± 0.12</td><td>Beta(0.867, 3.03)</td></tr><tr><td>ω<sub>c</sub></td><td>—</td><td>-0.3 ± 1.5</td><td>—</td></tr><tr><td>T<sub>eq,c</sub> (K)</td><td>215 ± 4</td><td>215 ± 4</td><td>—</td></tr></tbody></table>

the evolution of active regions contributes more significantly to the observed activity signal, whereas this effect appears to be less prominent in the RV and dDET time series.

As for the planetary signals, the minimum mass of these planets increases inside-out, ranging from 1.5 ± 0.4 $M_⊕$ for the candidate Gl 725 Bb to 3.4 ± 0.7 $M_⊕$ for Gl 725 Bc. According to the classification by Mishra et al. (2023), the architecture of this system falls into the “ordered” class. This class was identified as relatively rare in simulations, with a 1.5% occurrence rate. In contrast, their observation sample indicates a higher prevalence, with 37% of observed systems falling into this class. Therefore, our findings support the idea that ordered planetary systems may be more common than predicted by simulations. Based on Kopparapu et al. (2014) and using the stellar parameters from Table 3, the habitable zone (HZ) of Gl 725 B for an Earth-mass planet spans from 0.09–0.19 AU. Consequently, Gl 725 Bc falls within the habitable zone of its host star. We discuss the localization of Gl 725 Bc in the habitable zone in more detail in Sect. 6.3.

We did not have to take the effect of the binary companion Gl 725 A into account because both planets lie at a distance

that is smaller than the critical semimajor axis of the system, $a_c = 12.4$ AU. This value was calculated using the first formula from Holman & Wiegert (1999) and the orbital parameters from Cortés-Zuleta et al. (2025) for the binary system. Using these parameters, we also calculated an expected trend of 0.0238 ± 0.0059 m s⁻¹·d⁻¹, which is consistent within uncertainties with the trend obtained from the MCMC analysis of 0.016 ± 0.001 m s⁻¹·d⁻¹.

## 5.5. Analysis of the TESS photometric data

We used the transit least squares algorithm<fnref n="7" /> (TLS; Hippke & Heller 2019) to search for transit events in the blended light curves of Gl 725 A and B. No transits were identified by the TLS, including periodic events at the orbital period of planet b. Because of the long time span of the TESS data, a periodic event of 4 days or even longer than 30 days is expected to be easy to identify. The derived planetary parameters suggest a transit probability of only 3.7% for planet b and 0.9% for planet c, however, which might explain the no-transiting nature of the two planets.

To assess the transit detection sensitivity, we conducted an injection-recovery test in our TESS data for a range of orbital inclinations and planet radii following Cortés-Zuleta et al. (2025). For each planet, we performed 50 000 simulations of transit events generated by the batman Python package (Kreidberg 2015). The planet parameters were taken from Table 3, and the limb-darkening coefficients were taken from Table 15 in Claret (2017). For each run, the orbital inclination and planet radius were drawn randomly from uniform distributions between 85° and 90° for the inclination and 0.5 to $10R_⊕$ for the radius. Then, we applied the box least squares (BLS, Kovács et al. 2002) periodogram from Astropy<fnref n="8" /> to measure the power of the signal at the expected orbital period of planet b (4.765 days) and planet c (37.90 days).

We binned the orbital inclination and planet radius in bins of 0.1 units for each parameter and then computed the fraction of transit events that are detectable with the BLS periodogram (i.e., power >2000) in each bin. Since the light curves are blended, the true detectable planet radius needs to be scaled by the dilution factor, which has a value of 1.64. The results are shown in Fig. 6, where we found that in general, the detection limit of the TESS photometry is about $1.0R_⊕$ for aligned planets. For planet b, most of the planets with a radius greater than $1.0R_⊕$ and orbital inclinations between 87.5° and 90° could be detected. In the case of planet c, the transit would have been detected for orbital inclinations between 89.4° and 90° and for a planet radius larger than 1.2 $R_⊕$. For both cases, smaller planets are detectable for inclinations closer to 90°, as expected.

Using the spright tool (Parviainen et al. 2024), we predict the radius of Gl 725 Bb to be $1.2^{+0.2}_{-0.1} R_⊕$ and that of Gl 725 Bc to be $1.7^{+0.4}_{-0.3} R_⊕$ (see Fig. 6). The posterior distribution for planet c is bimodal, which is known as the “radius valley” (Fulton et al. 2017; Cloutier & Menou 2020; Gaidos et al. 2024), and it was approximated here as the sum of two Gaussian components centered at $1.41 ± 0.09 R_⊕$ and $1.82 ± 0.20 R_⊕$. The transit events of planet b are undetectable if the planet radius is $≤1.0 R_⊕$, which is unlikely from the posterior distribution of predicted planet radii of Fig. 6. Similarly, planet c is unlikely to be transiting because we predict a planet radius larger than $1.2 R_⊕$, which is detectable in the TESS photometry.

<footnote n="7">⁷ https://github.com/hippke/tls</footnote>

<footnote n="8">⁸ https://docs.astropy.org/en/stable/timeseries/bls.html</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Detectability of transit events of the candidate planet b (top) and the confirmed planet c (middle) in the TESS data of Gl 725 AB as a function of orbital inclination and planet radius. Bottom: posterior distributions of the planetary radii for Gl 725 Bb (blue) and Gl 725 Bc (brown) inferred using the spright tool (Parviainen et al. 2024).</figcaption></figure>

# 6. Discussion and conclusions

## 6.1. Addressing the RV effects of telluric absorption lines in the NIR using wapiti

To understand and mitigate the effect of telluric absorption lines on RV measurements for Gl 725 B, we performed controlled simulations by injecting realistic telluric residuals at the 3 ± 1% level into a clean stellar template. These simulations were designed to replicate the BERV constraints of Gl 725 B and showed that even modest telluric residuals induce systematic signals in the RVs extracted with the LBL algorithm. This reinforces the need for dedicated correction methods, such as wapiti, especially for targets with limited BERV coverage.

Using the per-line RV time series from our simulations, we conducted injection-recovery tests and found that wapiti successfully corrected for the effects of telluric contamination. In the narrow and wide BERV scenarios, it improved the detectability and enabled an accurate recovery of the planetary signal amplitudes.

While this is promising, it is important to recognize the limitations of these simulations. They represent a simplified

environment in which the star is modeled as a static template plus photon noise. In reality, additional complexities arise, including seasonal or chromatic variability in telluric residuals (α), night-sky emission lines (OH and Na), Moon contamination (Crepp et al. 2016), and variations in the blaze function due to instrumental or illumination changes (Halverson et al. 2015). These are not captured in our controlled setup, but would affect real observations.

In addition, a core limitation of wapiti (or any wPCA-based methods) is that wPCA is designed to capture variability rather than systematics. Although the wapiti method is implemented in a way that accounts for this limitation, it remains a key aspect to consider because many of its capabilities and limitations arise from it. For example, because of this tendency to capture variability, wPCA is highly sensitive to outliers and observations with a low signal-to-noise ratio, and it might be necessary to exclude some measurements in certain situations.

Finally, since wPCA constructs orthogonal components by design, real-world systematics that are not orthogonal might become entangled within the same vector. This would complicate the interpretation. Despite these limitations, wapiti remains a powerful tool for mitigating telluric-induced variability, as we demonstrate below in the analysis of real data from Gl 725 B.

## 6.2. Stellar age

Our multidimensional GP analysis retrieved a rotational period of 105.1 ± 3.3 d for Gl 725 B. We combined this value with the T<sub>eff</sub> derived by Cristofari et al. (2022) with the gyrochronology presented by Gaidos et al. (2023) to estimate an age of 8.7 ± 2.1 Gyr. This is consistent with the old but very uncertain age of Gl 725A (≥7 Gyr, Fouqué et al. 2023; Cortés-Zuleta et al. 2025). The correction for a nonsolar metallicity (different from the stars used to calibrate the gyrochronology) was estimated at +0.8 Gyr, but this is only based on theory. Three other caveats are that (1) the estimated age of Gl 725B is much older than the oldest calibrator (the 4 Gyr-old M67 cluster Dungee et al. 2022), and Skumanich-like spin-down (Skumanich 1972) is assumed after this point, (2) Gl 725B is near the cool end of the valid range, and (3) Gl 725B is in a binary system with a projected separation of 65 au, where its protostellar disk would have been truncated, likely leading to more rapid dissipation, less angular momentum exchange with the star, and hence, more rapid initial stellar rotation, as was observed in young clusters (e.g., Kraus et al. 2012; Messina 2019). Nevertheless, that the star and its planets are old is likely a reliable result. The kinematics of the system (U = -25.33, V = -12.95, W = 26.58 km sec<sup>-1</sup>) are consistent with the Galactic thin disk and are not particularly suggestive of nor inconsistent with an old age.

## 6.3. Planetary system Gl 725 B

Through the wapiti correction, we uncovered evidence of a compact ordered planetary system orbiting Gl 725 B. The outer planet, Gl 725 Bc, is firmly detected with a minimum mass of 3.4 ± 0.7 M<sub>⊕</sub> and an orbital period of 37.90 ± 0.17 days. In addition, we identified a shorter-period signal at 4.765 ± 0.004 days that might correspond to a possible inner planet candidate, Gl 725 Bb, with a minimum mass of 1.5 ± 0.4 M<sub>⊕</sub>. Further observations are required to confirm its planetary nature. The system thus appears to consist of at least one confirmed planet (Gl 725 Bc) and a promising close-in candidate.
Recent estimates of the luminosity of Gl725B derived from stellar radius measurements, empirical relations, or the fitting of synthetic spectra consistently converged to between $8.4 \times 10^{-3}$ and $9.2 \times 10^{-3}$ L⊙ (Cristofari et al. 2022; Boyajian et al. 2012; Mann et al. 2015). Based on the semi-major axes of 0.035 AU and 0.139 AU for planets Gl725Bb and Gl725Bc, respectively, their received stellar fluxes are approximately 7 and 0.45 times that of Earth. The latter value places Gl725Bc in an insolation regime analogous to that of Mars, and it positions it securely within the nominal circumstellar habitable zone.

Significant recent progress has been made in characterizing the inner boundary of the habitable zone, particularly in relation to the onset of the runaway greenhouse effect. This has been achieved using both 1D climate models (Kasting et al. 1993; Kopparapu et al. 2014; Yang et al. 2016; Chaverot et al. 2022) and 3D global climate models (Wordsworth et al. 2011; Yang et al. 2014; Kopparapu et al. 2017; Turbet et al. 2018; Wolf et al. 2020). Furthermore, recent investigations have explored the possibility of climate bistability or multistability near this inner edge (Turbet et al. 2023), highlighting the necessity of considering not only the current insolation of a planet, but also its temporal evolution when the long-term habitability is assessed. In contrast, fewer studies have addressed the outer edge of the habitable zone, likely because only a limited number of known planets are situated in this regime. Nonetheless, Ramirez & Kaltenegger (2017, 2018) investigated the evolution of the outer boundary in the habitable zone as a function of atmospheric composition using 1D climate models.

For M-dwarf stars with effective temperatures of about 3400 K, as is the case for Gl725 B, the conservative inner edge of the habitable zone is generally situated at a stellar flux of approximately 0.8 $S_⊕$, regardless of historical insolation trends (Turbet et al. 2023). The outer boundary, assuming a classical $N_2-CO_2-H_2O$ atmospheric composition, is situated around 0.25 $S_⊕$. This limit can vary slightly depending on the atmospheric makeup: A hydrogen-rich atmosphere (e.g., from volcanic outgassing) can push the boundary outward, while methane-rich atmospheres tend to shift it inward (Ramirez & Kaltenegger 2017). The outer edge, assuming a canonical $N_2-CO_2-H_2O$ atmospheric composition, lies near 0.25 $S_⊕$. This limit might be modulated by alternative atmospheric constituents: Enhanced $H_2$ abundances (e.g., via volcanic outgassing) can extend the outer edge outward, whereas elevated $CH_4$ levels might induce a modest inward shift (Ramirez & Kaltenegger 2017). In all plausible scenarios, Gl725 Bc remains well ensconced within the habitable zone, however.

The source Gl725 Bc is therefore the second closest known Earth or super-Earth-type planet located within the habitable zone of its star after Proxima b. Although its true mass is unknown, we can estimate its most probable value by applying the average correction factor for inclination ($\sqrt{4/3} \sim 1.15$), yielding an estimated mass of 4 $M_⊕$. According to Otegi et al. (2020), this mass range strongly suggests a rocky, terrestrial composition.

Atmospheric characterization of exoplanets is commonly performed through transmission spectroscopy, which requires transiting planets. This is an unlikely configuration for Gl 725 Bb and Gl 725 Bc. Alternative techniques have been suggested to study the atmospheres of nontransiting nearby exoplanets, however. One promising approach combines high spectral and spatial resolution and is known as high-dispersion coronagraphy (HDC; Lovis et al. 2017; Snellen et al. 2019).

Using the formula from Lovis et al. (2017) and assuming geometric albedos in the range $A_g = 0.1–0.6$ and planetary radii of

1.1–1.4 $R_⊕$ for Gl 725 Bb and 1.4–2.1 $R_⊕$ for Gl 725 Bc, we estimated that the planet-to-star flux ratio lies between $6 \times 10^{-8}$ and $3.3 \times 10^{-6}$ for Gl 725 Bb and between $5 \times 10^{-8}$ and $3.8 \times 10^{-6}$ for Gl 725 Bc. Their maximum angular separations, derived from Table 3, are 9.1 and 39.7 mas. Based on the diffraction limit of the ELT at 1300 nm ($2\lambda/D \approx 13.7$ mas; Suárez Mascareño et al. 2023), Gl 725 Bc would in principle fall within reach of HDC with a next-generation high-resolution spectrograph such as ANDES (Marconi et al. 2022), which will operate in the 400–1800 nm range at R ~ 100 000. Compared to Proxima b, the Earth-like planet in the habitable zone for which reflected-light characterization is currently considered the most accessible, the maximum angular separation of Gl725 Bc is slightly more favorable. The star–planet contrast might be up to twice as high, however, depending on the actual albedos of the two planets. This makes Gl725 Bc a particularly valuable target for advancing our understanding of the diversity of planets within the habitable zone.

Unfortunately, Gl 725 B lies in the northern hemisphere and is not accessible to southern facilities such as the ELT. The study of Gl 725 Bc would therefore require a comparable HDC instrument located in the north. The peculiar configuration of the BERV excursion for this star means that atmospheric characterization would be optimized using an instrument with a resolving power greater than 70 000 to enable a proper telluric correction. Alternatively, this might be conducted from space.

Space-based missions, on the other hand, offer an attractive alternative because they are not subject to telluric absorption. The proposed LIFE mission (Lichtenberg & Quanz 2023), a mid-infrared nulling interferometer, aims to detect and characterize the atmospheres of Earth-sized and super-Earth exoplanets through their thermal emission. LIFE would consist of four 2-meter apertures operating in the 4–18.5, µm range and would be particularly well suited to characterizing a planet such as Gl 725 Bc.

## Data availability

The data used in this work were recorded in the context of the SLS and are available to the public at the Canadian Astronomy Data Center. The analysis of these data is available in a *wapiti* tutorial at https://github.com/HkmMerwan/wapiti.

*Acknowledgements.* Based on observations obtained at the Canada-France-Hawaii Telescope (CFHT) which is operated from the summit of Maunakea by the National Research Council of Canada, the Institut National des Sciences de l'Université of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. The observations at the Canada-France-Hawaii Telescope were performed with care and respect from the summit of Maunakea which is a significant cultural and historic site. Based on observations obtained with SPIRou, an international project led by Institut de Recherche en Astrophysique et Planétologie, Toulouse, France. This study has been partially supported through the grant EUR TESS No. ANR-18-EURE-0018 in the framework of the Programme des Investissements d'AveNIR. This research made use of the following software tools: NumPy; a fundamental package for scientific computing with Python. *pandas*; a library providing easy-to-use data structures and data analysis tools. SciPy; a library for scientific computing and technical computing. *matplotlib*; a plotting library for the Python programming language. *tqdm*; a library for creating progress bars in the command line. *seaborn*; a data visualization library based on *matplotlib*. *wPCA*; a python package for weighted principal component analysis. *samsam*; an adaptive MCMC sampler.

## References

* Aigrain, S., Pont, F., & Zucker, S. 2012, MNRAS, 419, 3147

* Alonso-Floriano, F. J., Morales, J. C., Caballero, J. A., et al. 2015, A&A, 577, A128
Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437

Artigau, É., Astudillo-Defru, N., Delfosse, X., et al. 2014, SPIE Conf. Ser., 9149, 914905

Artigau, É., Cadieux, C., Cook, N. J., et al. 2022, AJ, 164, 84

Artigau, É., Cadieux, C., Cook, N. J., et al. 2024, AJ, 168, 252

Baize, P. 1976, A&AS, 26, 177

Baluev, R. V. 2008, MNRAS, 385, 1279

Bauer, F. F., Zechmeister, M., Kaminski, A., et al. 2020, A&A, 640, A50

Bean, J. L., Seifahrt, A., Hartman, H., et al. 2010, ApJ, 713, 410

Beriñas, Z. M., Amado, P. J., Anglada-Escudé, G., Rodríguez-López, C., & Barnes, J. 2016, MNRAS, 459, 3551

Berdyugina, S. V. 2005, Liv. Rev. Sol. Phys., 2, 8

Bertaux, J.-L., Lallement, R., Ferron, S., & Boonne, C. 2014, in 13th International HITRAN Conference, 8

Bonfils, X., Delfosse, X., Udry, S., et al. 2013, A&A, 549, A109

Boyajian, T. S., von Braun, K., van Belle, G., et al. 2012, ApJ, 757, 112

Buchner, J. 2021, J. Open Source Softw., 6, 3001

Cale, B. L., Reefe, M., Plavchan, P., et al. 2021, AJ, 162, 295

Carleo, I., Sanna, N., Gratton, R., et al. 2016, Exp. Astron., 41, 351

Carmona, A., Delfosse, X., Bellotti, S., et al. 2023, Near-IR and optical radial velocities of the active M dwarf star Gl 388 (AD Leo) with SPIRou at CFHT and SOPHIE at OHP: A 2.23 day rotation period and no evidence for a corotating planet

Cegla, H. M., Watson, C. A., Shelyag, S., Mathioudakis, M., & Moutari, S. 2019, ApJ, 879, 55

Charpentier, P., Moutou, C., Donati, J.-F., et al. 2025, A&A, 701, A17

Chaverot, G., Turbet, M., Bolmont, E., & Leconte, J. 2022, A&A, 658, A40

Claret, A. 2017, A&A, 600, A30

Cloutier, R., & Menou, K. 2020, AJ, 159, 211

Collier Cameron, A., Ford, E. B., Shahaf, S., et al. 2021, MNRAS, 505, 1699

Cook, N. J., Artigau, É., Doyon, R., et al. 2022, APERO: A Pipeline to Reduce Observations, Astrophysics Source Code Library [record ascl:2211.019]

Cortés-Zuleta, P., Boisse, I., Klein, B., et al. 2023, A&A, 673, A14

Cortés-Zuleta, P., Boisse, I., Ould-Elhkim, M., et al. 2025, A&A, 693, A164

Crepp, J. R., Crass, J., King, D., et al. 2016, SPIE Conf. Ser., 9908, 990819

Cretignier, M., Dumusque, X., Aigrain, S., & Pepe, F. 2023, A&A, 678, A2

Cristofari, P. I., Donati, J. F., Masseron, T., et al. 2022, MNRAS, 516, 3802

Crockett, C. J., Mahmud, N. I., Prato, L., et al. 2012, ApJ, 761, 164

Delchambre, L. 2015, MNRAS, 446, 3545

Delisle, J. B., Ségransan, D., Dumusque, X., et al. 2018, A&A, 614, A133

Delisle, J. B., Unger, N., Hara, N. C., & Ségransan, D. 2022, A&A, 659, A182

Donati, J. F., Kouach, D., Moutou, C., et al. 2020, MNRAS, 498, 5684

Donati, J.-F., Cristofari, P. I., Finociety, B., et al. 2023a, MNRAS, 525, 455

Donati, J.-F., Lehmann, L. T., Cristofari, P. I., et al. 2023b, MNRAS, 525, 2015

Dorn, R. J., Bristow, P., Smoker, J. V., et al. 2023, A&A, 671, A24

Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45

Dumusque, X. 2018, A&A, 620, A47

Dumusque, X., Udry, S., Lovis, C., Santos, N. C., & Monteiro, M. J. P. F. G. 2011, A&A, 525, A140

Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132

Dungee, R., van Saders, J., Gaidos, E., et al. 2022, ApJ, 938, 118

Faria, J. P., Suárez Mascareño, A., Figueira, P., et al. 2022, A&A, 658, A115

Fischer, D. A., Anglada-Escude, G., Arriagada, P., et al. 2016, PASP, 128, 066001

Fouqué, P., Martioli, E., Donati, J. F., et al. 2023, A&A, 672, A52

Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109

Gaia Collaboration 2020, VizieR Online Data Catalog: I/350

Gaidos, E., Mann, A. W., Kraus, A. L., & Ireland, M. 2016, MNRAS, 457, 2877

Gaidos, E., Claytor, Z., Dungee, R., Ali, A., & Feiden, G. A. 2023, MNRAS, 520, 5283

Gaidos, E., Ali, A., Kraus, A. L., & Rowe, J. F. 2024, MNRAS, 534, 3277

Halverson, S., Roy, A., Mahadevan, S., & Schwab, C. 2015, ApJ, 814, L22

Hara, N. C., Unger, N., Delisle, J.-B., Díaz, R. F., & Ségransan, D. 2022, A&A, 663, A14

Haywood, R. D., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517

Henry, T. J., Jao, W.-C., Subasavage, J. P., et al. 2006, AJ, 132, 2360

Hippke, M., & Heller, R. 2019, A&A, 623, A39

Hippke, M., David, T. J., Mulders, G. D., & Heller, R. 2019, AJ, 158, 143

Hirano, T., Kuzuhara, M., Kotani, T., et al. 2020, PASJ, 72, 93

Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621

Hsu, D. C., Ford, E. B., & Terrien, R. 2020, MNRAS, 498, 2249

Huélamo, N., Figueira, P., Bonfils, X., et al. 2008, A&A, 489, L9

Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108

Kervella, P., Arenou, F., Mignard, F., & Thévenin, F. 2019, A&A, 623, A72

Kirkpatrick, J. D., Henry, T. J., & McCarthy, Donald W., J. 1991, ApJS, 77, 417

Kopparapu, R. K., Ramirez, R. M., SchottelKotte, J., et al. 2014, ApJ, 787, L29

Kopparapu, R. K., Wolf, E. T., Arney, G., et al. 2017, ApJ, 845, 5

Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369

Kraus, A. L., Ireland, M. J., Hillenbrand, L. A., & Martinache, F. 2012, ApJ, 745, 19

Kreidberg, L. 2015, PASP, 127, 1161

Larue, P., Delfosse, X., Carmona, A., et al. 2025, A&A, 701, A216

Latouf, N., Wang, S. X., Cale, B., & Plavchan, P. 2022, AJ, 164, 212

Lichtenberg, T., & Quanz, S. P. 2023, in EGU General Assembly Conference Abstracts, EGU General Assembly Conference Abstracts, EGU-13634

Lovis, C., Snellen, I., Mouillet, D., et al. 2017, A&A, 599, A16

Mahmud, N. I., Crockett, C. J., Johns-Krull, C. M., et al. 2011, ApJ, 736, 123

Mann, A. W., Feiden, G. A., Gaidos, E., Boyajian, T., & von Braun, K. 2015, ApJ, 804, 64

Marchwinski, R. C., Mahadevan, S., Robertson, P., Ramsey, L., & Harder, J. 2015, ApJ, 798, 63

Marconi, A., Abreu, M., Adibekyan, V., et al. 2022, SPIE Conf. Ser., 12184, 1218424

Messina, S. 2019, A&A, 627, A97

Metcalf, A. J., Anderson, T. B., Bender, C. F., et al. 2019, Optica

Mignon, L., Delfosse, X., Meunier, N., et al. 2025, A&A, 700, A146

Mishra, L., Alibert, Y., Udry, S., & Mordasini, C. 2023, A&A, 670, A68

Moutou, C., Delfosse, X., Petit, A. C., et al. 2023, A&A, 678, A207

Otegi, J. F., Bouchy, F., & Helled, R. 2020, A&A, 634, A43

Ould-Elhkim, M., Moutou, C., Donati, J. F., et al. 2023, A&A, 675, A187

Ould-Elhkim, M., Moutou, C., Donati, J.-F., et al. 2025, A&A, 696, A152

Parviainen, H., Luque, R., & Palle, E. 2024, MNRAS, 527, 5693

Perrakis, K., Ntzoufras, I., & Tsionas, E. G. 2014, Computat. Statist. Data Anal., 77, 54

Plavchan, P., Latham, D., Gaudi, S., et al. 2015, arXiv e-prints [arXiv:1503.01770]

Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S., & Roberts, S. 2015, MNRAS, 452, 2269

Ramirez, R. M., & Kaltenegger, L. 2017, ApJ, 837, L4

Ramirez, R. M., & Kaltenegger, L. 2018, ApJ, 858, 72

Rayner, J., Tokunaga, A., Jaffe, D., et al. 2022, PASP, 134, 015002

Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, SPIE Conf. Ser., 9143, 914320

Robertson, P., Stefansson, G., Mahadevan, S., et al. 2020, ApJ, 897, 125

Sabotta, S., Schlecker, M., Chaturvedi, P., et al. 2021, A&A, 653, A114

Sameshima, H., Matsunaga, N., Kobayashi, N., et al. 2018, PASP, 130, 074502

Sebold, S., & Perktold, J. 2010, in 9th Python in Science Conference

Skumanich, A. 1972, ApJ, 171, 565

Snellen, I., Albrecht, S., Anglada-Escude, G., et al. 2019, arXiv e-prints [arXiv:1908.01803]

Suárez Mascareño, A., González-Álvarez, E., Zapatero Osorio, M. R., et al. 2023, A&A, 670, A5

The Astropy Collaboration 2018, astropy v3.1: a core python package for astronomy, Zenodo

Thibault, S., Brousseau, D., Poulin-Girard, A.-S., et al. 2022, SPIE Conf. Ser., 12184, 1218454

Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A&A, 612, A86

Turbet, M., Fauchez, T. J., Leconte, J., et al. 2023, A&A, 679, A126

Ulmer-Moll, S., Figueira, P., Neal, J. J., Santos, N. C., & Bonnefoy, M. 2019, A&A, 621, A79

Wang, S. X., Latouf, N., Plavchan, P., et al. 2022, AJ, 164, 211

Wolf, E. T., Haqq-Misra, J. D., Kopparapu, R., et al. 2020, in AGU Fall Meeting Abstracts, 2020, P020-01

Wordsworth, R. D., Forget, F., Selsis, F., et al. 2011, ApJ, 733, L48

Yang, J., Boué, G., Fabrycky, D. C., & Abbot, D. S. 2014, ApJ, 787, L2

Yang, J., Leconte, J., Wolf, E. T., et al. 2016, ApJ, 826, 222

Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577

Zechmeister, M., Reiners, A., Amado, P. J., et al. 2018, A&A, 609, A12
## Appendix A: The wapiti method step-by-step

In this appendix we present the wapiti method in detail, outlining each step in Fig. A.1 as applied to the RV data of Gl 725 B and a tutorial can be found online<fnref n="9" />. In summary, the wapiti method applies a wPCA to the LBL per-line RV time series, with the goal of finding the principal vectors that best serve as proxies for the systematic errors in the RV data.

Before applying wPCA to the data, per-line RV time series with more than 50% of missing data are excluded. Subsequently, the average RV time series is subtracted from each per-line RV time series to ensure that Keplerian or any achromatic signals, such as trends, remain unaffected by wPCA and are not misinterpreted as systematic errors. Finally, the time series are standardized by using their inverse-variance weighted mean and variance (Ould-Elhkim et al. 2023).

The wapiti method selects the relevant components by evaluating their contribution to the RV fit. Given a baseline model (e.g., a simple constant offset), we compute the BIC for this model and then the BIC after adding each principal component individually. The log BF comparing these two models allows us to quantify the improvement in fit. The robustness of the log BF value is assessed through a permutation test, in which each column of the RV component data is permuted, effectively shuffling the data to test the stability of the log BF under random configurations.

Specifically, for each permutation, a new wPCA model is fitted to the permuted data, and the log BFs are recalculated for each component. This process is repeated 1000 times, and the resulting distribution of log BFs is examined to determine the significance of the original fit for each component using a p-value test. A component is selected if its log BF value has a p-value below $10^{-3}$. In the case of Gl 725 B only the 1st, 2nd, 7th and 9th components are selected (see Fig. A.1).

Subsequently, we rearrange the selected components in descending order of their BIC values. During this process, we identify the component that, when included, minimizes the BIC, and designate it as the new primary component. We then iteratively determine the subsequent components: at each step, we add the next component that, when combined with the already selected components, further minimizes the BIC. This procedure produces a reordered sequence of principal vectors, denoted $\hat{V}_i$. Using the mathematical notation of Cretignier et al. (2023), this approach corresponds to calculating the BIC when fitting the following model to the data:

$$RV(t) = \gamma_{SPIRou} + trend \cdot t + \sum_{i=1}^{N} a_i V_i. \quad (A.1)$$

Here, we include an offset term $\gamma_{SPIRou}$, a trend term (due to the binary companion Gl 725 A), the epochs of observation $t$ and the coefficients $a_i$ to fit the principal vectors $V_i$ in the model. The parameter $N$ represents the number of components used in the reconstruction, and can be set to 0, corresponding to the data without applying the wapiti correction. We fit these models using weighted least-squares regression with the statsmodels module (Seabold & Perktold 2010), which also provided their BIC values.

This reordering method is similar to the $\chi^2$ reordering technique employed in Collier Cameron et al. (2021), the key difference being our use of the BIC as the criterion to penalize the inclusion of too many components. For Gl 725 B, the components in BIC-reordered order are the 2nd, 9th, 1st and 7th. To

<figure><img src="image_9.png" /><figcaption>**Fig. A.1.** Top: log BF of each individual wPCA component fitted separately to the RV data. Components that pass the permutation test with a p-value below $10^{-3}$ are highlighted in red. Middle: BIC values before and after reordering the significant components. Bottom: gain in model evidence associated with each additional component and the red dashed line marks the threshold of 5.</figcaption></figure>

determine how many components to include, we add the BIC-reordered components one by one and monitor the gain in model quality by computing the change in BIC between consecutive models. We stop adding components once the log BF gain per added component drops below a threshold of 5 (see Fig. A.1).

<footnote n="9">⁹ https://github.com/HkmMerwan/wapiti</footnote>
## Appendix B: Physical origin of the first principal vector in the case of telluric contamination

To understand the origin of the first principal vector $V_1$, we calculated what we refer to as an indicator of telluric contamination, denoted as $d1_⊕$. This indicator is determined by computing the Pearson correlation for all lines between the telluric template and the first derivative of the star template. As a reminder, in the LBL algorithm, the RV is determined by subtracting the star template from the observed spectrum and then projecting the residuals onto the first derivative of the star template (Artigau et al. 2022). Therefore, the use of Pearson correlation is a natural choice since it allows us to investigate the linear correlation of this derivative with the telluric lines. By averaging these individual time series, we derive the telluric indicator $d1_⊕$.

However, calculating this correlation for all lines would be unnecessary since not all lines contribute to the systematic effects with the same importance. To identify the lines most responsible for the systematics, we apply the method introduced in Ould-Elhkim et al. (2023), which consists in computing the z-scores $z_i$ for each component $i$ among all lines. Using the mathematical framework from Cretignier et al. (2023) this means performing a linear fit of each spectral line $j$ to the principal vectors $V_i$ and obtaining the corresponding coefficients $a_{i,j}$ such as

$$RV(t, \lambda_j) = \gamma_j + \sum_i a_{i,j} V_i(t), \qquad (B.1)$$

where $\gamma_j$ is a constant offset. Using these coefficients, we can then compute the z-scores $z_{i,j}$ of each line $j$ for each component $i$, defined as $a_{i,j} = \mu_i / \sigma_i$, where $\mu_i$ and $\sigma_i$ represent the weighted mean and standard deviation of the coefficients across the spectral lines.

Consequently, in the two cases of a narrow and wide BERV coverage, we computed $d1_⊕$ by only calculating the Pearson coefficients for the top 1% of lines based on the absolute $z_1$ values. We chose the top 1% lines to fasten the computation but using other values such as 5 or 10% had no impact on the conclusion of our results. Figure B.1 shows that $d1_⊕$ and $V_1$ exhibit a similar pattern and are, in fact, highly correlated with a Pearson correlation coefficient of $r = 0.98$ and $r = 0.86$ for, respectively, the narrow and wide BERV coverage cases. In conclusion, our interpretation of $V_1$ from these simulations is that it serves as a proxy to assess the impact of telluric absorption lines on the first derivative of the template.

## Appendix C: Telluric contamination of stellar activity indicators

Our simulations also show that stellar activity indicators can be contaminated by tellurics. This is notably the case for the dLW (Zechmeister et al. 2018): the dLW quantifies alterations in the average width of line profiles in relation to the template (Artigau et al. 2022). In our simulations, we observed that at the 3% level of telluric contamination, this indicator is affected in both a wide and narrow BERV coverage (see Fig. C.1). This is also evident in the real observations of Gl 725 B, which show similar contamination (Fig. C.2). This could pose challenges when attempting to identify and characterize stellar activity within the data, as it compromises the reliability of dLW as a stellar activity indicator.

Moreover, it is possible to extend the application of wPCA to the per-line dLWs, resulting in the computation of the first

principal vector, denoted as $W_1$. This principal vector has shown promise as a stellar activity indicator, as demonstrated when studying the young active M dwarf AU Mic (Donati et al. 2023a). However, $W_1$ also appears to be susceptible to telluric contamination, as seen in both our simulations (Fig. C.1) and real observations (Fig. C.2). Thus, the $W_1$ indicator is also vulnerable to tellurics and caution should be exercised when using it for activity characterization.

If stellar activity creates variability in the per-lines dLWs, it would be more advisable to check for activity signals in all the relevant principal vectors $W_i$ instead of just $W_1$. In the case of AU Mic, its success could be attributed to the high activity level of this young M dwarf, surpassing the influence of tellurics and contributing to the indicator's effectiveness by being responsible for the most variability in the data. Nevertheless, the impact of telluric lines on $W_1$ has been shown to be manageable through a filtering technique that effectively removes telluric contamination, making it a reliable proxy for the small-scale magnetic field (Charpentier et al. 2025).

## Appendix D: Supplementary material regarding the RV analysis

### Appendix D.1: Corner plots of the posterior distributions

See Figs. D.1 and D.2.

### Appendix D.2: Statistical significance of the planetary signals

To assess the reliability of the two signals detected at 4.77 d and 37.9 d, we computed a FIP periodogram (Hara et al. 2022) from the samples produced by the UltraNest sampler (Buchner 2021) used to estimate the Bayesian evidences of the $n_{max} = 0, 1, 2, 3$-planet models.

#### 1. Sampler and priors. We employ ReactiveNested Sampler from UltraNest, with a mixed-random Slice Sampler step, to explore the joint posterior of the parameters from Eq. (6):

$$\{\gamma, \text{trend}, (a_i)_{i=1...n_{\text{sys}}}, (A_j, B_j, P_j)_{j=1...n_{\text{max}}}\}.$$

The prior distributions are:

- $\gamma \sim \mathcal{U}(500, 1500) \text{ m s}^{-1}$

- trend $\sim \mathcal{U}(-0.1, 0.1) \text{ m s}^{-1} \text{ d}^{-1}$,

- $a_i \sim \mathcal{U}(-100, 100) \text{ m s}^{-1}$,

- $A_j, B_j \sim N(0, 5) \text{ m s}^{-1}$

- $P_j \sim \text{LogUniform}(1.5, 50) \text{ d}$.

For each $n_{max} \in \{0, 1, 2, 3\}$, we perform three independent runs with 400 live points each, and then pool the weighted samples.

#### 2. FIP computation and results (Fig. D.3)

Details on the FIP computation can be found in a tutorial online<fnref n="10" />, we plotted $-\log_{10}$ FIP (blue) and $\log_{10}$ TIP (red) as a function of period for $n_{max} = 1, 2, 3$.

- The posterior probability to have $n_{max} = 0$ planets is of $5 \times 10^{-79}$

- The posterior probability to have $n_{max} = 1$ planets is of $5 \times 10^{-45}$, only the 37.9 d signal appears and hints of the 4.77 d signal appears in the TIP.

- The posterior probability to have $n_{max} = 2$ planets is of $1 \times 10^{-20}$, two clear peaks emerge at 4.77 d and 37.9 d.

- The posterior probability to have $n_{max} = 3$ planets is of 1.

<footnote n="10">¹⁰ https://github.com/nathanchara/FIP/blob/main/fip_computation.ipynb</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. B.1.** Relationship between the first principal component $V_1$ and the telluric contamination indicator $d1_⊕$ in the simulated data. **Top:** How $V_1$ and $d1_⊕$ vary with BERV for narrow and wide coverage. **Bottom:** Periodogram of $V_1$ and the correlation between $V_1$ and $d1_⊕$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. C.1.** Telluric contamination of stellar activity indicators of the simulated data in the wide (green) and narrow (red) BERV coverage cases. **Left:** dLW versus the BERV and the respective periodograms below. **Right:** $W_1$ time series versus the BERV and the respective periodograms below.</figcaption></figure>

## Appendix D.3: Sensitivity of the detected signals to the choice of NZP

To assess the robustness of the periodicities against instrumental long-term drifts, we repeated the signal search over a grid of NZP time series. We considered 16 NZP variants and four SNR

thresholds applied to the time series (0, 100, 120, 150), yielding $16 \times 4 = 64$ NZP×SNR realizations of the RV dataset. For each realization, we ran the iterative Bayesian periodogram (Sect. 5), and we recorded the best period and its log BF at iteration 1. Whenever iteration 1 crossed the detection threshold (log BF >
<figure><img src="image_12.png" /><figcaption>**Fig. C.2.** Telluric contamination of stellar activity indicators of Gl 725 B. Left: dLW versus the BERV and the respective periodograms (below). Right: $W_1$ time series versus the BERV and the respective periodograms.</figcaption></figure>

5), we proceeded to iteration 2 and repeated the same procedure on the residuals.

For both iterations, we grouped the returned periods into clusters using a 2% relative tolerance in period and we computed summary statistics per cluster (consensus period, multiplicity across realizations, and evidence metrics). Below we only quote the clusters relevant to the planetary interpretation.

Iteration 1 (before adding any Keplerian) Across all NZP×SNR realizations, the most frequent first-iteration solution is a cluster centred on $P = 37.894$ d with a narrow spread (period range [37.585, 38.051] d) and strong evidence ($\log BF_{mean} = 8.05$, $\log BF_{median} = 8.30$, $\log BF_{max} = 12.49$; $n = 25$ realizations). Other first-iteration clusters appear less coherent and/or with lower evidence (e.g. a 13.80 d cluster with $\log BF_{mean} \approx 0.70$, $n = 19$; a 92.40 d cluster with $\log BF_{mean} \approx 5.19$, $n = 5$). A handful of realizations return a ~4.77 d peak already at iteration 1 (period range [4.763, 4.771] d, $\log BF_{max} = 8.41$), but with low average support ($\log BF_{mean} \approx -0.21$, $n = 5$), indicating that this short-period solution is not consistently preferred before accounting for the ~38 d signal.

Iteration 2 (conditioned on the iteration 1 detection). Among the realizations where iteration 1 was significant (39 in total), the iteration 2 search most frequently returns a cluster at $P = 4.770$ d (period range [4.763, 4.771] d) with moderate evidence on average ($\log BF_{mean} = 1.49$, $\log BF_{median} = 1.29$, $\log BF_{max} = 6.76$; $n = 23$). A weaker ~38 d family is sometimes re-selected at iteration 2 ($\log BF_{mean} \approx 2.00$, $n = 4$), as are seasonal/rotational harmonics near 90–135 d ($\log BF_{mean} \approx 4.33$ at ~91 d; $n = 4$).

**Summary.** The $P \sim 38$ d signal is (i) the most common outcome at iteration 1 across NZP×SNR choices, and (ii) consistently associated with strong $\log BF$ values, with a tight inter-realization period dispersion. This points to a high stability of the ~38 d solution with respect to the NZP modelling and supports its planetary nature.

Conversely, the $P \sim 4.77$ d signal becomes the most frequent solution after the ~38 d component is included, but its average evidence at iteration 2 remains modest (though it reaches $\log BF \approx 6.8$ in the best case). The moderate mean $\log BF$ and the residual sensitivity to the NZP choice argue for a conservative classification of the ~4.77 d signal as a candidate.
<figure><img src="image_13.png" /><figcaption>**Fig. D.1.** Corner plot of the posterior distributions for instrumental, stellar activity, and systematics parameters, including SPIRou offsets, jitter terms, wapiti coefficients, and hyperparameters such as rotation period ($P_{\text{rot}}$), timescale ($\rho$), and decay term ($\eta$).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. D.2.** Corner plot of the posterior distributions for the orbital parameters of planets b and c, including orbital period, phase ($\lambda_0$) and semi-amplitude ($K$).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. D.4.** log BF of the strongest periodic signal as a function of period, for all tested combinations of NZP realizations and SNR thresholds. **Top:** Results from the first iteration (before adding any Keplerian component). **Bottom:** Second iteration, i.e. search in the residuals after the first detected signal is included. A horizontal dashed line marks the detection threshold at log BF = 5.</figcaption></figure>

<caption>**Fig. D.3.** FIP periodograms of Gl 725 B for models with up to *n*<sub>max</sub> = 1, 2, 3 planets (top to bottom).</caption>