# Line-by-line Velocity Measurements: an
Outlier-resistant Method for Precision Velocimetry

Étienne Artigau, Charles Cadieux, Neil Cook, René Doyon, Thomas Vandal, Jean-François Donati, Claire Moutou, Xavier Delfosse, Pascal Fouqué, Eder Martioli, et al.

► To cite this version:

Étienne Artigau, Charles Cadieux, Neil Cook, René Doyon, Thomas Vandal, et al. Line-by-line Velocity Measurements: an Outlier-resistant Method for Precision Velocimetry. The Astronomical Journal, 2022, 164 (3), pp.84. 10.3847/1538-3881/ac7ce6. hal-03827457

HAL Id: hal-03827457

https://hal.science/hal-03827457

Submitted on 25 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Line-by-line Velocity Measurements: an Outlier-resistant Method for Precision Velocimetry

Étienne Artigau<fnref n="1" />,<fnref n="2" /> ID, Charles Cadieux<fnref n="1" /> ID, Neil J. Cook<fnref n="1" /> ID, René Doyon<fnref n="1" />,<fnref n="2" /> ID, Thomas Vandal<fnref n="1" /> ID, Jean-François Donati<fnref n="3" /> ID, Claire Moutou<fnref n="3" /> ID, Xavier Delfosse<fnref n="4" /> ID, Pascal Fouque<fnref n="3" />,<fnref n="5" /> ID, Eder Martioli<fnref n="6" />,<fnref n="7" /> ID, François Bouchy<fnref n="8" /> ID, Jasmine Parsons<fnref n="1" />,<fnref n="9" /> ID, Andres Carmona<fnref n="4" /> ID, Xavier Dumusque<fnref n="8" /> ID, Nicola Astudillo-Defru<fnref n="10" /> ID, Xavier Bonfils<fnref n="4" /> ID, and Lucille Mignon<fnref n="4" /> ID

Received 2022 February 18; revised 2022 June 23; accepted 2022 June 27; published 2022 August 8

## Abstract

We present a new algorithm for precision radial velocity (pRV) measurements, a line-by-line (LBL) approach designed to handle outlying spectral information in a simple but efficient manner. The effectiveness of the LBL method is demonstrated on two data sets, one obtained with SPIRou on Barnard's star, and the other with the High Accuracy Radial velocity Planet Searcher (HARPS) on Proxima Centauri. In the near-infrared, the LBL provides a framework for meters-per-second-level accuracy in pRV measurements despite the challenges associated with telluric absorption and sky emission lines. We confirm with SPIRou measurements spanning 2.7 yr that the candidate super-Earth on a 233 day orbit around Barnard's star is an artifact due to a combination of time sampling and activity. The LBL analysis of the Proxima Centauri HARPS post-upgrade data alone easily recovers the Proxima b signal and also provides a 2σ detection of the recently confirmed 5 day Proxima d planet, but argues against the presence of the candidate Proxima c with a period of 1900 days. We provide evidence that the Proxima c signal is associated with small, unaccounted systematic effects affecting the HARPS-TERRA template-matching radial velocity extraction method for long-period signals. Finally, the LBL framework provides a very effective activity indicator, akin to the FWHM derived from the cross-correlation function, from which we infer a rotation period of $92.1_{-3.5}^{+4.2}$ days for Proxima.

Unified Astronomy Thesaurus concepts: Exoplanets (498); Radial velocity (1332); Astronomy data analysis (1858)

## 1. Introduction

Within the last three decades, the search for planets orbiting stars other than our Sun has moved from a speculative endeavor to one of the most active fields in astronomy. The first successful technique to uncover exoplanets around stars was the radial velocity (RV) method (Lindgren & Dravins 2003), where one measures the reflex motion of a star as it is orbited by a low-mass companion. While that technique has been superseded in the raw number of discoveries by the transit method, RV measurements remain central in exoplanet studies for providing mass measurements. Since the watershed discovery of the planet orbiting 51 Pegasi, the state-of-the-art accuracy achieved by this technique has been lowered from ~15 m s⁻<fnref n="1" /> (Mayor & Queloz 1995) to better than 20 cm s⁻<fnref n="1" /> (Faria et al. 2022).

While early precision radial velocity (pRV) efforts focused on Sun-like stars, M dwarfs are now the topic of active research. Having a relatively low mass (0.078 to ~0.6 M☉; Baraffe et al. 2015) and small radius (0.1–0.6 R☉), a given planet will impart a larger RV signal and a much deeper transit than if it were orbiting a Sun-like star. When focusing on planets located in the habitable

zone, this effect is further enhanced: habitable zones being more compact for M dwarfs than for Sun-like stars. This allows for meters-per-second-level accuracy measurements to detect Earth-mass planets in the habitable zone of their parent stars (e.g., Bonfils et al. 2013; Delfosse et al. 2013; Anglada-Escudé et al. 2016), a feat yet to be achieved for such a planet around a Sun-like star, for which the RV signal would be at the level of 10 cm s⁻<fnref n="1" />.

The extension of pRV to late-type stars, combined with an interest in transit spectroscopy, which displays numerous atmospheric spectroscopic features in the near-infrared (e.g., He at 1.083 µm to the CO bands starting at 2.29 µm), has led to the development of pRV spectrographs covering the *YJHK* bandpasses (Artigau et al. 2012; Mahadevan et al. 2012; Kotani et al. 2014; Claudi et al. 2017; Quirrenbach et al. 2018).

Precision spectroscopy in the near-infrared (NIR) is significantly more challenging than in the optical domain for a number of reasons. First, telluric absorption bands are much more numerous in the *YJHK* domain, adding an imprint at near-zero velocity on all spectroscopic observations. Left uncorrected, biases due to telluric absorption would completely overwhelm planetary signals. A number of approaches have been proposed to correct telluric absorption in the NIR using both model-based (e.g., Gullikson 2014; Smette et al. 2015) and data-based methods (Artigau et al. 2014; Bedell et al. 2019; Cretignier et al. 2021), or the observation of hot stars shortly before or after the

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Université de Montréal, Département de Physique, IREX, Montréal, QC, H3C 3J7, Canada; etienne.artigau@umontreal.ca</footnote>

<footnote n="2">² Observatoire du Mont-Mégantic, Université de Montréal, Montréal, QC, H3C 3J7, Canada</footnote>

<footnote n="3">³ Université de Toulouse, CNRS, IRAP, 14 Avenue Belin, F-31400 Toulouse, France</footnote>

<footnote n="4">⁴ Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="5">⁵ Canada-France-Hawaii Telescope, CNRS, Kamuela, HI 96743, USA</footnote>

<footnote n="6">⁶ Laboratório Nacional de Astrofísica, Rua Estados Unidos 154, Itajubá, MG 37504-364, Brazil</footnote>

<footnote n="7">⁷ Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="8">⁸ Département d'astronomie, Université de Genève, Chemin des Maillettes 51, CH-1290 Versoix, Switzerland</footnote>

<footnote n="9">⁹ Department of Physics & McGill Space Institute, McGill University, 3600 Rue University, Montréal, QC, H3A 2T8, Canada</footnote>

<footnote n="10">¹⁰ Departamento de Matemática y Física Aplicadas, Universidad Católica de la Santísima Concepción, Alonso de Rivera 2850, Concepción, Chile</footnote>
science observation (Vacca et al. 2003). The NIR night sky also displays strong emission features due mostly to OH airglow (Rousselot et al. 2000).

At the typical resolution of pRV spectrographs ($\frac{\lambda}{\Delta\lambda} \sim 50,000-100,000$), emission lines occupy only a small fraction of resolution elements (typically percent level), but are strong enough to bias the RV measurement, especially for stars near the ecliptic pole with small barycentric radial velocity (BERV). Other effects that are specific to NIR spectrographs involve the significantly worse detector cosmetics compared to optical spectrographs. The state-of-the-art detectors used in NIR pRV spectrographs are Hawaii-4RG or Hawaii-2RG mosaics, which have dark current and cosmetics (e.g., bad pixel clusters, hot pixel clusters, dead columns, persistence) that lead to fixed structures in the extracted spectra which may bias pRV observations (Artigau et al. 2018a). Last but not least, NIR pRV requires a cryogenic spectrograph with thermal control at a level of 1 mK, which in itself has been a major engineering challenge for all instruments with similar science goals (e.g., Mirabet et al. 2014; Stefansson et al. 2016; Donati et al. 2020).

This paper presents a new method for RV measurements that is tailored to spectra having spurious excursions in flux from unaccounted for effects: the line-by-line (LBL; see Appendix A for code availability) pRV method. The technique has been proposed by Dumusque (2018) and further explored by Cretignier et al. (2020) in the context of optical pRV observations and mostly tailored toward filtering stellar activity in solar-type stars. Here, the concept is expanded to include the full RV content of molecular bands of M dwarfs and account for the needs of NIR pRV measurements where resilience against outliers is key. Section 2 provides the mathematical background pertaining to the discussion as well as a brief overview of other pRV numerical techniques. Section 3 presents the data sets used for the demonstration of the LBL method. The mathematical framework of the LBL algorithm is described in Section 4. Section 5 presents a demonstration of the LBL technique based on SPIRou and High Accuracy Radial velocity Planet Searcher (HARPS) data sets on Barnard's star (Barnard 1918) and Proxima Centauri (Innes 1915), respectively. We conclude with Section 6.

## 2. Precision Radial Velocity Measurement Methods

### 2.1. Bouchy et al.'s Mathematical Framework

Bouchy et al. (2001; hereafter B01) presents a framework that can be used to reduce the problem of pRV measurements in noisy data to simple arithmetic; this work expands on a formalism first introduced by Connes et al. (1996). For pixel $i$, if one subtracts a spectrum $A(i)$ with infinitesimal RV shift $\delta V$ (or, expressed as a fractional wavelength change, $\frac{\delta V}{c} = \frac{\delta \lambda}{\lambda}$) from a zero velocity template $A_0(i)$, the residual ought to be the derivative of the template times the velocity change. More explicitly, it corresponds to the derivative of the template with respect to the wavelength grid expressed as a velocity. In this context, Equation (2) in B01 is

$$A(i) - A_0(i) = \frac{\partial A_0(i)}{\partial \lambda(i)} \delta \lambda(i), \quad (1)$$

which can be expressed as a velocity change

$$\frac{\delta V(i)}{c} = \frac{A(i) - A_0(i)}{\lambda(i) \frac{\partial A_0(i)}{\partial \lambda(i)}}. \quad (2)$$

With an estimate of the local noise $\sigma(i)$ ($[A(i) - A_0(i)]_{\text{RMS}}$ in B01), which we obtain through a running standard deviation estimate of the $A(i) - A_0(i)$, one can determine the per-pixel RV accuracy $\sigma_V(i)$. In our implementation, the noise estimate is computed with a running box the size of 500 km s$^{-1}$ (same size as the high-pass filtering width; see Section 4.10) by taking half the difference between the 84th and 16th percentiles of residuals of $A(i) - A_0(i)$. This departs slightly from the B01 and Dumusque (2018) implementation of the method as we do not explicitly assume that the noise is solely due to the combined contributions of Poisson and readout noise in the spectrum. An empirical determination of the noise relative to the template may be important when the detector readout-noise contribution is not straightforward to determine (e.g., spatially correlated readout noise in IR arrays, integration-time-dependent readout noise in up-the-ramp readout schemes) or when modal noise for optical elements is a significant contribution (Blind et al. 2017). With the definition above for $\sigma(i)$, Equation (2) becomes

$$\sigma_V(i) = \frac{c \sigma(i)}{\lambda(i) \frac{\partial A_0(i)}{\partial \lambda(i)}}. \quad (3)$$

Considering that an optimal weighted sum is performed with weights proportional to the inverse of the square of uncertainties, one therefore has a mean velocity for the pixels of interests. That is,

$$\delta V = \frac{\sum \delta V(i) \sigma_V(i)^{-2}}{\sum \sigma_V(i)^{-2}}. \quad (4)$$

The uncertainty of a weighted sum becomes

$$\sigma_V = \frac{1}{\sqrt{\sum \sigma_V(i)^{-2}}}. \quad (5)$$

The B01 formalism is, by construction, optimal in terms of signal-to-noise (S/N) and provides a limit that cannot be surpassed in terms of RV accuracy for a given spectrum and noise.

A key point that will be used in the following sections is that summations in Equations (4) and (5) can be performed over an interval in wavelength that is of arbitrary length.

### 2.2. The Cross-correlation Function Method

One of the earliest methods to derive precise RVs from high-resolution spectroscopy is the cross-correlation function (CCF). The CCF method has its historical roots in pre-CCD spectroscopy (Baranne et al. 1979), where a mechanical mask was superposed onto a cross-dispersed spectrum and the total flux passing through holes corresponding to line positions was integrated with a single-pixel photoelectric sensor. By moving the mask mechanically, one could directly record the mean line profile. In modern spectrographs, the spectrum is cross correlated with a comb of delta functions of varying strength matching known stellar lines (Baranne et al. 1996; Pepe et al. 2002; Lafarga et al. 2020). One then derives a weighted mean line profile, the weight of each line scaling as its depth and
local S/N. The mean line S/N is far larger than that of the input spectrum and scales as the square root of the number of lines. The CCF provides a central velocity measurement and the higher moments of the line profiles inform on stellar activity (e.g., Santos et al. 2001). One significant advantage of the CCF compared to other methods is that it provides a velocity measurement with a single spectrum, assuming that one has a mask for the relevant spectral type. It also provides more complex information on line structures, such as the presence of double lines in a spectroscopic binary, that are not explicitly retrieved by other pRV measurement methods.

As detailed in Section 4.3, pRV in the NIR domain is strongly hampered by telluric absorption residuals and its time-varying impact on a given line profile. These effects limit the usefulness of the CCF method at these wavelengths.

### 2.3. Template Matching

While simple, robust, and powerful, the CCF technique does not fully exploit the RV content (Anglada-Escudé & Butler 2012). A number of groups have devised numerical methods where a high-S/N template is adjusted to a given spectrum to measure its velocity. Given the proper weights, as described in Section 2.1, this method can fully exploit the RV content of the observation.

Although the numerical details vary between implementations, template matching is now the most widely used method for pRV measurements of late-type stars with rich molecular bands. Virtually all recent measurements of low-mass planets around M dwarfs have employed this type of technique. The main template-matching implementations widely used in the literature are HARPS-TERRA (Anglada-Escudé & Butler 2012), NAIRA (Astaridillo-Defru et al. 2017), WOBBLE (Bedell et al. 2019), and SERVAL (Zechmeister et al. 2018).

## 3. Demonstration Data Sets

In describing the LBL pRV measurement algorithms below, one needs sample data sets to perform the demonstration. As the initial motivation behind our implementation of the LBL method is to improve pRV measurements' resilience to outliers in IR observations, we used SPIRou data on Barnard's star, known to be stable within a few meters per second. To compare the performances of the LBL method in a domain with few telluric lines with a publicly available benchmark analysis, we also performed LBL analysis using HARPS data on Proxima Centauri.

### 3.1. SPIRou Data on Barnard's Star

Observations on Barnard's star were obtained with SPIRou as part of the SPIRou Legacy Survey (Donati et al. 2020) between 2018 September and 2021 October. The typical Barnard's star observation has a per-pixel S/N of 100, 150, 200, and 200 in the Y, J, H, and K photometric bandpasses, respectively.

The SPIRou data presented here were reduced with version 0.7.194 of APERO (A PipelinE to Reduce Observations; N. J. Cook et al. 2022, in preparation). APERO performs all the steps from the raw 2D science frames to science-ready, wavelength-calibrated data. The main data reduction steps are as follows. APERO starts with science frames that are constructed from the detector control software; these 4096 × 4096 frames represent the per-pixel accumulation rate of incoming flux (ramp images) and are corrected for nonlinearity. APERO performs a set of quality checks on these frames to identify spurious pixels and filters

correlated noise between amplifiers before the more standard extraction procedures. Frames are registered to a common position to account for slow bulk shifts/rotations of the science array, and the 49 diffraction orders are identified and extracted optimally using the Horne (1986) formalism while accounting for the shape of the SPIRou pupil slicer (Micheau et al. 2018). The product of the APERO extraction is E2DS (4088 × 49 pixels), an array consisting of all extracted spectral orders. Once science data and calibrations have been extracted, wavelength calibration is performed similarly to the scheme described in Hobson et al. (2021). As SPIRou is a spectropolarimeter with two science fibers (A and B) as well as a calibration fiber (C), one can extract per-fiber data for spectropolarimetric work or combined AB extraction for intensity-only work. The calibration fiber allows one to observe a reference Fabry-Pérot simultaneously to science observations and calibrate intranight drifts; on typical nights, these drifts are at the <2 m s⁻¹ level. In the present work, we used combined AB extractions and their corresponding C fiber drift files. In addition to accurate wavelength calibration, it is essential to remove telluric absorption signatures to minimize their imprint in science data. APERO handles telluric absorption in a two-step process. A first-order correction is done by adjusting a TAPAS (Transmissions of the Atmosphere for Astrophysical data; Bertaux et al. 2014) absorption model that only leaves the H₂O and dry absorption (combining CH₄, O₂, CO₂, N₂O, and O₃) optical depths as free parameters. This leaves percent-level residuals for deep (>50%) absorption features and has the advantage of simplicity. This pre-cleaning is also performed on a set of rapidly rotating hot stars, from which we derive a library of residuals that are correlated against H₂O and dry absorption optical depth to create a residual model with 3° of freedom per pixel (one term proportional each to the optical depth of H₂O and dry components and a constant). This telluric residual model is subtracted from science frames after the TAPAS pre-cleaning. This technique leaves comparable residuals to principal-component-analysis-based methods (Artigau et al. 2014), but has the advantage of being applied locally; spurious pixels in either the hot star or science observations only affect a handful pixels rather than the entire fit. A representative telluric-cleaned spectrum before and after telluric correction is shown in Figure 1. Once pRV measurements have been performed, a two-step correction of the RV measurements is necessary. One first accounts for the drifts as measured with the C fiber simultaneously with science observations to account for intranight drifts. In addition, a zero-point correction is applied to the data from the observation of a set of RV-stable calibrators, similar to what is done with the SOPHIE (Spectrographe pour l'Observation des Phénomènes des Intérieurs stellaires et des Exoplanètes) spectrograph (Courcol et al. 2015). The details of zero-point correction will be presented in a future contribution (T. Vandal et al. 2022, in preparation). The data set comprises 871 observations. The median per-point error is 2.57 m s⁻¹ with a standard deviation of 2.59 m s⁻¹. We did not bin the individual frames, keeping the four exposures within each sequence as individual measurements in the reported time series.

### 3.2. High Accuracy Radial velocity Planet Searcher Data on Proxima

We used HARPS (Pepe et al. 2002) publicly available data on Proxima Centauri to assess the performance of the LBL technique in the optical domain. We retrieved the 2D (per-order) spectra and CCF data (M2V mask) from the ESO science archive (Delmotte et al. 2006), and concentrated only on epochs
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Sample spectrum of Barnard's star over a short-domain wavelength within the *H* band. The telluric-corrected observation is overplotted with the template registered to that star's line-of-sight velocity accounting for the barycentric radial velocity at the time of observation.</figcaption></figure>

following the fiber upgrade (Lo Curto et al. 2015), that is, after 2015 May 29. A simple telluric correction was performed by adjusting a TAPAS (Bertaux et al. 2014) model to the data. The LBL velocities are then compared in Section 5.2 to the same data analyzed with the HARPS-TERRA and NAIRA template-matching algorithms (Damasso et al. 2020; Astudillo-Defru et al. 2017). The FWHM information, used as an activity indicator and discussed in Section 5.3, were drawn from the CCF measurements provided by the ESO science archive.

# 4. The Line-by-line Method

## 4.1. Motivation

One of the main issues with any pRV measurement is the handling of spurious data lines in the spectrum that, if unaccounted for, can bias the measurement well above the typical uncertainties due to the white-noise contribution of the data. In the absence of spurious outliers, the RV measurement is conceptually straightforward: one obtains a high-S/N template from the ensemble of observations and performs the summations in Equations (4) and (5) to derive a velocity and the corresponding uncertainty. In essence, spurious regions in the spectrum are problematic as at some line-of-sight velocities they will randomly align with sharp spectroscopic features such as line wings (i.e., lines of large $|\partial A_0(i)/\partial\lambda(i)|$). Regardless of the method, one needs to identify and reject these spurious lines.

Precision radial velocity in the NIR is further complicated by a number of issues. First of all, the Earth's atmosphere displays numerous absorption bands that are superimposed on the science data and must be corrected. Irrespective of the telluric subtraction methods used, one expects an increased noise in the telluric-cleaned regions from the lower effective throughput, as well as systematic effects due to the limitations of models and unaccounted for effects in the calibration data. One can reject the domain surrounding deep absorption, but in order to preserve the same domain observed throughout the year, one needs to reject a domain that is up to ~64 km s⁻¹ larger than the line itself because of Earth's orbit, which adds a line-of-sight change of ±32 km s⁻¹ for ecliptic targets. The loss is exemplified in Reiners et al. (2010), which reports a 55% and 46% effective loss in wavelength domain due to telluric absorption in the *J* and *H* bands, respectively. As we show in Section 4.3, this loss can be reduced to nearly zero with the proper data-analysis framework.

## 4.2. The Method

The LBL method is conceptually rather simple: one performs the B01 summation on the spectra in order to derive a velocity in a way that is optimal considering the noise in the data. The subtlety being that this sum is done step-wise; one performs the sum on relatively small domains called “lines.” One then performs a consistency check between the velocity of all lines, and finally averages line velocities into a single per-spectrum velocity with the corresponding uncertainty, propagated from the per-line uncertainty. As shown in Dumusque (2018) and Cretignier et al. (2020), the LBL algorithm can be efficient in identifying and accounting for lines affected by stellar activity for Sun-like stars in the optical domain. Here, the two-step process allows one to perform a very efficient rejection of lines polluted by a number of artifacts, both known (telluric absorption artifacts, detector defects, sky emission lines) and unknown effects. The LBL method is not immune to other effects such as the propagation of uncertainties in the wavelength solution. This type of uncertainty is different in nature from RV uncertainties induced by artifacts in the spectra and will affect any pRV method equally.

## 4.3. Gains from Domain Splitting

When obtaining a single RV measurement for a given spectrum, it is problematic to determine whether the measurement has been significantly biased by unknown residuals, either from uncertainties in telluric absorption correction, sky emission subtraction, or unknown effects. As RV effects from planets are expected to be achromatic, one can divide the spectrum into a number of subdomains, each having its own velocity measurement, and verify the statistical consistency between these measurements. As a thought experiment, consider a spectrum with a 1 m s⁻¹ theoretical RV accuracy, but having a single uncorrected glitch in the data that leads to a 2 m s⁻¹ error in the global RV measurement. If the spectrum is subdivided in *N* domains with equal RV content, the RV accuracy of each domain will be reduced by a factor $\sqrt{N}$, but the velocity error induced by the glitch will increase by a factor *N* within its segment. If the domain is split in *N* = 100 segments, the RV uncertainty for each segment will be 10 m s⁻¹, but the segment affected by the glitch will have its velocity offset by 200 m s⁻¹, or a 20σ outlier that can readily be identified. Furthermore, the cost in wavelength coverage of rejecting the segment with a spurious velocity will be small,
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Sample region of the Barnard's star template spectrum with the limits of *line* domains shown as dashed gray lines. Color-coding is arbitrary and highlights individual lines.</figcaption></figure>

scaling down as $N^{-1}$. There is an additional benefit with the increase in *N*: by setting a certain threshold in RV offset relative to all other segments, one sets an upper limit on the discrepancy of a given segment to all others in the spectrum, then the maximum global RV error injected from a glitch scales down as $1/\sqrt{N}$. In the toy model, a $5\sigma$ cutoff for segments would set an upper limit at $50\ \text{m s}^{-1}$ on any segment relative to others, which, in turn, cannot inject an error larger than $0.5\ \text{m s}^{-1}$ in the average RV. For a global RV accuracy of $\sigma_g$, $N$ segments and a rejection at the $n_{\sigma c}$ level, the largest error a glitch could inject in the global RV is $n_{\sigma c}\sigma_g/\sqrt{N}$. The idea of splitting a spectrum in numerous subdomains to identify outliers has also been used before in the context of iodine cell pRV measurements (see Section 6.5 in Hatzes 2019).

## 4.4. Mathematical Framework

For the purpose of the LBL algorithm, a line is simply defined as being the domain between two consecutive local maxima in the template spectrum. This has the advantage of simplicity, allowing the unambiguous empirical definition of lines in dense molecular bands. This method also benefits from being independent of resolution once stellar features are resolved. One can numerically differentiate the template, and find points having a zero derivative and negative second derivative (i.e., local maxima). The noise in the template should be significantly smaller than the point-to-point flux change in order to avoid identifying noise excursions as lines. Figure 2 illustrates the positions of lines for a representative region of the SPIRou Barnard's star spectrum template. To avoid having very narrow lines that may be spurious in a nearly flat region of the spectrum, we smooth the template with a boxcar kernel of 1 FWHM prior to finding the local maxima.

Given a template and a line domain, one computes a per-line velocity using Equation (4) and its corresponding uncertainty with Equation (5). One then derives the *velocity spectrum* (see Figure 3); in the case of Barnard's star, there are ~16,000 velocities with widely varying uncertainties; the best lines have uncertainty at the $50\ \text{m s}^{-1}$ level while the worst ones have uncertainties at the many kilometers per second level. The global RV precision for the sample spectrum discussed in Figures 1, 2, 3, and 4 is $2.53\ \text{m s}^{-1}$. There are a number of interesting features in the velocity spectrum. It is notable that a significant fraction of the RV content is retrieved from the telluric absorptions *outside* of the nominal *YJHK* bandpasses. The lower-middle panel of Figure 3 shows lines between *H* and *K* band, at ~1820 nm, that have a good RV accuracy (some at <300 $\text{m s}^{-1}$/line) without obvious outliers to the mean

velocity. Uncertainties are strongly chromatic, being significantly smaller in the *H* and *K* bands than blueward of 1400 nm. The fact that the best domain in the NIR lays longward of 1400 nm was documented by Artigau et al. (2018b) and Figueira et al. (2016). By performing a running weighted mean (see Appendix B) on the velocity spectrum over a $\Delta\lambda/\lambda = 5\%$, one can directly and precisely measure this RV content, as shown in Figure 4.

## 4.5. Iterations of the Line-by-line Algorithm

The B01 equations are only valid in the regime for which the higher-order derivatives are vanishingly small compared to the first-order derivative in Equation (1). This is analogous to the Newton algorithm to find the roots of a function by representing the function as a Taylor series expansion of the first order. This assumption is not always valid for pRV series, and one should use the LBL algorithm iteratively to satisfy this requirement. For a given observation, we spline the template $A_0$ to the expected systemic velocity of the star plus its BERV offset, derive $\delta V$ and $\sigma_V$ values, and, if $\delta V$ is not consistent with zero (in practice $|\delta V| < 0.1 \sigma_V$), we update the input velocity of $A_0$ and recompute $\delta V$ and $\sigma_V$. The LBL algorithm has a “capture velocity” that is about the FWHM of the effective line profile, which is typically $5\ \text{km s}^{-1}$ for mid-Ms in the NIR (Artigau et al. 2018b). The algorithm fails to converge if the input velocity is off by tens of kilometers per second, e.g., not accounting for the BERV in the initial splining of $A_0$. If the systemic velocity of the target is unknown, a first estimate should be determined beforehand, at the sub-FWHM level, with other numerical methods that can scan over large velocity ranges (e.g., CCF). For the LBL, the absolute velocity measurement of a template is described in Section 4.9. For RV-stable targets (with planetary and activity signals at the few meters per second level), the algorithm converges in one or two steps. To test the convergence in the presence of a large RV signal, we revisited the TOI-1278 data set presented in Artigau et al. (2021) where an $18\ M_{\text{Jup}}$ brown dwarf imparts a $K \sim 2300\ \text{m s}^{-1}$ RV signal on its M0V host. The RV excursions for the system are at about half of the FWHM of stellar lines. For TOI-1278, the algorithm typically converges in three steps, with a typical residual uncertainty that decreases by a factor ~30 between steps. The per-epoch RV uncertainty with the LBL is $9.4\ \text{m s}^{-1}$, which is a two-fold improvement compared to the previous published values, which were measured with a CCF (see Table 1 in Artigau et al. 2021). The TOI-1278 Keplerian fit values are consistent with those published albeit with significantly reduced uncertainties; this target is currently the subject of ongoing monitoring with SPIRou and its extended RV monitoring will be presented in a future publication.

## 4.6. Outliers

One key point in the LBL analysis is the filtering of problematic lines, and therefore the identification and abundance of outliers in the velocity spectra. Assuming that uncertainties are Gaussian in the input spectrum and properly accounted for, the B01 framework should lead to a normal distribution of uncertainties for all lines. In practice, we show that this is nearly the case out to $2\sigma$, with a small fraction of >5$\sigma$ outliers (see Figure B1). It is notable that outliers in the velocity spectrum do not specifically appear in telluric bands as one could expect (see Figures 3, B1, and Appendix B). This is
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Sample line-by-line velocity spectrum for Barnard's star showing the scatter of velocity measurements over a wide velocity range (±10 km s⁻¹ of mean velocity; upper panel); lines with uncertainties smaller than 300 m s⁻¹ are shown in dark blue to improve readability. Lines with uncertainties <300 m s⁻¹ are mostly concentrated within the *H* band and a large part of the *Y* and *J* bands are dominated by lines with kilometers-per-second-level radial velocity (RV) uncertainties. Lower panels show the same plot zoomed-in on one of the domains with the highest density of high-accuracy lines (1500–1550 nm; compare with RV accuracy in Figure 4), a domain strongly affected by telluric absorption but with significant RV content (1800–1850 nm), and the CO band head (past 2290 nm). Parts of the domain, such as the redder half of the *K* band and most of the *J* band, have very few, if any, high-quality lines. These regions correspond to domains where lines are either shallow (*K*) or broader (e.g., *J*; see Artigau et al. 2018b) than they are in the *H* band. This is also seen in Figure 4. These domains contribute very little to the overall RV budget. The large number of high-quality lines past 2290 nm (lower-left inset) correspond to the CO band heads.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Radial velocity (RV) accuracy achieved for a Δλ/λ = 5% running window (blue line). As demonstrated in Artigau et al. (2018b), the *H* band contains the domain with the highest RV content for the near-infrared domain. Remarkably, a single bin centered at 1.75 µm would allow a 5 m s⁻¹ RV accuracy, only a factor of ~2 worse than the full domain of SPIRou (2.53 m s⁻¹) for the observation considered. The full spectrum of Barnard's star is shown as an orange overplot.</figcaption></figure>

suggestive that the *APERO* telluric subtraction method is effective. The presence of a small number of outliers in the line velocities implies that one must use a form of sigma clipping. We adopt a simple mixture model (Bishop 2006; Hogg et al. 2010) to estimate the mean velocity, in which there are two

distributions, one for valid lines that follow a normal distribution corresponding to the uncertainties as derived from the B01 framework and having a common velocity, and the other, broader distribution corresponding to outliers. Appendix B provides the details on the numerical implementation. The resulting per-epoch velocity is $\bar{v}$ with the corresponding uncertainty $\bar{\sigma}$.

The LBL algorithm comes with a notable advantage relative to a number of other pRV methods regarding the RV content: it does not require the rejection of any wavelength domain to account for yearly BERV variations. In order to preserve the wavelength domain through the year, we accept that the domain will change with the yearly line-of-sight velocity variations and reject lines if they are found to be discrepant. This significantly increases the effective RV content of a given observation and basically sets us in the condition 3 of Figueira et al. (2016; only increase in variance due to the lower transmission is accounted for) rather than the very conservative condition 2 therein (rejection of all lines within ±30 km s⁻¹ of any absorption) and gets close to the RV accuracy that would be possible without any telluric absorption over the same domain (condition 1). Furthermore, unlike several template-matching methods, the LBL algorithm does not require any free parameters; there is no threshold in the equivalent of sigma clipping, this being handled in a statistically consistent way (see Appendix B).
## 4.7. Debiasing

One potential issue with the LBL method is that the actual list of valid lines will vary slightly from spectrum to spectrum. If a line has a mean velocity through the entire time series that is slightly offset from the velocity spectrum mean, then having it rejected at some epochs will lead to a bias. One needs to account for the varying line list. This step is done by constructing a 2D matrix, with *N* lines along one axis and *M* velocity spectra along the other axis; this matrix has gaps that vary between spectra. One needs to construct two 1D vectors, one along the *N*-axis and one along the *M*-axis, that correspond to the per-line bias and per-epoch velocity. Debiasing is done iteratively; one first computes a per-epoch mean velocity using the mixture model in Appendix B and then computes a per-line mean velocity after subtracting the per-epoch mean velocity. This allows one to determine the mean per-line velocity offset relative to the mean velocity of all lines. This offset is then subtracted from each line over the entire time series prior to computing the final $\bar{v}$ and $\bar{\sigma}$. This process could be iterated further, but in practice a single iteration is sufficient for velocity to converge well within $\bar{\sigma}$. The debiasing implies that the addition of new data will very slightly affect the RV measured in previous observations as the debiasing model will improve with time. In practical cases, this is a second-order effect that is smaller than the per-epoch RV uncertainty. The debiasing provides a significant improvement in the time series point-to-point scatter; without it, the rms of the Barnard's star sequence (2.59 m s⁻¹), presented in Section 5, increases by 50 cm s⁻¹.

## 4.8. Template Construction and Usage

The template spectrum is constructed from the median combination of BERV-registered spectra after telluric correction. The template construction is done with coadded, stitched spectral orders rather than per-order spectra. This allows for a proper accounting of wavelength overlap between consecutive orders. Figure 1 shows a small domain of the Barnard's star SPIRou spectrum, with the telluric correction applied to generate the telluric-corrected and template spectrum. Template creation is part of the *APERO* outputs and a standalone template creation is provided with the LBL tools (See Appendix A).

Considering that numerous lines and spectral features of interest have typical depths at the >5% level relative to the surrounding mean flux (a true continuum cannot be defined in the rich molecular bands of M dwarfs; see Figure 2), one needs to construct a template that has a sufficiently good S/N for identified lines to be true spectroscopic features rather than noise excursions. In practice, a template with S/N > 200 is sufficient to unambiguously identify spectral features and define “lines.” Considering that pRV time series typically have tens of observations, each with S/N > 50, this condition is generally met as the combined S/N of the templates scales as the square root of the number of observations. One further requirement in constructing the template is that any large RV shift from companions should be accounted for if the $\frac{\partial A_0(i) \lambda(i) \delta V}{\partial \lambda(i) c} \ll \sigma(i)$ condition is not met. This condition is generally met for planetary companions with signals below the ~100 m s⁻¹ level. If no template is available or a strong planetary signal is present, one can use a template from a star of a similar spectral type and rotational velocity. We computed the LBL sequence on Barnard's star detailed in Section 5.1 with the template constructed from GJ 1012 observations; both stars

are M4V and slow rotators, and the latter had a large number (449) of observations obtained through the SPIRou Legacy Survey (Donati et al. 2020). The resulting time series is similar to the one obtained with the Barnard's star's template, albeit with a ~20% worse point-to-point scatter. While we have not systematically explored the degradation of pRV performances for all subtype combinations, time series obtained with a model matching within ±2 stellar subtypes generally have performances similar to those obtained with an object's own template.

## 4.9. Absolute Velocity Measurement

The LBL algorithm is by nature differential as it measures a velocity difference between a template and a given science spectrum. It is often desirable to have an absolute velocity for a given observation, so we measure the systemic velocity of our templates and add it to the velocity difference between an observation and the corresponding template. To measure the velocity of templates, we find the lines using the method described above and obtain a CCF by cross correlating this line list with a PHOENIX synthetic spectrum matching the star's temperature; these spectra are drawn from the Goettingen spectral library (Husser et al. 2013). The peak of the CCF using an empirical line position onto a synthetic spectrum is expected to be opposite of the star's systemic velocity. By using different templates on the same time series, we can estimate the typical errors on the systemic velocity. We used the Barnard's star sequence described here and measured its velocity with templates from 11 stars within 100 K (GJ1012, GJ1103, GJ1105, GJ1151, GJ1289, GL15B, GL169.1A, GL445, GL447, Barnard's star, and GL725B). The systemic velocity of Barnard's with its own template is −109,974 m s⁻¹, while the dispersion of systemic velocities derived using all templates is 37 m s⁻¹. We therefore estimate that the typical uncertainties on systemic velocities from the LBL are at that level for high-S/N observations. We note that some systemic effects, such as uncertainties in the line positions of models or limits to our knowledge of the gravitational redshift, are common to all stars and not accounted for in this error estimate.

## 4.10. High-pass Filtering Data

One problem in subtracting a template from a science observation is the proper scaling of the template in amplitude. One can leave each line with an arbitrary scaling factor, as is done in Dumusque (2018), or adopt a coherent amplitude over a larger spectral domain. Furthermore, IR astronomical arrays display persistence; this is an effect wherein a ghost image of the previous observation is detected at reduced amplitude in subsequent images, even after detector reset (e.g., Artigau et al. 2018a; Baril & Albert 2008). In NIR high-resolution spectroscopy, persistence mostly displays as a dilution of the line depths as it adds a continuum contribution that is distinct from that of the target. Modeling persistence is challenging as it requires knowledge of the entire time series of the night, including calibrations. In order to avoid problems linked to the dilution of the continuum in NIR observations both the science spectrum and the template are high-passed by subtracting a low-pass filter over a 500 km s⁻¹ running window. The high-pass filtering is performed by subtracting the running mean of a window of a given width centered on the wavelength of interest. We explored windows ranging down to 50 km s⁻¹, and
this parameter has minimal impact on the end result. This effectively removes the persistence contribution to the continuum level, an issue that has been shown to be problematic in the case of observing fainter M dwarfs with SPIRou, in particular TRAPPIST-1, immediately after either telluric standards or transit spectroscopy targets that compound a bright target with a long staring, worsening persistence effects. After high-passing, the template is multiplied by a scaling factor to match the amplitude of the science spectrum through a least-squares adjustment. This is done per-order, allowing for slight large-scale chromatic losses between science observations. Chromatic losses could arise from varying scattering properties of the atmosphere, coating degradation in telescope optics or varying fiber coupling efficiency under varying seeing, the seeing profile being mildly chromatic (Fried 1966). In typical SPIRou observations, the median J-to-K flux ratio typically varies by up to 15% on monthly timescales.

### 4.11. Activity Indicator

The B01 framework corresponds to a projection of the residuals between a given spectrum and the corresponding template onto the first derivative of the template with wavelength. One can extend this framework, by letting the spectrum A(i) line profile width vary by an infinitesimal amount δw, relative to the mean width w of the template A₀(i):

$$A(i) - A_0(i) = \frac{\partial A_0(i)}{\partial \lambda(i)} \delta\lambda(i) + \frac{\partial A_0(i)}{\partial w(i)} \delta w(i). \quad (6)$$

In this context, a spectrum as being the sum of well-separated Gaussians having the same widths leads to an interesting analytical solution. Admittedly M-dwarf spectra are much more complex with numerous overlapping features, but this framework does provide insight into the physical interpretation of the projection of residuals onto higher-order derivatives. For a Gaussian line having a width w_g centered at λ_0, the equivalence between the first derivative against w_g and the second derivative with respect to wavelength is given by

$$\frac{\partial}{\partial w_g} \left( \frac{e^{-0.5 \left( \frac{\lambda - \lambda_0}{w_g} \right)^2}}{w_g} \right) = w_g \frac{\partial^2}{\partial \lambda^2} \left( \frac{e^{-0.5 \left( \frac{\lambda - \lambda_0}{w_g} \right)^2}}{w_g} \right). \quad (7)$$

In other words, changes in the line profile width (δw) from spectrum to spectrum may be approximated as a second derivative term in wavelength, as long as δw/w ≪ 1. In the case of rich molecular bands, this is an approximation as a line width cannot be rigorously defined, but it would be strictly true for isolated Gaussian lines with a w_g width. By taking Equations (6) and (7), we have

$$A(i) - A_0(i) = \frac{\partial A_0(i)}{\partial \lambda(i)} \delta\lambda(i) + \frac{\partial^2 A_0(i)}{\partial \lambda(i)^2} w_g(i) \delta w_g(i). \quad (8)$$

This correspondence between the second derivative of the spectrum and the width for Gaussian lines was previously discussed in detail in Zechmeister et al. (2018; Section 3.2 and Appendix B therein); they propose an activity indicator, the differential line width (dLW ≡ σΔσ), where σ corresponds to the equivalent width of Gaussian lines expressed in velocity unit, and Δσ, the variation of this width. Taking Equation (8) for processes with $\frac{\partial A_0(i)}{\partial \lambda(i)} = 0$ (e.g., symmetric line profile

changes such as a broadening), we have

$$\frac{\mathrm{dLW}(i)}{c^2} = \frac{A(i) - A_0(i)}{\lambda(i)^2 \frac{\partial^2 A_0(i)}{\partial \lambda(i)^2}}. \quad (9)$$

Following the same formalism as for Equations (3) and (4), one can derive a LBL dLW, which can then be averaged out in an outlier-resilient way (as described in Appendix B). The dLW indicator has m² s⁻² units and scales with the variation in the line profile FWHM, a metric widely used in the literature to characterize the activity level of a star. Considering that a Gaussian line has a FWHM = $2\sqrt{2 \ln 2} w_g \approx 2.35w_g$, one obtains

$$\mathrm{dLW} = \frac{\mathrm{FWHM} \times \Delta\mathrm{FWHM}}{8 \ln 2}. \quad (10)$$

However, one caveat in converting a dLW measurement into ΔFWHM is that knowledge of the effective FWHM of the template is required. One can derive this value by using a CCF to allow a direct comparison with FWHM changes detected through other methods. If dLW is used as an activity indicator alone, one can use it directly for period searches and detrending. This is demonstrated on the Proxima Centauri data in Section 5.3. The weighted mean in velocity and dLW are computed independently for the sake of simplicity; this implies that a given point may be an outlier within one distribution but not in the other. We tested a joint rejection, and the impact is minimal, typically well below 0.1σ for either value.

One notable characteristic of the dLW is that a variation in the instrument line-spread-function profile has the same impact regardless of the stellar FWHM. This is not the case for changes in the FWHM as measured with a CCF; the larger the FWHM, the smaller the impact of a change in the instrumental profile. This implies that long-term monitoring of activity can be corrected for instrumental profile changes directly by using the dLW, as measured on a set of calibrator stars or with calibration data such as the FP spectrum used to determine the wavelength solution.

## 5. Demonstrations of the Line-by-line Method

The LBL algorithm allows the measurement of Keplerian signals at a few meters per second level with meters-per-second-level uncertainties in SPIRou data sets, as presented in recent work. It has been used to measure TOI-1759b's mass through a measurement of a pRV signal of K = 2.3 ± 0.7 m s⁻¹ (Martioli et al. 2022). Similarly, it has been used to measure the mass of the sub-Neptunes orbiting TOI-2136 and TOI-1452 (Gan et al. 2022; Cadieux et al., submitted), respectively, corresponding to signals at the K = 4.2 ± 1.4 m s⁻¹ and 3.5 ± 0.9 m s⁻¹ level. These discoveries are based on individual measurements at the 6 to 10 m s⁻¹ accuracy level with typically four measurements (one polarimetric sequence) per epoch, which are limited by the modest S/N obtained on these relatively faint targets (H = 9.6–10.0) and do not represent a fundamental noise floor for the technique. Here, we present a detailed analysis of the multiyear time series obtained on Barnard's star at IR wavelength and Proxima in the optical (respectively with SPIRou and HARPS), focusing on the statistical properties of these data sets in comparison with other pRV measurement methods.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Radial velocity (RV) time series of Barnard's star with the line-by-line (LBL) method and the standard SPIRou cross-correlation function. No activity correction has been performed. The median LBL RV uncertainty is 2.57 m s⁻¹, only slightly smaller than the point-to-point dispersion of 2.59 m s⁻¹. The Ribas et al. (2018) circular orbital solution is shown. The data strongly favors the scenario in which the RV is constant compared to that with the candidate planet with the reported solution (see Section 5.1).</figcaption></figure>

## 5.1. Line-by-line versus Cross-correlation Function on Barnard's Star

Barnard's star has been observed regularly since the first light of SPIRou, and it serves as a benchmark object for the development of SPIRou data-analysis tools. The star is known to display little activity due to its slow rotation and thick-disk age. It is also among the brightest M dwarfs in the sky due to its proximity (1.83 pc; Gaia Collaboration et al. 2021). Furthermore, the recent claim of a super-Earth on a 233 day orbit (Ribas et al. 2018) and the subsequent attribution of the pRV signal to activity (Lubin et al. 2021) make it a compelling target. We defer the full analysis of the Barnard's star RV time series to a future publication (T. Vandal et al. 2022, in preparation), and here present a summary.

We computed Barnard's star's complete RV time series with the LBL and the CCF methods. The CCF RVs were extracted using the *spirou-ccf*<fnref n="11" /> package (see Martioli et al. 2022 for details). The resulting RV curves are presented in Figure 5, with the circular Keplerian solution from Ribas et al. (2018) overplotted. As shown in the velocity histogram in Figure 6, the LBL algorithm provides twice less point-to-point scatter (rms = 2.59 m s⁻¹), compared to the CCF-derived velocities (rms = 7.53 m s⁻¹). The LBL algorithm returns RV uncertainty estimates (2.57 m s⁻¹) that are consistent with the noise level as determined by subtracting a template from the data, but this noise level is, itself, significantly higher than what would be expected from Poisson statistics in the *H* and *K* bands. The underlying reason for this excess noise is still being investigated, but modal noise (Blind et al. 2017; G. Gaisné et al. 2022, in preparation) is one of the main plausible causes, which would explain the worsening at longer wavelengths

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Histogram of radial velocities relative to the median velocity for Barnard's star observed with SPIRou. Line-by-line values have a standard deviation of 2.59 m s⁻¹ for a 2.57 m s⁻¹ median uncertainty. Cross-correlation function values have a 7.53 m s⁻¹ standard deviation.</figcaption></figure>

q(i.e., a fiber has fewer modes at longer wavelengths and a fractionally larger modal noise). With only Poisson noise, one would expect an RV accuracy at the 1.0 m s⁻¹ level.

Given the 2.59 m s⁻¹ dispersion in our Barnard's star LBL time series spanning over 2.7 yr and typical uncertainty of 2.57 m s⁻¹, one would expect the candidate planet unveiled by Ribas et al. (2018; *K* = 1.20 ± 0.12 m s⁻¹) to dominate the excess noise budget. We conducted a statistical analysis to objectively determine whether this signal is detected in our LBL data. First, we computed the Bayesian log-evidence (ln*Z*) using the *nestle* (Barbary 2021) nested sampling algorithm for a flat model, i.e., systemic velocity and extra white noise terms only, and for either of Ribas et al.'s (2018) circular and eccentric solutions. Both planetary models are strongly rejected by our data with a difference in log-evidence of Δln*Z* = -23.9 for the eccentric and Δln*Z* = -33.0 for the circular solution compared to the flat solution. This result implies that other sources of noise (e.g., instrumental, stellar jitter, planets) are needed to explain the observed dispersion. These will be addressed in T. Vandal et al. (2022, in preparation), with no evidence of the claimed super-Earth in our data set. Then, we explored the possibility of a 233 day planet, but with a different semiamplitude *K* than the one reported in Ribas et al. (2018). For this, we generated Keplerian models with *radvel* (Fulton et al. 2018) and sampled the posterior distribution of *K* with a Markov Chain Monte Carlo (MCMC) Bayesian framework using *emcee* (Foreman-Mackey et al. 2013). We employed a uniform prior distribution *U*(0, 5) for *K*, while the orbital parameters (period and time of periastron) followed a Gaussian prior with standard deviation corresponding to their uncertainty reported in Ribas et al. (2018). The plausible values for *K* (posterior distribution) were incompatible with *K* = 1.20 ± 0.12 m s⁻¹ at the >5σ significance level. Our analysis therefore strongly disfavors the existence of a 233 day planet around Barnard's star with a minimal mass similar to the candidate super-Earth revealed by Ribas et al. (2018).

## 5.2. Line-by-line versus Template Matching on Proxima Centauri

While the LBL algorithm has been developed with IR measurements in mind, we tested it on publicly available HARPS data on Proxima Centauri. In addition to being our nearest stellar neighbor, the star is known to host an Earth-mass

<footnote n="11">¹¹ github.com/edermartioli/spirou-ccf</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Radial velocity measurements of the HARPS post-fiber change on Proxima Centauri with the line-by-line (LBL) and HARPS-TERRA algorithms. The top panel shows both data sets with Proxima b's and d's Keplerian contribution (solution from Faria et al. 2022). The bottom panel shows the difference between LBL and HARPS-TERRA; within each sequence the agreement is excellent (rms difference of 0.7 and 1.2 m s⁻¹ in 2016 and 2017, respectively), but a systematic difference between sequences at the ~3 m s⁻¹ level is seen. Uncertainties in the difference between LBL and HARPS-TERRA are derived from the per-sequence rms. One cannot simply add the errors from each method in quadrature as the noise sources in the two are highly correlated; the two used the same input data sets with the same noise, albeit with different weighting schemes.</figcaption></figure>

companion close to its habitable zone (Anglada-Escudé et al. 2016). More recently, a second close-in terrestrial planet was identified by Faria et al. (2022). This well-studied star thus provides an interesting benchmark for pRV measurement techniques.

The LBL RVs of Proxima were computed using a template generated from all HARPS spectra collected after the fiber change. We compare below the LBL analysis of Proxima with the Damasso et al. (2020) publicly available RVs derived with the HARPS-TERRA template-matching algorithm. We also computed the velocity of the time series with NAIRA (see Astudillo-Defru et al. 2017, Section 3.1), an algorithm to compute RVs through the template-matching approach.

The NAIRA algorithm Doppler shifts all spectra to the stellar frame to compute a median spectrum (stellar template), which is then used to obtain a telluric template from the median of spectra divided by the stellar template shifted to the Earth's frame. We iterate the procedure to properly disentangle the stellar and telluric components, which in turn depends on the number of spectra acquired in different BERVs. The telluric spectrum is used to mask only the zones affected by tellurics, instead of discarding ±30 km s⁻¹ around a given telluric line. The stellar template is shifted in a range of velocities around the star systemic velocity to compute the likelihood for each RV step. The likelihood is maximized to obtain the optimal RV of the star.

The main difference between HARPS-TERRA and NAIRA comes from the weights applied to each pixel and individual RVs from each aperture. NAIRA performs a simple chi-square minimization assuming Poisson statistics, HARPS-TERRA uses weights that are proportional to the residuals; NAIRA computes the minimization on all apertures (orders) at the same time, while HARPS-TERRA obtains an order-by-order minimization to apply weights from their dispersion. Furthermore, HARPS-TERRA masks tellurics deeper than 1% in a synthetic spectrum, and NAIRA masks tellurics for which the Doppler information is significantly larger than the Doppler content of a stellar line located at the same wavelength.

This analysis on a common data set allows a direct comparison of the LBL method with state-of-the-art template-matching algorithms, and an understanding of eventual

differences. The differences seen in the velocity put the detection of the candidate planet at long period identified by Damasso et al. (2020) in a new light.

The analyzed Proxima Centauri data consist of two densely sampled sequences, one in 2006 late January to early April and the other in 2017 July to September. A third sequence, obtained in 2019 March to July, was analyzed with the LBL and NAIRA algorithms but has not been analyzed by Damasso et al. (2020).

Figure 7 presents the main differences between LBL and HARPS-TERRA RVs during the common epochs of 2016 and 2017. If the two extraction methods produce almost perfectly matching RV in a timescale of months, the HARPS-TERRA deviate by ~3 m s⁻¹ between 2016 and 2017. This trend is absent from the LBL time series. Similarly, Figure 8 provides a comparison between the LBL and NAIRA methods, showing an agreement between the two measurements at the ~70 cm s⁻¹ level without any significant difference between sequences.

The most notable difference between the two campaigns is the BERV change: the 2016 sequence corresponds to the largest yearly line-of-sight value (from 16 to 21 km s⁻¹), while the 2017 sequence corresponds with the lowest BERV excursion (−14 to −21 km s⁻¹ over these dates). Any RV extraction method free of systematic effects should yield residuals that are uncorrelated with the BERV once the Keplerian and activity signals are removed. This can easily be tested on the Proxima data using four RV extraction methods: LBL, CCF, HARPS-TERRA, and NAIRA. To quantify this systematic effect, we define the quantity $K_{BERV}$ as the semiamplitude of the yearly peak-to-peak signal induced by the correlation between RV residuals and the BERV for that specific target and time sequence. We first removed Proxima b's and d's Keplerian signals in all four data sets using the Faria et al. (2022) orbital solution. Modeling the activity signal for the four extraction methods was out of the scope of this simple analysis. Instead, we estimated how much the activity can contribute in $K_{BERV}$ due to random and finite sampling of Proxima's rotation and BERV. We measured the slope of 1000 random realizations of the activity model of the LBL RV time series presented in Section 5.4 with BERV. From this, we derived $K_{BERV}$ values strictly due to activity with a
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Same as Figure 7, but for the NAIRA time series. The measurements are in remarkable agreement, with a 0.7 m s⁻¹ rms difference between values. The ~3 m s⁻¹ bulk offset between the 2016 and 2017 data sets, as seen in the line-by-line to HARPS-TERRA comparison, is not present. Overall, the scatter after subtraction of the Proxima b and d ephemerides is the same at the 1.9 m s⁻¹ level for both time series.</figcaption></figure>

mean of 0.04 m s⁻¹ and standard deviation of 0.46 m s⁻¹. We thus added this uncertainty term in quadrature to all four $K_{BERV}$'s presented below. As shown in Figure 9, both the CCF and HARPS-TERRA RV residuals show a trend with BERV, i.e., $K_{BERV} = 0.9 \pm 0.5$ and $1.8 \pm 0.5$ m s⁻¹, respectively, but this slope is consistent with zero for LBL and NAIRA ($K_{BERV} = 0.3 \pm 0.5$ and $0.4 \pm 0.5$ m s⁻¹). Figure 9 strongly suggests that the LBL and the NAIRA template-matching methods are largely immune to the systematic effects affecting HARPS-TERRA and, to a lesser extent, the CCF measurements. Determining the underlying causes leading to the correlation of velocities with BERV in CCF and HARPS-TERRA data is beyond the scope of this paper, but we note that this behavior is expected for suboptimal telluric-line corrections. While developed in the context of NIR observations, this analysis illustrates the benefit of the outlier statistical rejection that the LBL allows at optical wavelengths. The differences between these models cannot be attributed to differences in zero-point or wavelength calibration as they all use the same wavelength grid provided by the HARPS pipeline.

While pRV time series analysis requires a proper accounting of stellar activity for accurate determination of Keplerian orbits (see Section 5.4), the RV measurements prior to any activity correction are worth discussing. Figure 10 shows the three sequences analyzed with the LBL and the two with HARPS-TERRA compared with the Keplerian motion expected from the candidate planet, Proxima c, proposed by Damasso et al. (2020). One sees that the LBL, and thus NAIRA (see Figure 8), data are incompatible with the velocity slope expected for the Keplerian motion of Proxima c. The HARPS-TERRA values do display a velocity change between 2016 and 2017 that could be interpreted as a genuine Keplerian signal, but it was shown above that this variation can be explained by a systematic effect. A more quantitative analysis is presented later that further questions the existence of Proxima c.

## 5.3. Line-by-line Activity Indicator On Proxima

Activity indicators are central in the interpretation of a pRV sequence as they sample the physical processes that create activity jitter, which can mimic Keplerian signals. The activity indicator variations are modulated at the stellar rotation period, and therefore consist of independent observations useful for detrending the RV curve from activity and extract meters-per-

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Radial velocity residuals as a function of barycentric radial velocity after subtracting Proxima b and d signals using the Faria et al. (2022) orbital fit. The HARPS-TERRA values display a $K_{BERV}$ signal at the 1.8 m s⁻¹ level while it is consistent with zero for the line-by-line and NAIRA analysis. The cross-correlation data show an intermediate $K_{BERV}$.</figcaption></figure>

second-level planetary signals. One way to do this is to first model the activity tracer (e.g., photometry, FWHM, bisector tilt, dLW) with a Gaussian process (GP), and then use this “trained” GP to correct the RVs or to jointly fit the RV and activity indicator. While the rotation period is often referred to with a single value, it is important to keep in mind that one expects this signal to have a time-dependent periodicity. In the case of Proxima, one expects that differential rotation and secular differences on the inferred period may simply reflect latitudinal changes in the distribution of active regions with time. Suárez Mascareño et al. (2016) and Wargelin et al. (2017) point to a ~7 yr activity cycle, implying that the latitudinal distribution of active regions is likely to evolve on that timescale. Furthermore, different methods have differing latitude-weighted sensitivity. FWMH and bisector measurements are fractionally more affected by equatorial spots than photometric measurements; one therefore does not expect a perfect correspondence between measurements taken at varying epochs and methods. We here present a summary of the key literature values for the rotation period of Proxima:

1. Benedict et al. (1998) used HST guider photometry to measure a 83.5 day rotation period (no uncertainty reported) between 1992 and 1997 in a bandpass centered at 583 nm.

2. Wargelin et al. (2017) in a 15 yr long compilation of All Sky Automated Survey (Pojmanski 1997) V-band photometry (2001–2015) found a period ranging from
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Residual time series for Proxima Centauri with the line-by-line (LBL) and HARPS-TERRA data sets after subtracting the Proxima b and d Keplerian solutions of Faria et al. (2022). The HARPS-TERRA time series shows a decrease of 3.2 m s⁻¹ between the two runs, consistent with the reflex motion due to Proxima c. The LBL time series, extending over larger set of observations in mid-2019, shows no such trend and is inconsistent with the Damasso et al. (2020) orbital solution for Proxima c.</figcaption></figure>

77.1 days to 90.1 days, pointing to an evolving lightcurve and providing a strong case for differential rotation. It is noteworthy that this analysis gives the largest span of measured rotation periods and is also the longest one, covering more than two activity cycles. This compilation supersedes Kiraga & Stepien (2007) and Suárez Mascareño et al. (2016), which used earlier releases of the same time series; the rotation periods determined by these authors are consistent with Wargelin et al.'s (2017) analysis (see Section 2 therein).

3. Newton et al. (2018), using photometry obtained from 2014 to 2018 with MEarth (Irwin et al. 2009), measured a $P_{\text{rot}} = 89$ days with a typical 10% uncertainty for their whole sample.

4. Klein et al. (2020) identified a period of 89.8 ± 4 days through optical spectropolarimetry in a sequence contemporaneous with our 2017 data set.

Using the LBL framework, we can infer strong constraints on the Proxima rotation period over the timespan of the data set by using the dLW activity indicator. This sequence is compared with the CCF FWHM from the HARPS pipeline. As for the first-order derivative measurement at the core of the LBL method, the dLW has the notable advantage over the CCF FWHM of being computed in an outlier-resilient way.

We modeled the two time series with a GP using the quasiperiodic covariance kernel, a widely used kernel for modeling variability generated by stellar rotation (Angus et al. 2018). The details of our GP prescription are given in Appendix C and Table C1 provides the complete list of fitted hyperparameters.

Figure 11 summarizes the comparison between dLW and CCF FWHM for modeling Proxima Centauri activity. The two indicators are highly correlated, particularly during the 2016 and 2019 sequences. The quasiperiodic GP models yield a rotation period of $P_{\text{rot}} = 92.1_{-3.5}^{+4.2}$ days for the LBL dLW and

$P_{\text{rot}} = 79.6_{-4.2}^{+9.7}$ days for the CCF FWHM. The posterior distribution obtained with dLW is more symmetric and has a smaller $1\sigma$ confidence interval than the CCF-derived one. While both indicators yield periods consistent with literature values, we note that the LBL dLW provides a rotation period closer to the contemporaneous measurements made by Newton et al. (2018) and Klein et al. (2020).

## 5.4. Keplerian and Activity Model of Proxima

In this section, we present a joint analysis (Keplerian and activity) of the complete Proxima LBL RV time series (2016, 2017, and 2019 sequences). We adopt the naming scheme wherein Proxima b is the 11.18 day planet (Anglada-Escudé et al. 2016; Suárez Mascareño et al. 2020; Faria et al. 2022), Proxima c is the 1900 day candidate planet reported by Damasso et al. (2020), and Proxima d is the 5.15 day signal first identified as a possible companion with a $0.29 \pm 0.08 M_{\oplus}$ minimum mass in Suárez Mascareño et al. (2020), later confirmed in Faria et al. (2022). We follow this convention, acknowledging that the naming may evolve. We seek to quantify the level of each of the proposed planetary signals in RVs produced by the LBL algorithm.

We investigated the following suite of models ($M$):

1. $M_{\text{dLW}}$: No planet, activity GP only.

2. $M_{\text{b+dLW}}$: One planet, Proxima b, with activity GP.

3. $M_{\text{bd+dLW}}$: Two planets, Proxima b, d, with activity GP.

4. $M_{\text{bc+dLW}}$: Two planets, Proxima b, c, with activity GP.

5. $M_{\text{bcd+dLW}}$: Three planets, Proxima b, c, d, with activity GP.

For the GP activity component, we utilized a quasiperiodic kernel trained on the dLW time series, that is, with prior distributions on shared hyperparameters corresponding to the posteriors from the GP regression presented in Section 5.3. For the Keplerian component of the joint fit, we allowed nonzero values for the eccentricity of Proxima b, while fixing $e = 0$ for Proxima c and Proxima d when included in the model. We used noninformative priors for the orbital period ($P$) and the time of inferior conjunction ($T_0$) of Proxima b (see Table C2). For the models with Proxima c, we used Gaussian priors on $P_c$ and $T_{0,c}$ based on Damasso et al. (2020) as our data do not cover a full period of the candidate planet. The priors on the period and phase of Proxima d were uniform, but fairly narrow around the value from Faria et al. (2022; see Table C3). This is because our data set does not include ESPRESSO data, which was mainly responsible for the detection. This model should thus be viewed as an informed search for planet d in our data, not a blind search, similarly to a transiting planet with a known period and phase. Appendix C provides additional details on the joint fits and Table C2 summarizes prior and posterior distributions.

We computed the Bayes Information Criterion (BIC; Schwarz 1978) for each $M$ considered. The smallest BIC, corresponding to the preferred model, was obtained for the single planet $M_{\text{b+dLW}}$, with a BIC $b+dLW = 739.1$. Compared to the no planet, activity GP only solution, we measure a $\Delta$BIC $\equiv$ BIC $b+dLW$ $-$ BIC $dLW = -37.6$, providing compelling evidence of the Proxima b Keplerian signal in our data. The two planets model $M_{\text{bd+dLW}}$ produced a $\Delta$BIC $= -5.8$, which advocates against the necessity of adding the $\sim 0.4$ m s⁻¹ Proxima d signal to explain our data. Taken without further input, our data and analysis would therefore be insufficient to
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Proxima Centauri activity modeling with the dLW and the CCF FWHM. Left panels: changes in the dLW line shape parameter compared to the CCF FWHM. The green and red curves represent the mean prediction of the best-fit quasiperiodic Gaussian process (GP; described in Appendix C), each surrounded by a 68% confidence envelope. Middle panel: direct comparison between the dLW (expressed in FWHM) and the CCF FWHM. One does not expect a perfect correspondence between the two variables, as the equivalence between the dLW and the FWHM is only strictly true for Gaussian lines. Right panel: posterior distributions on $P_{\text{rot}}$ from the quasiperiodic GP model of the dLW and CCF FWHM time series. The 68% confidence intervals on $P_{\text{rot}}$ are represented with green (LBL) and red (CCF) regions. The black dashed lines show the 10% relative uncertainty interval for Proxima's rotation period ($P_{\text{rot}} = 89$ days) from Newton et al. (2018). There is no specific uncertainty for the Newton et al. (2018) rotation period, but the typical uncertainties in their sample are at the 10% level.</figcaption></figure>

identify this signal. This is unsurprising as our set of observations consists of a subsample of Suárez Mascareño et al. (2020) and does not include ESPRESSO RVs. We nonetheless point out that the LBL method yields a $P_d$ and $K_d$ (see Table C3 and Figure C1) fully consistent with the results of Suárez Mascareño et al. (2020) and Faria et al. (2022), with $K_d$ detected at the $2\sigma$ level.

Other models including Proxima c are strongly rejected according to their BIC ($\Delta$BIC = -15.1 for $M_{bc+dLW}$, $\Delta$BIC = -21.2 for $M_{bcd+dLW}$), providing quantitative evidence against the existence of the candidate planet signal in the LBL data. This conclusion is further strengthened by the fact that the NAIRA analysis of the same data set provides nearly identical results to the LBL and shows no significant RV difference between the 2016 and 2017 sequences. We note that Damasso et al. (2020) did not clearly identify this signal in their analysis of the HARPS data only (not considering UVES observations), but the inclusion of the 2016 and 2017 sequences had a large impact on the final period and increased by 25% the semiamplitude (Figure 2, therein). Combined with the absence of a clear detection of Proxima c in imaging and in astrometry (Gratton et al. 2020; Kervella et al. 2020; Benedict & McArthur 2020), our results cast serious doubts on the existence of this planet.

## 6. Conclusions

The results described here illustrate the usefulness of the LBL algorithm in handling systematic uncertainties in pRV. It outperforms, in terms of raw RV uncertainties, the CCF method as it fully utilizes the RV content of the input spectra. At IR wavelengths, the LBL algorithm opens the door to velocimetry at the meters per second level, as exemplified by the constraints added to the candidate planet around Barnard's star presented here and mass measurements with SPIRou of planets uncovered through Transiting Exoplanet Survey Satellite photometry.

At optical wavelengths, the method appears on par with the template-matching method in terms of accuracy, which is not overly surprising as both methods use the full RV content in spectra. However, LBL appears to display lower levels of

systematic drifts over timescales of many months as it can efficiently reject the part of the spectra affected by telluric residuals. Benchmarking over tens of targets would be necessary to disentangle the exact causes of differences in BERV trends.

A number of byproducts of the LBL algorithm could be used to further constrain activity signals. The velocity spectrum, as shown in Figure 3, provides the opportunity to measure the chromatic velocity signal in a self-consistent manner with self-consistent uncertainties. This has proved to be a useful tool in disentangling activity signals for very active objects (e.g., Zechmeister et al. 2018; Cale et al. 2021). Sorting lines by their strength may allow for the probing of different depths in the stellar photosphere and could be affected differentially by activity. This avenue has not yet been explored systematically.

As detailed in Section 4.11 and in Zechmeister et al. (2018), the projection on the second derivative is a very useful tool for activity monitoring, providing a more robust tool than the CCF FWHM for the determination of the rotation period in the Proxima Centauri data (see Figure 11). Projections onto higher-order derivatives could be interesting as they would provide information on changes in higher moments of the line profiles. This would be conceptually analogous to the approach of Collier Cameron et al. (2021), wherein activity is measured by a decomposition of the spectrum autocorrelation function into principal components. Numerically, we have found that the projection of the third derivative of the line profile is directly linked to the bisector tilt that traces changes in line symmetry, though we have been unable to find an analytical correspondence similar to the one presented in Section 4.11 linking the change in the second-order derivative with changes in the FWHM. Assessing the performance of the third-order derivative as an activity indicator, in comparison with indicators such as bisector slope changes, is the subject of ongoing work. Similarly to chromatic activity, a future avenue to account for stellar activity is to measure velocity with a selected subset of lines, such as is done in Bellotti et al. (2022). This has been done in the case of Sun-like stars (Dumusque 2018), but has yet to be explored for M dwarfs and in the NIR.

Our analysis of the HARPS Proxima Centauri data is a warning that long-period, low-mass planets characterized by
small semiamplitudes at the level of a few meters per second may be resulting from unaccounted for systematic effects and warrants a reanalysis of existing data using the LBL framework.

This work is based on observations obtained at the Canada-France-Hawaii Telescope, which is operated from the summit of Maunakea by the National Research Council of Canada, the Institut National des Sciences de l'Univers 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. This work is party supported by the Natural Science and Engineering Research Council of Canada, the Fonds Québécois de Recherche (Nature et Technologie), and the Trottier Family Foundation through the Institute for Research on Exoplanets.

The authors thank Rodrigo F. Díaz for suggestions regarding the use of the mixture model formalism in the LBL framework. E.A. thanks Thomas Davy for the fruitful discussions regarding the underlying logic of the pRV measurements with SPIRou that ultimately led to the current work. J.-F.D. acknowledges funding from the European Research Council under the H2020 Research & Innovation program (grant #740651 NewWorlds). X.D. and A.C. acknowledge funding from ANR of France under contract number ANR18CE310019 (SPlaSH). This work is supported by the French National Research Agency in the framework of the Investissements d'Avenir program (ANR-15-IDEX-02), through the funding of the "Origin of Life" project of the Université Grenoble Alpes. N. A.-D. acknowledges the support of FONDECYT project 3180063.

**Softwares:** This research made use of the following software tools: 1. Astropy, a community-developed core Python package for Astronomy (The Astropy Collaboration et al. 2013, 2018). 2. bottleneck (bottleneck.readthedocs.io). 3. ds9, a tool for data visualization supported by the Chandra X-ray Science Center (CXC) and the High Energy Astro-physics Science Archive Center (HEASARC) with support from the JWST Mission office at the Space Telescope Science Institute; (Joye & Mandel 2003). 4. emcee (Foreman-Mackey et al. 2013). 5. george (Ambikasaran et al. 2016). 6. IPython (Perez & Granger 2007). 7. matplotlib, a Python library for publication quality graphics (Hunter 2007). 8. Numba (Lam et al. 2015). 9. NumPy (Harris et al. 2020). 10. radvel (Fulton et al. 2018). 11. SciPy (SciPy 1.0 Contributors et al. 2020). 12. tqdm (Da Costa-Luis et al. 2021).

**Facilities:** CFHT(SPIRou), ESO La Silla (HARPS), ESO Science Archive.

## Appendix A
Code Accessibility

The most up-to-date implementation of the LBL algorithm is available to the community at tbl.exoplanets.ca, with a link to the relevant GitHub repository and sample data sets. Tools to generate the templates and masks are also provided. The rdb format output of the LBL codes is compatible with the Data &

Analysis Center for Exoplanets<fnref n="12" /> online tools to analyze pRV time series.

## Appendix B
Weighted Mean in the Presence of Outliers

The LBL method requires the combination of a large number of velocities for which uncertainty estimates are known, but allowing for the possibility that a small fraction of values are drawn from a population of spurious values for which we have no prior knowledge of the expected distribution. The fraction of lines expected to belong to the spurious population can be estimated from the fraction of high-$\sigma$ outliers in the velocity distribution. Assuming that the outlier represents a fraction $f_0$ of all lines, that each $i$ point has an uncertainty estimate $\sigma_i$, a value $v_i$, a weighted mean velocity $\bar{v}$, the normalized distance to mean $\sigma_\mu = \frac{v_i - \bar{v}}{\sigma_i}$, the likelihood that one point is a valid point is

$$p_i = \frac{(1 - f_0)e^{-\sigma_\mu^2/2}}{(1 - f_0)e^{-\sigma_\mu^2/2} + f_0}.$$

The mean velocity is assumed to be the same for all lines and can be determined through a weighted sum that takes into account both the likelihood of the point being valid and the uncertainties as derived from the LBL algorithm:

$$\bar{v} = \frac{\sum_i v_i p_i \sigma_i^{-2}}{\sum_i p_i \sigma_i^{-2}}.$$

The corresponding uncertainty for the mean velocity is

$$\bar{\sigma} = \sqrt{\frac{1}{\sum_i p_i \sigma_i^{-2}}}.$$

The above equations are coupled but converge rapidly with an iterating algorithm. One starts with an approximation of $\bar{v}$ by using the median value of $v_i$ to determine $p_i$, and uses these value to refine the $\bar{v}$ value. Assuming that $f_0$ is very small (typically <0.1% in SPIRou data), $\bar{v}$ converges to a small fraction of $\bar{\sigma}$ within two iterations. We compiled the number of outliers in the Barnard's star data, consisting of $4.6 \times 10^6$ individual line measurements (247 spectra and ~19,000 lines per spectrum) and found that $5.8 \times 10^{-4}$ of the lines are >5$\sigma$ outliers. This is a much larger fraction than would be expected from a normal distribution (~6 $\times 10^{-7}$) and suggests that these lines belong in vast majority to the "outlier" population.

To determine the experimental $f_0$ value, we performed a Gaussian fit on the density distribution of $\frac{v_i - \bar{v}}{\sigma}$ (see Figure B1), subtracted the fit from the distribution and found an excess of $2 \times 10^{-3}$ of values beyond $\pm 3\sigma$ compared to a normal distribution. We therefore adopt $f_0 = 2 \times 10^{-3}$ for all our analyses, though the resulting pRV time series are largely insensitive to changes in $f_0$ ranging from $10^{-2}$ to $10^{-4}$. An $f_0 = 2 \times 10^{-3}$ will lead to a $p_i = 0.5$ for a $3.5\sigma$ event and $p_i \sim 0.02$ for a $4.5\sigma$ event. In practice, this summation scheme is equivalent to a ~4$\sigma$ clipping with soft edges rather than binary rejection. This weighted mean is used both for the determination of the mean velocity and dLW with the same $f_0$ value.

<footnote n="12">¹² dace.unige.ch/dashboard/</footnote>
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure B1.** Left: for a sample Barnard’s star observation, the histogram of the per-line velocity difference from the mean velocity divided by the corresponding line uncertainty closely follows a normal distribution (dashed line) with only a handful of >5σ outliers. Right: spectral distribution of >5σ outliers in a random Barnard’s star observation.</figcaption></figure>

## Appendix C
Material Regarding Gaussian Process Models

In Section 5.3, we modeled Proxima stellar activity in the LBL dLW and the HARPS pipeline CCF FWHM time series with a GP. We employed the following quasiperiodic covariance kernel:

$$k_{i,j} = A^2 \exp \left[ -\frac{|t_i - t_j|^2}{2\ell^2} - \Gamma^2 \sin^2 \left( \frac{\pi|t_i - t_j|}{P_{\text{rot}}} \right) \right] + s^2 \delta_{i,j}, \quad (C1)$$

where $|t_i - t_j|$ is the time interval between data i and j, A is the amplitude of the GP, $\ell$ is the correlation timescale, $\Gamma$ scales the periodic component of the GP, $P_{\text{rot}}$ is the stellar rotation period, and s is an extra white noise term added in the diagonal of $k_{i,j}$. Our GP modelization consisted of five hyperparameters, all fitted in logarithmic form {$\ln A, \ln \ell, \ln \Gamma, \ln P_{\text{rot}}, \ln s$}. We used george (Ambikasaran et al. 2016) to compute GP predictions and calculate the log-likelihood. The parameter posterior distributions were sampled using a MCMC algorithm with emcee (Foreman-Mackey et al. 2013). We used 100 walkers, ran 25,000 steps, and rejected (burn-in) the first 2000 samples. The number of steps was greater than 50 times the autocorrelation timescale (computed by emcee), which is generally an indication of convergence (Sokal 1997; Foreman-Mackey et al. 2019). Table C1 lists the prior distributions adopted for the hyperparameters of the dLW and CCF FWHM GP regressions, as well as their posterior’s median, 16th, and 84th percentiles. Figure 11 shows the mean GP predictions for the dLW and the CCF FWHM, using the best-fit hyperparameters.

In Section 5.4, we jointly fit Keplerian signals and activity in the LBL RVs for five different models ($M_{dLW}$, $M_{b+dLW}$, $M_{bd+dLW}$, $M_{bc+dLW}$, and $M_{bcd+dLW}$). We generated Keplerian orbits with radvel (Fulton et al. 2018) and used the same quasiperiodic kernel as described above for the GP activity component. We adopted as prior distributions on the hyperparameters $\ln \ell$, $\ln \Gamma$, and $\ln P_{\text{rot}}$, the posteriors of the GP regression of the dLW time series (Table C1). This is justified because the stellar activity processes affecting the RVs and the dLW are expected to share common properties (correlation timescale, periodicity, etc.) described by the hyperparameters. We sampled the posterior distributions for each M with emcee, using 100 walkers, 50,000 steps, and a burn-in of 5000. The model favored according to a BIC comparison is the single planet $M_{b+dLW}$. We describe the free parameters of model $M_{b+dLW}$ and report the prior and posterior distributions in Table C2. The $M_{bd+dLW}$ model informed on the mass determination significance of Proxima d in the LBL data. In Table C3, we compare our fitted values of P and K for Proxima b and d with the ones obtained in Suárez Mascareño et al. (2020) and Faria et al. (2022). The phase-folded best-fit Keplerian components for planet b and d are presented in Figure C1.

The results obtained with the preferred model $M_{b+dLW}$ are compatible with literature values. We measure a semiamplitude for Proxima b ($K_b = 1.27 \pm 0.17 \text{ m s}^{-1}$) that is fully consistent with Damasso et al. (2020; $K_b = 1.2 \pm 0.1 \text{ m s}^{-1}$), Suárez Mascareño et al. (2020; $K_b = 1.38 \pm 0.10 \text{ m s}^{-1}$), and Faria et al. (2022; $K_b = 1.24 \pm 0.07 \text{ m s}^{-1}$). The increased uncertainty on our $K_b$ is explained by the fact that the LBL time series consists of a subsample of Damasso et al.’s (2020) and Suárez Mascareño et al.’s (2020) data set, focusing on post-fiber change observations with HARPS.
<figure><img src="image_16.png" /><figcaption>**Figure C1.** Phase-folded radial velocity (RV) curves (top panels) and orbital fit residuals (bottom panels) for Proxima b and d (model $M_{bd+dLW}$), with systemic velocity and activity Gaussian process removed. Binned RVs (0.1 phase bin) are represented with black data points. The best-fit Keplerian models for Proxima b and d are shown with solid blue curves, with 68% confidence interval envelopes in light blue.</figcaption></figure>

<table><caption>**Table C1**
Fitted Hyperparameters of the dLW and the CCF FWHM Gaussian Process Regressions</caption><thead><tr><th>Parameter</th><th>Prior</th><th>LBL dLW</th><th>Prior</th><th>CCF FWHM</th></tr></thead><tbody><tr><td>ln A</td><td>U(-5, 15)</td><td>9.0 ± 0.2</td><td>U(-5, 5)</td><td>2.7<sup>+0.3</sup><sub>-0.2</sub></td></tr><tr><td>ln ℓ</td><td>U(3, 10)</td><td>5.3<sup>+0.4</sup><sub>-0.5</sub></td><td>U(3, 10)</td><td>4.9<sup>+0.5</sup><sub>-0.7</sub></td></tr><tr><td>ln Γ</td><td>U(-5, 5)</td><td>0.43<sup>+0.16</sup><sub>-0.17</sub></td><td>U(-5, 5)</td><td>0.27<sup>+0.23</sup><sub>-0.24</sub></td></tr><tr><td>ln P<sub>rot</sub></td><td>U(ln 60, ln 120)</td><td>4.523<sup>+0.045</sup><sub>-0.039</sub></td><td>U(ln 60, ln 120)</td><td>4.38<sup>+0.11</sup><sub>-0.05</sub></td></tr><tr><td>ln s</td><td>U(-5, 15)</td><td>7.15 ± 0.11</td><td>U(-5, 5)</td><td>1.73 ± 0.06</td></tr></tbody></table>

**Note.** The LBL dLW has units of m²/s², compared to m s⁻¹ for the CCF FWHM. The adopted priors for ln A and ln s differ for both indicators because of the different units. For the other parameters, ℓ is in days, Γ is dimensionless, and P<sub>rot</sub> is in days for both data sets.
<table><caption>**Table C2**
Fitted Orbital Parameters for Proxima b and Activity GP (Model $M_{b+\text{dLW}}$)</caption><thead><tr><th>Parameter</th><th>Prior<sup>a</sup></th><th>Posterior</th><th>Description</th></tr></thead><tbody><tr><td colspan="4" style="text-align: center; font-style: italic;">Fitted parameters</td></tr><tr><td>γ (m s<sup>-1</sup>)</td><td>U(-22002, -21982)</td><td>-21992.2 ± 0.5</td><td>Systemic velocity component</td></tr><tr><td>T<sub>0,b</sub> (BJD-2450000)</td><td>U(7895, 7910)</td><td>7898.5 ± 0.4</td><td>Time of inferior conjunction</td></tr><tr><td>P<sub>b</sub> (days)</td><td>U(<b>10.0</b>, <b>20.0</b>)</td><td>11.1888<sup>+0.0060</sup><sub>-0.0059</sub></td><td>Orbital period</td></tr><tr><td>K<sub>b</sub> (m s<sup>-1</sup>)</td><td>U(0, 5)</td><td>1.27 ± 0.17</td><td>RV semiamplitude</td></tr><tr><td>h<sub>b</sub> = &radic;<span style="text-decoration: overline;">e<sub>b</sub></span> cos ω<sub>b</sub></td><td>U(-1, 1)</td><td>-0.06<sup>+0.25</sup><sub>-0.23</sub></td><td>Parameterization of e and ω</td></tr><tr><td>k<sub>b</sub> = &radic;<span style="text-decoration: overline;">e<sub>b</sub></span> sin ω<sub>b</sub></td><td>U(-1, 1)</td><td>0.08<sup>+0.29</sup><sub>-0.28</sub></td><td>Parameterization of e and ω</td></tr><tr><td colspan="4" style="text-align: center; font-style: italic;">Derived parameters</td></tr><tr><td>M<sub>p,b</sub> sin i (M<sub>⊕</sub>)<sup>b</sup></td><td>...</td><td>1.07 ± 0.14</td><td>Minimal mass</td></tr><tr><td>e<sub>b</sub></td><td>...</td><td>0.10<sup>+0.13</sup><sub>-0.07</sub></td><td>Orbital eccentricity</td></tr><tr><td>ω<sub>b</sub> (rad)</td><td>...</td><td>0.8<sup>+1.5</sup><sub>-2.7</sub></td><td>Argument of periastron</td></tr><tr><td colspan="4" style="text-align: center; font-style: italic;">Activity GP parameters</td></tr><tr><td>ln A</td><td>U(-5, 5)</td><td>0.35<sup>+0.25</sup><sub>-0.24</sub></td><td>Amplitude of the GP</td></tr><tr><td>ln ℓ</td><td>GP(dLW)</td><td>5.3<sup>+0.4</sup><sub>-0.5</sub></td><td>Coherence length of the GP</td></tr><tr><td>ln Γ</td><td>GP(dLW)</td><td>0.52 ± 0.14</td><td>Periodic component scale</td></tr><tr><td>ln P<sub>rot</sub></td><td>GP(dLW)</td><td>4.515<sup>+0.032</sup><sub>-0.037</sub></td><td>Stellar rotation period</td></tr><tr><td>ln s</td><td>U(-5, 5)</td><td>0.12 ± 0.09</td><td>Extra white noise term</td></tr></tbody></table>

**Notes.**

<sup>a</sup> GP(dLW) is the posterior distribution on the hyperparameter from the dLW GP regression.

<sup>b</sup> Using Mann et al.'s (2015) estimate of Proxima's mass ($M_{\star} = 0.1221 \pm 0.0022 M_{\odot}$) calculated in Suárez Mascareño et al. (2020).

<table><caption>**Table C3**
Comparison Between The LBL Detection of Proxima b and d (Model $M_{bd+\text{dLW}}$) with Published Literature Values</caption><thead><tr><th>Parameter</th><th>Prior</th><th><sup>a</sup>LBL</th><th><sup>b</sup>SM2020</th><th><sup>c</sup>Fa2022</th></tr></thead><tbody><tr><td>P<sub>b</sub> (days)</td><td>U(<b>10.0</b>, <b>20.0</b>)</td><td>11.1881<sup>+0.0061</sup><sub>-0.0058</sub></td><td>11.1842 ± 0.0007</td><td>11.1868<sup>+0.0029</sup><sub>-0.0031</sub></td></tr><tr><td>K<sub>b</sub> (m s<sup>-1</sup>)</td><td>U(<b>0</b>, 5)</td><td>1.22 ± 0.17</td><td>1.37 ± 0.10</td><td>1.24 ± 0.07</td></tr><tr><td>P<sub>d</sub> (days)</td><td>U(<b>5.0</b>, <b>5.3</b>)</td><td>5.167<sup>+0.047</sup><sub>-0.091</sub></td><td>5.168<sup>+0.051</sup><sub>-0.069</sub></td><td>5.122<sup>+0.002</sup><sub>-0.036</sub></td></tr><tr><td>K<sub>d</sub> (m s<sup>-1</sup>)</td><td>U(<b>0</b>, 5)</td><td>0.38<sup>+0.19</sup><sub>-0.20</sub></td><td>0.35<sup>+0.10</sup><sub>-0.11</sub></td><td>0.39 ± 0.07</td></tr></tbody></table>

**Notes.**

<sup>a</sup> This work, from publicly available HARPS data (post-fiber change, N = 181).

<sup>b</sup> Suárez Mascareño et al. (2020) values, Table 3 (GP+2 Signals) therein, from UVES ($N_{\text{UVES}} = 63$), HARPS ($N_{\text{HARPS}} = 196$), and ESPRESSO ($N_{\text{ESPRESSO}} = 77$) data.

<sup>c</sup> Faria et al. (2022) values, from their Table C.1 using the template-matching RVs from ESPRESSO (N = 114).

## ORCID iDs

Étienne Artigau https://orcid.org/0000-0003-3506-5667

Charles Cadieux https://orcid.org/0000-0001-9291-5555

Neil J. Cook https://orcid.org/0000-0003-4166-4121

René Doyon https://orcid.org/0000-0001-5485-4675

Thomas Vandal https://orcid.org/0000-0002-5922-8267

Jean-François Donati https://orcid.org/0000-0001-5541-2887

Claire Moutou https://orcid.org/0000-0002-2842-3924

Xavier Delfosse https://orcid.org/0000-0001-5099-7978

Pascal Fouqué https://orcid.org/0000-0002-1436-7351

Eder Martioli https://orcid.org/0000-0002-5084-168X

François Bouchy https://orcid.org/0000-0002-7613-393X

Jasmine Parsons https://orcid.org/0000-0002-6013-4655

Andres Carmona https://orcid.org/0000-0003-2471-1299

Xavier Dumusque https://orcid.org/0000-0002-9332-2011

Nicola Astudillo-Defru https://orcid.org/0000-0002-8462-515X

Xavier Bonfils https://orcid.org/0000-0001-9003-8894

Lucille Mignon https://orcid.org/0000-0002-5407-3905

## References

Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O'Neil, M. 2016, *ITPAM*, 38, 252

Anglada-Escudé, G., & Butler, R. P. 2012, *ApJS*, 200, 15

Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, *Natur*, 536, 437

Angus, R., Morton, T., Aigrain, S., Foreman-Mackey, D., & Rajpaul, V. 2018, *MNRAS*, 474, 2094

Artigau, É., Doyon, R., Hernandez, O., et al. 2018a, *Proc. SPIE*, 10709, 107091P

Artigau, É., Malo, L., Doyon, R., et al. 2018b, *AJ*, 155, 198

Artigau, É., Bouchy, F., Delfosse, X., et al. 2012, *Proc. SPIE*, 8451, 84513I

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

Artigau, É., Hébrard, G., Cadieux, C., et al. 2021, *AJ*, 162, 144

Astudillo-Defru, N., Forveille, T., Bonfils, X., et al. 2017, *A&A*, 602, A88

Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, *A&A*, 577, A42
Baranne, A., Mayor, M., & Poncet, J. L. 1979, VA, 23, 279

Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373

Barbary, K. 2021, nestle: Nested sampling algorithms for evaluating Bayesian evidence, Astrophysics Space Code Library, ascl:2103.022

Baril, M. R., & Albert, L. 2008, Proc. SPIE, 9149, 702121

Barnard, E. E. 1918, AJ, 3, 191

Bedell, M., Hogg, D. W., Foreman-Mackey, D., Montet, B. T., & Luger, R. 2019, AJ, 158, 164

Bellotti, S., Petit, P., Morin, J., et al. 2022, A&A, 657, A107

Benedict, G. F., & McArthur, B. E. 2020, RNAAS, 4, 86

Benedict, G. F., McArthur, B., Nelan, E., et al. 1998, AJ, 116, 429

Bertaux, J. L., Lallemand, R., Ferron, S., Boonne, C., & Bodichon, R. 2014, A&AS, 564, A46

Bishop, C. M. 2006, Pattern Recognition and Machine Learning, Information Science and Statistics (New York: Springer),

Blind, N., Conod, U., & Wildi, F. 2017, arXiv:1711.00835

Bonfils, X., Lo Curto, G., Correia, A. C. M., et al. 2013, A&A, 556, A110

Bouchy, F., Pepe, F., & Queloz, D. 2001, A&A, 374, 733

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

Claudi, R., Benatti, S., Carleo, I., et al. 2017, EPJP, 132, 364

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

Connes, P., Martic, M., & Schmitt, J. 1996, Ap&SS, 241, 61

Courcol, B., Bouchy, F., Pepe, F., et al. 2015, A&A, 581, A38

Cretignier, M., Dumusque, X., Allart, R., Pepe, F., & Lovis, C. 2020, A&A, 633, A76

Cretignier, M., Dumusque, X., Hara, N. C., & Pepe, F. 2021, A&A, 653, A43

Da Costa-Luis, C., Larroque, S. K., Altendorf, K., et al. 2021, tqdm: A fast, Extensible Progress Bar for Python and CLI, Zenodo, doi:10.5281/ ZENODO.5517697

Damasso, M., Del Sordo, F., Anglada-Escudé, G., et al. 2020, SciA, 6, eaax7467

Delfosse, X., Bonfils, X., Forveille, T., et al. 2013, A&A, 553, A8

Delmotte, N., Dolensky, M., Padovani, P., et al. 2006, in ASP Conf. Ser. 351, Astronomical Data Analysis Software and Systems XV (San Francisco, CA: ASP), 690

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

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

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

Figueira, P., Adibekyan, V. Z., Oshagh, M., et al. 2016, A&A, 586, A101

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Foreman-Mackey, D., Farr, W., Sinha, M., et al. 2019, JOSS, 4, 1864

Fried, D. L. 1966, JOSSA, 56, 1372

Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504

Gaia Collaboration, Smart, R. L., Sarro, L. M., et al. 2021, A&A, 649, A6

Gan, T., Soubkiou, A., Wang, S. X., et al. 2022, MNRAS, 514, 4120

Gratton, R., Zurlo, A., Le Coroller, H., et al. 2020, A&A, 638, A120

Gullikson, K. 2014, Telfit: Fitting the telluric absorption spectrum, Astrophysics Space Code Library, ascl:1405.002

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hatzes, A. P. 2019, The Doppler Method for the Detection of Exoplanets (Bristol: IOP Publishing)

Hobson, M. J., Bouchy, F., Cook, N. J., et al. 2021, A&A, 648, A48

Hogg, D. W., Bovy, J., & Lang, D. 2010, arXiv:1008.4686

Horne, K. 1986, PASP, 98, 609

Hunter, J. D. 2007, CSE, 9, 90

Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

Innes, R. T. A. 1915, CiUO, 30, 235

Irwin, J., Charbonneau, D., Nutzman, P., & Falco, E. 2009, in IAU Symp. 253, Transiting Planets (Cambridge: Cambridge Univ. Press), 37

Joye, W. A., & Mandel, E. 2003, in ASP Conf. Ser. 295, Astronomical Data Analysis Software and Systems XII, Vol. 295 ed. H. E. Payne et al. (San Francisco, CA: ASP), 489

Kervella, P., Arenou, F., & Schneider, J. 2020, A&A, 635, L14

Kiraga, M., & Stepien, K. 2007, AcA, 57, 149

Klein, B., Donati, J.-F., Hébrard, É. M., et al. 2020, MNRAS, 500, 1844

Kotani, T., Tamura, M., Suto, H., et al. 2014, Proc. SPIE, 1947, 914714

Lafarga, M., Ribas, I., Lovis, C., et al. 2020, A&A, 636, A36

Lam, S. K., Pitrou, A., & Seibert, S. 2015, in SC15: The International Conference for High Performance Computing, Networking, Storage and Analysis (New York: ACM), 1

Lindegren, L., & Dravins, D. 2003, A&A, 401, 1185

Lo Curto, G., Pepe, F., Avila, G., et al. 2015, Msngr, 162, 9

Lubin, J., Robertson, P., Stefansson, G., et al. 2021, AJ, 162, 61

Mahadevan, S., Ramsey, L., Bender, C., et al. 2012, Proc. SPIE, 8446, 84461S

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

Martioli, E., Hébrard, G., Fouqué, P., et al. 2022, A&A, 660, A86

Mayor, M., & Queloz, D. 1995, Natur, 378, 355

Micheau, Y., Challita, Z., Gallou, G., et al. 2018, Proc. SPIE, 10702, 107025R

Mirabet, E., Carvas, P., Lizon, J.-L., et al. 2014, Proc. SPIE, 9151, 91513Y

Newton, E. R., Mondrik, N., Irwin, J., Winters, J. G., & Charbonneau, D. 2018, AJ, 156, 217

Pepe, F., Mayor, M., Rupprecht, G., et al. 2002, Msngr, 110, 9

Perez, F., & Granger, B. E. 2007, CSE, 9, 21

Pojmanski, G. 1997, AcA, 47, 467

Quirrenbach, A., Amado, P. J., Ribas, I., et al. 2018, Proc. SPIE, 10702, 107020W

Reiners, A., Bean, J. L., Huber, K. F., et al. 2010, ApJ, 710, 432

Ribas, I., Tuomi, M., Reiners, A., et al. 2018, Natur, 563, 365

Rousselot, P., Lidman, C., Cuby, J.-G., Moreels, G., & Monnet, G. 2000, A&A, 354, 1134

Santos, N. C., Mayor, M., Naef, D., Queloz, D., & Udry, S. 2001, in IAU Symp. 202, Planetary Systems in the Universe (Cambridge: Cambridge Univ. Press), 121

Schwarz, G. 1978, AnSta, 6, 461

SciPy 1.0 Contributors, Virtanen, P., Gommers, R., et al. 2020, NatMe, 17, 261

Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77

Sokal, A. 1997, in Functional Integration, ed. C. DeWitt-Morette, P. Cartier, & A. Folacci, Vol. 361 (New York: Springer), 131

Stefansson, G., Hearty, F., Robertson, P., et al. 2016, ApJ, 833, 175

Suárez Mascareño, A., Rebolo, R., & González Hernández, J. I. 2016, A&A, 595, A12

Suárez Mascareño, A., Faria, J. P., Figueira, P., et al. 2020, A&A, 639, A77

The Astropy Collaboration, Robitaille, E. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

The Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Vacca, W., Cushing, M., & Rayner, J. 2003, PASP, 115, 389

Wargelin, B. J., Saar, S. H., Pojmański, G., Drake, J. J., & Kashyap, V. L. 2017, MNRAS, 464, 3281

Zechmeister, M., Reiners, A., Amado, P. J., et al. 2018, A&A, 609, A12