© The Authors 2025

1 Introduction

Doppler spectroscopy, also known as the radial velocity (RV) method, has been one of the most successful techniques used in exoplanet detection and characterisation (e.g. Mayor et al. 2014). Due to the advent of extremely high-precision instrumentation such as the Echelle SPectrograph for Rocky Exoplanet and Stable Spectroscopic Observations (ESPRESSO) (Pepe et al. 2021), this method is currently capable of detecting planets with RV semi-amplitudes below the 1 m/s level (e.g. Faria et al. 2022). Although instrumental precision has reached the level of ~10 cm/s (similar to the signal caused by Earth orbiting the Sun), this level is particularly difficult to achieve in practice due to the effects caused on RV by stellar variability, or stellar noise, which cannot be reduced with better instrumentation.

Sources of stellar noise can be segregated into different timescales. Stellar p-mode oscillations are characterised by the propagation of acoustic waves in the convective envelopes of late-type stars, which translates into RV variability (e.g. Kjeldsen & Bedding 1995; Bedding et al. 2001; Arentoft et al. 2008; Kunovac Hodžić et al. 2021). In solar-type stars this phenomenon is on timescales of minutes (Christensen-Dalsgaard 2004; Kjeldsen et al. 2005), and their effects can be easily attenuated using exposure times with longer timescales or by averaging multiple observations during a given night (Dumusque et al. 2011). The amplitude of oscillations decreases with decreasing stellar effective temperature, rendering this effect negligible in cooler late-type stars. Granulation, and its larger counterpart supergranulation, are caused by the convection pattern in the stellar photosphere and have timescales from hours to a few days (De Rosa & Toomre 2004; Brandt & Getling 2008; Meunier et al. 2015; Cegla et al. 2018). Granulation is also temperature dependent and its influence on RV decreases towards the later-type stars (Meunier et al. 2017a,b; Liebing et al. 2021). The effects of granulation can also be reduced by careful observational strategies (Dumusque et al. 2011; Meunier et al. 2017b), although the impact of supergranulation on RV is still poorly understood. Both p-mode oscillations and granulation induce RV variability up to the m/s level (e.g. Arentoft et al. 2008; Bedding et al. 2001; Meunier et al. 2015).

Stellar activity, in the form of active regions such as spots and faculae, affect RV on timescales typical of the rotation period and its harmonics (Saar & Donahue 1997; Santos et al. 2000; Desort et al. 2007; Lagrange et al. 2010; Meunier et al. 2010; Boisse et al. 2011; Reiners et al. 2013). The signals depend on the number, the latitudinal and longitudinal distribution of the active regions, and the stellar inclination, making them quasi-periodic (e.g. Angus et al. 2018) and possibly emulating planetary signals (e.g. Bonfils et al. 2007; Huélamo et al. 2008; Figueira et al. 2010; Santos et al. 2014; Robertson et al. 2014; Carmona et al. 2023). Activity-induced RV variability can reach hundreds of m/s (e.g. Donati et al. 2024), well above the instrumental precision. Furthermore, the evolution of active regions on timescales of years and decades, in the form of activity cycles (Baliunas et al. 1995), can also produce periodic RV variability of the order of tens of m/s in solar-type stars (Lovis et al. 2011) and M dwarfs (Gomes da Silva et al. 2012).

Spectral lines can be affected by stellar activity in different forms, depending on the type of active region, the magnetic field strength, the disk position, and the star spectral type and activity level. For example, the photometric effect of spots, caused by their temperature contrast to the quiet photosphere, when modulated by stellar rotation, distorts the disk-integrated spectral lines by changing the balance between the blueshifted approaching half of the disk and the redshifted receding half, causing a periodic line asymmetry that affects the centroid position, and therefore the measured RV (e.g. Reiners et al. 2010; Hébrard et al. 2014). The photometric effect is also known to be wavelength-dependent. There is evidence that its influence on RV decreases towards redder wavelengths where the spot and faculae contrasts with the photosphere are expected to decrease (e.g. Reiners et al. 2010; Carmona et al. 2023). However, this dependence on wavelength is primarily driven by contrast, with lower contrasts exhibiting a more pronounced decrease in amplitude with wavelength.

Another effect caused by active regions is the inhibition of convective blueshift. Convective motions in the photosphere contribute to the granulation pattern where hot plasma rises from the stellar interior in bright granules, cools down, and sinks in the intergranular lanes surrounding the granules. Since the areas in the Sun covered by granules are larger than the intergranular lanes, the blueshifted flux from these bright regions is greater than the redshifted flux from the intergranular lanes. As a result, the overall effect is a net blueshift of the spectral lines (Dravins et al. 1981). This effect affects spectral lines depending on their depth. Shallower lines that formed in hotter regions deeper in the photosphere are more strongly affected (Reiners et al. 2016; Meunier et al. 2017a). Active regions are marked by intense magnetic fields (e.g. Schrijver & Zwaan 2000), which disable the convection motions and attenuate the convective blueshift. This disruption alters the balance across the stellar disk, leading to RV variability (e.g. Meunier et al. 2010). Since the convective blueshift decreases with the effective temperature of stars, this effect is stronger for solar-type stars, but is reduced for M dwarfs (Liebing et al. 2021).

Generally, the identification of activity signals is performed using activity-sensitive spectral lines, such as Ca II H&K, Hα, Na I, and others, whose rotational modulation can be analysed using tools such as periodograms (e.g. Zechmeister & Kürster 2009; VanderPlas 2018). The signals observed in these lines can also be used as anchors when modelling the activity contribution to RV variability (e.g. Rajpaul et al. 2015; Faria et al. 2020, 2022). Other techniques, such as the cross-correlation function (CCF), can harvest the information of multiple spectral lines into an average line. The parameters of the CCF profile, such as RV, full width at half maximum (FWHM), contrast, and the bisector inverse slope (BIS), are all sensitive to line profile distortions caused by stellar activity (e.g. Queloz et al. 2001; Santos et al. 2000, 2003, 2010; Boisse et al. 2011; Gomes da Silva et al. 2012; Figueira et al. 2014; Dumusque et al. 2014). However, these methods do not use all the information available in the stellar spectrum. Although activity-sensitive lines might be more affected by the strong magnetic fields in active regions, all spectral lines are affected by stellar activity variability. Furthermore, methods that combine spectral lines into an average line profile end up averaging the activity effects on those lines if they affect the lines differently. Therefore, to better understand the activity effects on RV, it is preferable to include most of the information available in the stellar spectrum by analysing how the activity affects most of the lines (Giguere et al. 2016).

Recent studies have used the full spectrum information to analyse how line profiles and their RVs are affected by activity. Most of them explored the spectra of GK dwarfs, mainly that of alpha Cen B. Thompson et al. (2017) identified line deformations caused by activity; several lines showed depth variations and FWHM variability modulated by stellar rotation and correlated with the $R_{HK}^{\mathrm{\prime }}$ chromospheric emission ratio. A list of approximately 40 absorption lines showing line depth correlation with the Ca II H&K activity index and rotation modulation was provided by Wise et al. (2018) after studying line profile variability. An extension of this work by Ning et al. (2019) identified lines sensitive to activity via an automated method based on Bayesian variable selection using the activity indices based on the Ca II H&, Na I, and Hα lines and the CCF BIS and FWHM. Equivalent widths and line asymmetries were measured by Lisogorskyi et al. (2019) and compared to the Ca II H&K indicator, who found approximately 350 lines with strong correlations. Dumusque (2018) showed that by selecting activity-sensitive or -insensitive lines to compute the overall RV, the activity signal could be enhanced or attenuated. Using line-by-line (LBL) RVs, a relationship between line depth and RV activity sensitivity was detected by Cretignier et al. (2020), where deeper lines show weaker RV activity effects. This can be explained by the shallower lines closer to the photosphere that have higher sensitivity to the inhibition of convection caused by active regions. Al Moulla et al. (2022, 2024) found that measuring RV at different line heights or formation temperatures provides different correlations with activity and different periodicities, meaning that activity influence on RV depends on the line formation temperature. At higher temperatures, the RV is strongly correlated with inhibition of convection and the Ca II H&K index, while at cooler temperatures (higher in the atmosphere) the RV is anticorrelated to the global RV and the Ca II H&K index. Recently, Artigau et al. (2024) introduced a method for computing stellar temperature perturbations on individual lines using a library of high-resolution High-Accuracy Radial-velocity Planet Searcher (HARPS) spectra across various spectral types.

In the only study to date focusing on M dwarfs, Lafarga et al. (2023) used a similar approach to Dumusque (2018), but applied it to a sample of six early to mid-M dwarfs observed with the visual arm of the high-resolution spectrograph Calar Alto High-Resolution Search for M dwarfs with Exoearths with Near-infrared and Optical Échelle Spectrographs (CARMENES) (Quirrenbach et al. 2016, 2018). They found that there is a different activity dependence for the same lines in different stars of similar spectral type, activity level, and rotation period. They further used LBL RVs insensitive to activity to decrease the RV scatter by around two to five times in their sample stars.

For this work, we analysed the influence of activity on the near-infrared (NIR) spectra of two M dwarfs, Gl 581 and Proxima. From the perspective of detecting and characterising exoplanets, these cool, small stars offer many advantages over their larger solar-type counterparts. M dwarfs are abundant, comprising approximately 75% of all stars in the solar neighbourhood (Henry et al. 1996). This means that most of the planets closer to the solar neighbourhood are orbiting these cool stars. Their intrinsic low mass makes it easier to detect the RV signals induced by planets, thereby increasing the precision to which lower mass or longer orbit planets can be detected. Furthermore, their lower luminosity also contributes to a habitable zone closer to the star, enabling the easier detection and characterisation of potentially habitable planets. Another advantage over FGK dwarfs is that the RV contribution from convective blueshift (and its attenuation by active regions) is very low (Beeck et al. 2013; Baroch et al. 2020; Liebing et al. 2021) as is the effect of stellar oscillations. Moreover, activity induced RV amplitude decreases with wavelength, and is expected to be lower at redder and infrared wavelengths (Desort et al. 2007; Reiners et al. 2010) where these stars have a higher fraction of their fluxes.

This article is organised as follows. In Sect. 2, we present the sample stars. A description of the HARPS and Near-infrared Planet Searcher (NIRPS) data is provided in Sect. 3. In Sect. 4, the known activity indicators are analysed along with the detection of the rotation periods. Our methodology, including the extraction of line profile parameters and the calculation of different pseudo-equivalent widths, is presented in Sect. 5 and our results are discussed in Sect. 6. We conclude our findings in Sect. 7.

2 Target stars

In this work, our aim is to study the effects of stellar activity and to detect activity-sensitive lines in the spectra of M dwarfs in the NIR regime. We are interested in finding lines that can be used as activity indicators for stars with different spectral types and activity levels. This type of analysis is also computer intensive because we deal with modelling thousands of spectral lines; we therefore limited our study to two stars. We are interested in: (1) stars at the two sides of the limit to fully convective regime, at around spectral type M3–4 (Siess et al. 2000), representing two cases of M dwarfs with different internal structures; (2) stars with different activity levels, which we identified visually as having the Hα line in emission or absorption; (3) stars with a fair amount of observations and obvious rotation modulations in the activity indicators provided by the Data & Analysis Center for Exoplanets (DACE)¹, namely in the Hα index and the FWHM of the NIRPS CCF.

We therefore selected the planet-hosting stars Gl 581 (M3V) and Proxima (M5.5V), since both present a high number of observations with NIRPS and HARPS, a good signal-to-noise ratio (S/N) in the NIRPS spectra (>100), and a clear rotational modulation signature in the activity indices (see Sect. 4). Proxima is a fully convective, active star with emission in both Hα and Na I lines while Gl 581 shows absorption in both lines, meaning that its activity level is below that of Proxima. The Hα line is known to first increase in absorption and then in emission as the activity level increases (Cram & Mullan 1979; Stauffer & Hartmann 1986). This behaviour could cause confusion when assessing the stellar activity levels. However, since at high activity levels the emission is above the continuum, we are able to discriminate between high activity levels (strong emission) and low to medium activity levels (absorption). Nevertheless, to confirm the activity levels of both stars, we estimated their log $R_{HK}^{\mathrm{\prime }}$ using the Ca II H&K lines extracted from the HARPS spectra and the calibration provided by Suárez Mascareño et al. (2015), which includes M dwarfs. We note that, for M dwarfs, the flux in the Ca II H&K is very low, and after background subtraction, can reach unphysical negative values, which will bias the index towards lower activity levels. We obtained values of −5.65 ± 0.01 dex and −5.52 ± 0.01 dex for Gl 581 and Proxima, respectively, showing that Proxima has a higher activity level. The basic stellar parameters of both stars are presented in Table 1.

Along with finding spectral lines whose profile variability shows correlations with known activity indices, we are also interested in finding the rotation period of the stars using those lines, and therefore we need to know their rotation periods. Several estimates of Gl 581’s rotation period can be found in the literature. Robertson et al. (2014) estimated a rotation period of 130 ± 2.0 days for Gl 581 using the time series of the Hα index. Suárez Mascareño et al. (2015) used 251 Ca II H&K HARPS (Mayor et al. 2003) observations of this star, with a time baseline of almost 8 years, and obtained a rotation period of 132.5 ± 6.3 days. More recently, Engle & Guinan (2023) found a rotation period of 148.1 ± 0.9 days using photometric observations. In the case of Proxima, Klein et al. (2021) detected a rotation period of 90 ± 4 days using 10 observations of the longitudinal magnetic field, and the Hα, Na I, and He I activity indices extracted with the HARPS-POL (the extension of HARPS for polarimetry observations) spectropolarimeter (Mayor et al. 2003; Snik et al. 2011). Faria et al. (2022) used Gaussian processes to simultaneously model ESPRESSO (Pepe et al. 2021) RVs and CCF FWHMs, deriving a rotation period of ${85.1}_{-0.8}^{+0.9}$ days. A rotation period of 86.4 ± 2.5 days was measured by Engle & Guinan (2023) using photometry. More recently, Suárez Mascareño et al. (2025) measured a rotation period of 83.2 ± 1.6 days using the same dataset presented here.

The above rotation periods were derived using different activity proxies, number of observations, time spans, at different epochs, and using different methods, which contribute to the different values obtained by each reference. We use these reported rotation periods as reference in our search of rotational variability with the known activity indicators, based on activity-sensitive lines in the visual and the CCF parameters (see Sect. 4).

3 Observations and data

3.1 HARPS and NIRPS spectroscopy

We used data from the NIRPS spectrograph (Bouchy et al. 2025) obtained simultaneously with visual spectra from HARPS (Mayor et al. 2003), both installed on the ESO 3.6 m telescope, located at the La Silla Observatory in Chile. NIRPS is a NIR (0.97–1.81 μm), fibre-fed, stabilised, high-resolution echelle spectrograph assisted by adaptive optics. HARPS is an optical (330–690 nm), fiber-fed, stabilised high-resolution echelle spectrograph.

Both Gl 581 and Proxima were observed as part of Work Package one (WP1) of the NIRPS consortium guarantee time of observation (GTO) under programmes 111.251P.001, 112.25NZ. 001 and 112.25NZ. 002 (PI: F. Bouchy). Proxima was observed to demonstrate the potential of NIRPS and to highlight the advantages of using NIR data to search for planets around M dwarfs, and to refine the orbital parameters of its planetary system (Suárez Mascareño et al. 2025). Meanwhile, Gl 581 was selected as a NIRPS target to assess the improvement in RV precision for characterising a planetary system whose central star is known to have an activity signal that is significant for detecting planets in the habitable zone (Larue et al., in prep.). A total of 82 (Gl 581) and 137 (Proxima) observations obtained with the high-efficiency (HE) mode (R~75 200) were carried out between 2024/04/02 and 2024/09/15 with a median S/N of 156 (Gl 581) and 187 (Proxima) for NIRPS and 20 (Gl 581) and 6 (Proxima) for HARPS. These observations were reduced using an adaptation of the ESPRESSO Data Reduction Software (DRS) pipeline (Pepe et al. 2021) for NIRPS and HARPS with versions 3.2.0 and 3.0.0 for Gl 581 and 3.2.0 for both instruments in the case of Proxima. We used order by order telluric corrected spectra as provided by the DRS. The CCF was computed by passing the spectra through a weighted cross-correlation (Baranne et al. 1996; Pepe et al. 2002) with a numerical mask with a step size equal to instrument pixel width. The cross-correlation masks used for Gl 581 were M4V for NIRPS and M3V for HARPS while for Proxima the M5V mask was used for both instruments. For a description of the telluric correction process, we refer to Allart et al. (2022) and Bouchy et al. (2025). All spectra were shifted to the stellar rest frame using the CCF RV and the Earth barycentric radial velocity (BERV) provided by the DRS. The observation log is provided in Table 2. All spectra were downloaded from DACE.

3.2 Activity indicators

Stellar variability caused by active regions can be estimated via the effects they cause on spectral lines by measuring their flux variations or variability in their widths, depths and line bisectors. To assess the activity variability of our stars, we extracted the CCF parameters FWHM, contrast, and BIS (Queloz et al. 2001) from the HARPS and NIRPS DRS, along with the activity indices Ca II H&K, Hα06, Hα16, and Na I, using ACTIN v2.0_beta14² (Gomes da Silva et al. 2018, 2021). The CCF parameters are known to be sensitive to activity variations of cool main sequence FGK (Lovis et al. 2011), M (Gomes da Silva et al. 2012), and evolved stars (Delgado Mena et al. 2018, 2023). We extracted two versions of the Hα index, with 0.6 Å (Hα06) and 1.6 Å (Hα16) central bandpasses to compare their behaviour for M-dwarfs. It was shown that Hα06 is better correlated with the Ca II H&K index for FGK dwarfs than Hα16 (Gomes da Silva et al. 2022), but there are indications their behaviour is similar for M-dwarfs (Meunier et al. 2022). The Ca II H&K, Hα16, and Na I indices were analysed in Gomes da Silva et al. (2011) for a sample of 30 M-dwarfs using HARPS spectra from the M-dwarf survey (Bonfils et al. 2007). It was shown that Na I has a strong correlation with Ca II H&K, while Hα16 is only strongly correlated for the more active stars, with the correlation dropping to zero or negative values as the Ca II H&K activity level decreases. As will be discussed in the following section, we ended up not using the Ca II H&K indicator in our analysis due to the very low S/N of the HARPS observations, which affects mostly the blue end of the spectrum of M dwarfs, where the Ca II H&K are located.

For each time series, a sequential 3σ-clipping from the median for each parameter and respective errors was carried out to remove outliers due to poor telluric correction, weather and instrumental issues, or stellar flares that could complicate our analysis. All data were binned per day to average out high-frequency variability.

4 Activity analysis with known indicators

To assess the activity sensitivity of the NIR lines, we need to analyse the known activity indicators behaviour, inter-correlations, and ability to detect the rotation modulation of each star with the available data cadence, time span, and S/N. By ‘known indicators’ we mean the activity proxies most analysed for M dwarfs, which includes the ones based on the Ca II H&K, Hα, Na I, and He I lines (e.g. Gomes da Silva et al. 2011) and the CCF parameters FWHM, contrast and BIS (e.g. Queloz et al. 2001; Gomes da Silva et al. 2012). The correlation between these different chromospheric indices and the RV activity signal is studied in more detail in a companion paper by the NIRPS consortium (Larue et al., in prep.).

We used the generalised Lomb–Scargle (GLS) periodogram (Zechmeister & Kürster 2009) as implemented by astropy (Astropy Collaboration 2022) to search for periods longer than 2.5 days (to exclude daily aliases) and shorter than the time span of observations of each star. The highest periodogram peaks, more significant than the 0.1% false alarm probability (FAP) level, and close to the literature rotation periods (see Sect. 2) were selected as a rotation period detection. We do not expect to find exactly the same periods as the reported ones since the data, methods and epochs of detection used are not the same as the ones used here.

Figures A.1a and A.1b show the GLS periodograms for the activity indicators. The time series will be provided in Larue et al., in prep. For Proxima, the first season of observations shows a period of 100 days, while in the second season the period is 80.6 days. This is observed in the GLS, where we have two significant peaks at ~78 and ~98 days surrounding the literature periods. Faria et al. (2022) showed that the evolution timescale of active regions in Proxima is ${151}_{-27}^{+38}$ days, a duration corresponding to around two rotation periods. More recently, as part of the NIRPS consortium and using the same dataset used here, Suárez Mascareño et al. (2025) obtained an evolution timescale of ${58.6}_{-1.2}^{+1.2}$ days, less than the rotation period of the star. Our data covers 532 days or ~6 rotation periods, so we expect to observe a change in the temporal evolution of active regions during our observing time. This variability could also contribute to complex behaviour in our activity time series. Thus, the emergence and disappearance of active regions, at different disk latitudes, can be the cause of the period change observed between the first and second seasons in the time series of this star (see Suárez Mascareño et al. 2025). However, the lower S/N in the HARPS observations for Proxima, with their increased noise, could also explain the double peak observed in the spectral indices for this star.

We found a strong correlation between the NIRPS CCF FWHM and BERV for Gl 581, with a Pearson coefficient of ρ = 0.6. When comparing the NIR spectral line parameters with the NIRPS CCF FWHM, we need to be careful to identify if the correlation is being driven by BERV (see Sect. 6.3). This trend was not observed for Proxima, nor for the other activity indicators.

4.1 Comparison of activity indices for Gl 581 and Proxima

For these M dwarfs, and in agreement with Meunier et al. (2022), the two versions of the Hα index have the same behaviour, both detecting similar rotation periods of 77.4 and 77.5 days for Proxima and 140.6 and 139.7 days for Gl 581, while having Pearson correlation coefficients of ρ = 0.99 and ρ = 0.98 for each star. This indicates that the two indices are being affected by the same activity phenomena crossing the stellar disk, contrary to what happens for FGK dwarfs where very often the indices show no correlation or are anti-correlated (Gomes da Silva et al. 2014, 2022). Because of this, and if confirmed for more M dwarfs, the Hα06 should be the Hα version of the index used for detecting activity for M and FGK dwarfs. Using the same bandpass will enable the direct comparison between the activity levels of these two groups of stars (after applying the respective bolometric corrections). Hereafter we use the Hα06 index when we refer to Hα.

It is known that Gl 581 exhibits an anti-correlation between the Hα and Na I emission line (ρ = −0.75; Gomes da Silva et al. 2011) but a positive correlation between NaI and the HARPS CCF FWHM (ρ = 0.60; Gomes da Silva et al. 2012), and we confirm that we observe it in this state. Concerning Proxima, and consistent with what was observed by Gomes da Silva et al. (2011) with the HARPS CCF FWHM, we detect a strong positive correlation between the NIRPS CCF FWHM, Hα and Na I. This different behaviour between the indicators could be a consequence of either the activity level (see Gomes da Silva et al. 2011) or the star internal structure, including the existence of a tachocline for Gl 581 but not for Proxima, which could be a source of different dynamos operating in the stars resulting in distinctive magnetic field topology in the stellar disks (e.g. Siess et al. 2000; Ossendrijver 2003).

The HARPS CCF parameters of both stars have high noise levels due to the low S/N obtained for this instrument. This is even more evident in the case of BIS, since its sensitivity decreases for slow rotating stars as shown by Santos et al. (2003).

The fact that we detect slightly different rotation periods for several activity indicators can be due to different activity indices being affected by noise, but also due to different activity phenomena, such as faculae and plages, spots, and filaments (e.g. Meunier & Delfosse 2009; Gomes da Silva et al. 2022), by different amounts. If these activity phenomena are not spatially correlated in latitude and the star has differential rotation, it is expected that different activity tracers will deliver slightly different rotation periods.

4.2 Gl 581 rotation

According to the discussion in Sect. 2, the rotation period is expected to be in the range of 130–148 days. We detect significant (FAP ≤ 0.1%) periods in this region with several indicators: P(I_Ha06) = 140.6 days, P(I_Ha16) = 139.7 days, P(I_NaI) = 139.0 days, P(F W H M_CCF,HARPS) = 144.8 days, and P(F W H M_CCF,NIRPS) = 137.4 days. Based on the period detection with the activity indices, we chose a broad range of period space to use in our rotation period detection with the NIRPS line parameters. Our detected periods range between 137.4 and 144.8 days. To ensure that our search range includes all periods reported in the literature as well as our own detections, with some margin for uncertainty, we decided to limit our range to 130–150 days This range will be used from now on to identify a rotation period detection for this star.

4.3 Proxima rotation

Our literature compilation in Sect. 2 indicates the rotation period of Proxima is around 80–90 days. We detect significant periods near this range for several indices: P(I_Ha06) = 77.4 days, P(I_Ha16) = 77.5 days, P(I_NaI) = 98.9 days, P(F W H M_CCF,NIRPS) = 79.0 days, P(Cont_CCF,NIRPS) = 78.5 days, and P(BIS_CCF,NIRPS) = 89.7 days. There are no detections with HARPS CCF parameters, possibly due to the low S/N for the HARPS data of this star, as shown in Table 2. Almost all the indicators that detect a period close to the literature rotation period also show a second peak near 100 days, which in the case of Na I is the highest peak. Although understanding the nature of these two signals is beyond the scope of this study, they can be produced by active regions transiting the stellar disk at different latitudes (see Suárez Mascareño et al. 2025). The periods detected by our activity indicators range between 77.4 and 98.9 days. Similarly to the case of Gl 581, we chose a period range capable of including our detected and the reported periods, resulting a range of 70–110 days. This is the range we use to search for rotational modulation in the line profile parameters.

5 Methodology

Spectral lines were first identified in a high-S/N master spectrum created for each star. The master spectra was produced as explained in Appendix B. To identify spectral lines we employed the method of Cretignier et al. (2020) as discussed in Appendix C.1. This method also provides a central bandpass for each identified line that we used to integrate pEW and fit the line profiles, as discussed in this section.

To assess line variability, we extracted several parameters of the line profiles. We extracted parameters using two independent methods to minimise biased results: estimating the line pseudo-equivalent width (pEW) and modelling the line profiles to extract the model parameters. In the next subsections these two methods are explained.

5.1 Pseudo-equivalent widths

We calculated the equivalent width of each identified line using two different approaches: (a) a normalisation dependent method, where the pseudo-continuum is assumed to be unity, which we denominate as pEW_N and (b) a normalisation independent method where the pseudo-continuum is estimated at the 0.9 quantile of the flux inside the 2 Δλ_C20 (the bandpass extracted via the Cretignier et al. (2020) method; see Appendix C.1) region centred at the line core, which we designate pEW_Q. Furthermore, and because we do not have an idea of the central bandpass that maximises the activity information that can be extracted from each line, which should be line dependent, we used several central bandpass sizes for each pEW variant. We calculated five central bandpass sizes by multiplying the bandpass, Δλ_C20, by 2, 1, 1/2, 1/3, and 1/4. Thus, we calculated ten pEW variants, pEW_Nk and pEW_Qk with k = 2, 1, 1/2, 1/3, 1/4 using central bandpasses of 2 Δλ_C20, Δλ_C20, 1/2 Δλ_C20, 1/3 Δλ_C20, and 1/4 Δλ_C20 for the normalisation dependent and normalisation independent cases, respectively. This methodology will help us avoid not detecting activity variability due to normalisation issues and poor bandpass selection, as well as providing information on the optimal bandpass to measure such variability.

The pseudo-equivalent width of line with centre λ for discrete data is given by $$pEW=\sum _i^N(1-\frac{F\left({\lambda }_i\right)}{F_c})\delta {\lambda }_i,$$(1)

where N is the number of i points (pixels) inside the bandpass Δλ, F(λ_i) is the flux of pixel i inside the bandpass, F_c the flux in the pseudo-continuum, and δλ_i the size of pixel i, with $\mathrm{\Delta }\lambda ={\sum }_i^N\delta {\lambda }_i$. During the computation of all variants of pEW we performed a cubic interpolation inside the 2 Δλ_C20 region, i.e. the maximum bandpass size considered for each line. This interpolation was carried out to take into account the finite pixel size, i.e. to align the bandpass edges with the wavelength grid and force the integrated regions inside the bandpass to remain constant for all observations. A grid step similar to the wavelength step size was used to maintain the wavelength resolution. The uncertainties in pEW were estimated via error propagation in the flux.

5.2 Line fitting

An advantage of fitting the spectral lines over the use of methods such as equivalent widths or core-flux measurements, is that we gain access to line profile parameters such as line centre, width and depth, which provide specific information about which components of the line are being affected by variability. Deformations of spectral lines by activity phenomena will affect the parameters of the line profile.

We started by fitting all identified lines in the master spectrum³ using three different functions: Gaussian, Lorentzian, and Voigt. We used lmfit⁴ (Newville et al. 2024) to fit these models to each line, using the flux inside the Δλ_C20 central bandpass⁵. The best model was selected as the one with the lowest ${\chi }_{\text{red }}^2$ value. The master spectrum fit parameters of the best model, line centre, depth and FWHM were extracted to serve as reference line profile parameters. These models were then used to fit the observation spectra, and extract time series of the profile parameters for all identified lines.

Before fitting the models, we applied a cubic interpolation to align the pixels with the central bandpass limits. A wavelength grid with a step similar to the wavelength step of the spectrum was used to maintain the resolution.

5.3 Correlations with activity

For each of the line profile parameters and pEW time series obtained from the observations spectra, we calculated the Pearson correlation coefficient, ρ, with the three activity indicators F W H M_CCF,NIRPS, Hα06, and Na I. Significant correlations were chosen as the ones with |ρ| ≥ 0.4 and p–value ≤ 0.01%. This p-value was selected to minimise the number of correlations that can arise ‘by chance’ from calculating hundreds of thousands of correlations. For example, for Proxima, we identified 2767 lines, have 13 line parameters and three activity indicators, meaning we calculated 2767 × 13 × 3 = 107 913 correlation coefficients. For the chosen p-value, this means we could have 107 913 × 0.0001 ≈ 11 correlations arising ‘by chance’ (for Gl 581 this value is lower, since we have a lower number of identified lines).

As mentioned in Sect. 4, we detected a strong correlation between one of our ‘anchor’ activity indices, F W H M_CCF,NIRPS, and the NIRPS BERV for Gl 581. Correlations between line profile parameters and BERV can arise from telluric contamination or correction residuals affecting the profiles. In an effort to separate the correlations between the line parameters and activity indicators from those with BERV, we added another criteria for the selection of significant correlations: the correlations between a parameter and activity indices must always be higher that the correlations between the parameter and BERV, |ρ| > |ρ_BERV|.

5.4 Rotation period search

Our main interest in this work is to find activity-sensitive lines from an exoplanet detection and characterisation point of view. Stellar activity generally affects RVs at timescales connected to the rotation period of the star and its harmonics (Boisse et al. 2011), and the magnetic cycles periods with timescales that range from years to decades (e.g. Baliunas et al. 1995; Lovis et al. 2011; Gomes da Silva et al. 2012). Usually, indicators of activity sensitive to rotation modulation are also able to detect the long-term cycles (e.g. Gomes da Silva et al. 2022), since these cycles are characterised mainly by an increase and decrease (along with latitude migration) of the active regions to which the indices are sensitive. Therefore, by detecting lines sensitive to rotation modulation we are assuming they will also be sensitive to magnetic cycle variations, as is the case of the indices based on the Ca II H&K, Hα, and Na I lines (Gomes da Silva et al. 2011, 2012, 2014, 2022).

Another motivation for the period search is that phase shifts between activity indicators can lead to a degradation or a lack of correlations between them (Collier Cameron et al. 2019) even though the two signals could be following activity variability and be able to detect the stellar rotation. To account for this, we searched for periodicities in the line parameters that could be close to the star’s rotation period (see Sect. 2). We used the GLS periodogram (Zechmeister & Kürster 2009; see our Sect. 4) to search for periodicities in the line parameters time series close to the rotation periods detected with the known indices. Significant periods were selected as having a GLS power higher than the 0.1% FAP threshold, computed analytically following Baluev (2008). We do not report uncertainties on the detected GLS periods because period errors from GLS is a debatable issue (see VanderPlas 2018). The period search range for each star was the one described in Sect. 4 and the results are presented in Sect. 6.4.

6 Results

6.1 Cleaned line list and profile fit accuracy

In this section, we discuss the line list resulted from the fitting of observations, cleaning of outliers and accuracy of the fitting process, and median distribution of the line profile parameters.

We were able to fit the observations for 2297 lines for Gl 581 and 3032 lines for Proxima. Some extracted parameters included very large median errors and dispersions for a few lines that were several orders of magnitude higher than the median values; they resulted from poor fitting, and we considered them to be outliers. We tested using σ and median absolute deviation (MAD) clipping to clean the data using the median uncertainties and found that applying a 20–MAD clipping from the median was the best solution to remove the outliers for all parameters without the risk of removing too many lines with errors closer to the median of the distribution. As a result, 63 lines from the Gl 581 and 44 from the Proxima line lists were removed. We also performed a similar MAD clip to the median standard deviations, removing a further 127 lines for Gl 581 and 109 for Proxima.

We compared the median values of the time series from the observations fitting with those obtained from fitting the master spectrum to assess the accuracy of the time series (see Appendix D). Parameters that deviate significantly from the master fit values indicate disagreement between the profiles in the master and observations spectra (possibly caused by normalisation) or difficulties with the fitting process for the lower-S/N observation spectra. We are accurate to ~200 m/s in line centre position, ~750 m/s in FWHM, and ~2% in line depth, at the 1σ level. All lines with a difference between observation and master fit values higher than ± 3σ from the median were removed, affecting 86 lines for Gl 581 and 112 for proxima.

Finally, our cleaned line list contains 2021 lines for Gl 581 and 2767 for Proxima, with 1358 lines common to both line lists⁶. These are the lines whose parameters we use to compare with the activity indicators and search for rotational modulation.

Fig. 1 Cumulative distributions of the median uncertainties and standard deviation for extracted fit parameters and pEW for Gl 581 (top panels) and Proxima (bottom panels). The red lines are the distribution of uncertainties and the blue lines the distribution of standard deviations. The median value, standard deviation, and error are marked for each parameter.

6.2 Precision and dispersion of line profile parameters

The distributions of median uncertainties and standard deviations for the extracted fit parameters and pEW are provided in Fig. 1. Only two of the 10 different pEW variants, one with the wider bandpass (pEW_2N) and the other with the narrower bandpass (pEW_1/4N), are shown since the distribution of errors and dispersion for the different pEW are very similar. Overall, we are capable to detect line variability down to the level of the uncertainties of our fitted parameters. We are able to measure line positions with a median precision of ~130 m/s and median dispersions ~150 m/s. Our FWHM median uncertainties are at the level of ~400 m/s, with median dispersion of ~600 m/s. The median precision on the line depths is 0.005, and we measure median dispersions of 0.007. The minimum dispersion we could measure on individual lines are 27 m/s in line position, 105 m/s in FWHM, and 0.3% in line depth. The best precision we could achieve for our lines are 24 m/s in line centre, 71 m/s in FWHM, and 0.2% in line depth. This line centre precision is around half of the 50 m/s precision reported by Artigau et al. (2022) for the RVs of their best lines, achieved using the LBL technique on Barnard’s star with Spectrograph for Planetary Research and Observations (SPiRou), and around two times the precision reported by Dumusque (2018) on some spectral lines of α Cen B with HARPS (visual wavelengths) using also the LBL method. For pEW, we obtain median errors between 0.5 and 2 mÅ and median dispersion between 1 and 6 mÅ. This showcases our ability to detect low-amplitude variability in the line profiles of NIR spectra with NIRPS.

6.3 Correlations with activity indices

In this section, we compare the number of lines with significant correlations between the line profile parameters, the pEW parameters, and the activity indices.

We found significant correlations between at least one line parameter and one activity indicator for 312 lines (15%) for Gl 581 and 760 lines (27%) for Proxima. These lines include 78 lines in common between the two stars, which shows that most of the lines with correlations are specific to each star and cannot be used to follow activity, as measured by our indicators, for both stars.

The number of lines with detected significant correlations between our line parameters, line centre, FWHM, depth, and pEW, discriminated by the three activity indicators is provided in Fig. 2. The parameter with more detections of significant correlations with activity for both stars is the pEW. This could be a result of the difficulty of fitting and measuring precise small variations in the line profile parameters for lines in the NIR with the NIRPS resolution, S/N, and/or telluric residuals present in these spectra. But we should also note that the pEW parameter includes 10 different versions that can best be adapted for each line, thus increasing the number of detections.

We can also confirm that all line parameters have higher number correlations with the NIRPS CCF FWHM than the chromospheric activity indices, for both stars. This can be a result of the different quality of the spectra between NIRPS and HARPS, from where the chromospheric indices were computed, with the S/N of the NIRPS observations being always higher than those of HARPS. This could also mean that we are mostly detecting line variability that is not caused by chromospheric emission in the lines. Since we do not expect that the inhibition of convective blueshift has a strong effect on the distortion of lines of M dwarfs (Liebing et al. 2021), we could be measuring the flux effect caused by active regions with different temperature contrasts to the quiet photosphere transiting the stellar disk of both stars. But we could also be measuring the Zeeman splitting of the line profiles caused by the presence of magnetic fields in the active regions (e.g. Fuhrmeister et al. 2022). For the two stars, the Hα06 chromospheric indicator always performs slightly better than Na I, at having correlations with line profile parameters.

In Fig. 3, we show the distribution of the lines with significant correlation between at least one parameter and one activity index as a function of line centre, FWHM, and depth, for Gl 581 and Proxima. Lines with significant correlation in common between the two stars are also shown. From the distributions of FWHM and line depth, we can observe that Proxima has a larger number of wider and deeper lines than Gl 581. In terms of wavelength position, the lines with correlations are in general distributed along all the NIRPS spectral range, although for Gl 581 there are more of these lines in the H band, redder than 14 000 Å, while for Proxima, most of the lines are located in the Y band, between around 10 000 and 12 000 Å. The lines with significant correlations detected on both stars are also well distributed along the full spectral range. The distribution of lines with correlations along the FWHM is mostly concentrated in the range around 4–23 km/s, with the median of the distribution at 6.8 km/s for Gl 581 and 8.3 km/s for Proxima. The lines with correlations in common between the two stars have a tendency to be slightly wider for Proxima, with the median of the distribution at 9.2 km/s, while for Gl 581, the median of the distribution stands at 6.8 km/s. In relation to line depth, in general, the lines of Proxima are deeper than those of Gl 581, and the lines with significant correlation follow this trend, with the median of the distribution of lines with correlations located at a depth of 0.15 for Gl 581 and 0.23 for Proxima. For Gl 581, the lines with correlations in common for the two stars have a median depth of 0.19, while for Proxima they have a slightly higher median depth of 0.29. Both stars show increasing fraction of lines with correlations for increasing line depth. For Gl 581, the fraction of lines with correlations is almost constant for shallow lines but increases significantly for depths higher than 0.4, reaching around 80%, while in the case of Proxima there is a steady increase with line depth. Since most of the lines with detected significant correlations are shallow lines with depths in the range 0.05–0.2, we are confident we are able to measure significant line variability in shallow lines, and that the increasing fraction of correlated lines with line depth is not due to a methodological bias.

Fig. 2 Number of lines with significant correlations between line profile parameters, pEW, and the activity indices for Gl 581 (top panel) and Proxima (bottom panel).

Fig. 3 Distribution of lines with significant correlations (|ρ| ≥ 0.4, |ρ| > |ρBERV|, and p-value ≤0.01%) with activity in at least one parameter as a function of line profile parameters. The black histogram includes all the analysed lines, blue represents the lines with significant correlations between at least one parameter and one activity index, and red shows the lines with significant correlations common to both stars. The purple line indicates the fraction of lines with correlations in at least one profile parameter and at least one activity indicator to the full line list; the errorbars are given by Poisson statistics.

6.4 Rotation period detections

As discussed in Sect. 5.4, we carried out a rotation period search using all extracted parameters time series. Overall, we have 36 lines with significant P_rot detection in at least one parameter for Gl 581 (out of 2021, 1.7%) and 874 lines (out of 2767, 32%) for Proxima. There are 17 of these lines that show profile variability consistent with the rotation period of both stars, which can be used as an activity proxy. These lines are further discussed in Sect. 6.5.

The distribution of median line centre, FWHM, and depth for the lines with detected rotation period, for both stars, are provided in Fig. 4. The wavelength distribution of lines with detected P_rot is mostly located in the J band for Gl 581 and Y band for Proxima, with the lines with detections common to the two stars, located mainly in a region between 12 000–13 500 Å. The lines with rotation detection have FWHM in the region between 4 and 20 km/s for Gl 581, with the lines of Proxima being more spread in FWHM, between 3 and 33 km/s. The median values of the FWHM distribution for the lines in common for Gl 581 and Proxima are 8.8 and 13.6 km/s, respectively. In terms of line depth, for Gl 581, the lines with P_rot detection appear distributed in two regions, between 0.05 and 0.2 and deeper than 0.4, with the median of the distribution at 0.24. However, there is a jump in the fraction of lines with detections between 0.4 and 0.6, reaching values of 75% at around 0.55. In the case of Proxima, the median of the depth distribution is located at 0.21, and the fraction of lines with a detection increases steadily with line depth, reaching 100% for depths higher than 0.6, similar to the case of the correlation analysis.

The most striking differences between the lines with significant correlations with activity and the lines with rotation detection is that for Gl 581 there are only 36 lines with rotation modulation but 312 with correlations with activity, while for Proxima, we found more lines with detected rotation period than lines with significant correlations with activity.

In the case of Gl 581, our first guess is that, for a line to have detectable rotation modulation on a given parameter, the effects of activity have to have a higher amplitude, than for the line parameters to have correlation with activity, which can vary between low and high correlation coefficients, depending on the effect amplitude on the parameter. We compared the distribution of fitted parameters standard deviations of line centre, FWHM and line depth for the lines with correlations and the lines with rotation detection. Contrary to our expectations, the lines with rotation detection are all mainly concentrated at the lowest dispersion range for all parameters (σ (centre) = 0.04–0.5 km/s; σFWHM = 0.18–1 km/s; σ(depth) = 0.003–0.01), while the lines with correlations with activity range between the full range of the dispersion distribution. This implies that, having a higher dispersion on the parameters does not ensure that it will be easier to detect the rotation period of the star.

For the case of Proxima, the discrepancy between lines with correlations and lines with rotation detection comes from our thresholds. Our enforcement that selected correlations must have higher absolute correlation coefficients than the absolute correlation coefficient with BERV is not used for the period analysis, since we can discriminate between the rotation period and periodicities caused by BERV, namely the harmonics of the 1-year signal. This causes some lines with high correlations with activity not to be selected due to BERV correlation. Furthermore, we have a very low p-value threshold of 0.01%, which decreases the number of selected lines with correlations. This is also why we have five lines with detected rotation period for Gl 581 and 346 for Proxima that do not pass the correlation selection. The line lists of selected activity-sensitive lines are very sensitive to the imposed selection criteria and different selections can be used, depending on the aims.

Nevertheless, the results described here show that several line profile parameters can be used to measure the rotation period of M dwarfs in the NIR and that, for more active late-type M dwarfs such as Proxima, several hundred lines have their profiles affected by rotational modulation and can be used to infer rotation periods.

Fig. 4 Distribution of lines with significant Prot detection (FAP ≤ 0.1%) as a function of line profile parameters. The colours and symbols are the same as in Fig. 3.

6.5 Activity-sensitive line lists

The aim of this work is twofold: to measure the effect of stellar activity on the NIR spectra of M dwarfs and to find new activity-sensitive lines in this wavelength domain. We identified 2021 lines for Gl 581 and 2767 for Proxima, which we used to analyse their profile parameters in search for activity effects.

The correlation analysis in Sect. 6.3 provides information about the correlation between different line profile parameters and the three activity indicators F W H M_CCF,NIRPS, I_Hα006, and I_NaI. The analysis in Sect. 6.4 provides information about the periodicities found in the profile parameters. We deliver these line lists in two CDS tables, one for each star. A description of columns is provided in Table E.1. These tables include 312 lines for Gl 581 and 760 lines for Proxima with significant correlations (|ρ| ≥ 0.4, |ρ| > |ρ_BERV|, and p-value ≤0.01%) between at least one line parameter and at least one activity index, and 36 lines for Gl 581 and 874 lines for Proxima with significant rotation period detection (130 ≤ P_Gl581 ≤ 150 days, 70 ≤ P_Proxima ≤ 110 days, FAP ≤0.1%). These line lists could be useful to select lines for RV extraction that are not too much affected by activity variability for these two stars.

There are 17 lines with rotation detection, common to both stars, whose basic parameters are provided in Table F.1. To identify the species responsible for each line, we cross-matched the line centres with the VALD3 database (Ryabchikova et al. 2015), selected for the wavelength range of NIRPS and an M-dwarf with the characteristics of Proxima. In this last table, there is a line previously analysed in the literature: K I λ12435.644 (Fuhrmeister et al. 2022; Terrien et al. 2022; Cortés-Zuleta et al. 2023). Although, in a study involving 324 M dwarfs, Fuhrmeister et al. (2022) found that the K I line rarely correlates with the strength of the Hα line, Terrien et al. (2022) found that for Gl 699 (M4V) and Teegarden’s star (M7V), the equivalent width and width of this line is modulated by the rotation of both stars, and that Zeeman splitting could be the cause of such modulation. In the analysis of Gl205 (M1.5V) with SPiRou, Cortés-Zuleta et al. (2023) did found the rotation period signature in the pEW of K I, however without enough strength to be considered significant. Our analysis for Gl 581 (M3V) and Proxima (M5V) indicate that both the line’s FWHM and depth are modulated by rotation.

Our 17 lines with rotational modulation for both stars all have stronger detections for Proxima, as demonstrated by the FAP values in Table F.1. The lines that most significantly detect the rotation of Gl 581 are KI λ12435.644, Mn I λ12903.288, K I λ11772.858, and Ti I λ9790.368, and they should be considered as potential activity indicators in the NIR for M-dwarf stars.

7 Conclusions

The aim of this study was to analyse the effects of stellar activity in the NIR, and to search for activity-sensitive lines in this wavelength region. For the first time we identified and analysed for activity effects thousands of NIR lines for two M dwarf stars of the NIRPS GTO with different activity levels and stellar types, Gl 581 (M3V) and Proxima (M5.5V).

We achieved a median precision on our line profile parameters of ~120 m/s in line centre position, ~400 m/s in FWHM, and ~0.5% in line depth, at the 1σ level. We measured dispersions in the parameters down to the precision level. We found hundreds of lines with strong correlations with known activity indicators for both stars. For Gl 581 15% of the analysed lines showed strong correlations (|ρ| ≥ 0.4, |ρ| > |ρ_BERV|, and FAP ≤ 0.01%) with the activity indices, while for Proxima we found a higher number of 27% lines. These lines can be found along the Y, J, and H bands, with a tendency for more lines with correlations in the H band for Gl 581, while for Proxima the tendency is for the majority of the lines to be located in the Y and J bands. The fraction of lines with correlations tends to increase with line depth for both stars. Most of the lines with correlations were detected using pEW, although we note that we used ten different versions of this parameter, which can be better adjusted to each line profile, thus contributing to maximising the number of detections.

We also found hundreds of lines for Proxima capable of detecting the star rotation period (32%), but only 36 lines for Gl 581 (1.7%). This discrepancy might be a result of the different activity levels of the two stars (Gl 581 is less active, with no Hα emission), but also their internal structure (Proxima with a fully convective interior).

Only a fraction of lines with correlations with activity (78 lines) or rotation period detection (17 lines) are common to both stars, further indicating the spectra are being affected differently for these stars and suggesting that the activity influence on the spectrum of M-dwarf stars should be analysed on a case-by-case basis.

Other results from this work include the following:

• Both versions of the Hα index have the same behaviour for these M dwarfs, regardless of activity level, spectral type, and internal structure. In the light of the results of Gomes da Silva et al. (2022), we recommend the use of Hα06 for the analysis of FGK and M dwarfs in order to have homogeneous activity indices, which can be used to compare activity between different spectral types (after applying a bolometric correction);

• Although the NIRPS CCF FWHM, Hα06, and Na I indices are very sensitive to rotation modulation, the first shows a higher number of lines with correlations with line profile parameters and detected rotation periods than the spectral indices based on the Hα and Na I lines. This is probably due to the higher S/N in the NIRPS observations when compared to HARPS.

This work shows that stellar activity affects hundreds of spectral lines in the NIR, mainly at higher activity levels, as is the case of Proxima. We expect that the activity phenomena distorting the lines arise from the flux effect caused by the different temperature contrasts between active regions and the quiet photosphere, or Zeeman splitting of the line profiles (e.g. Fuhrmeister et al. 2022), although we cannot rule out line-core emission from the chromosphere affecting some lines. This poses a challenge for exoplanet detection and characterisation using the RV method in the NIR. Methods for correcting line variability should be developed at the spectral level to mitigate these effects and increase RV precision.

The line lists with the correlation and period detection analysis are provided in CDS for the community. We curated a list of 17 spectral lines with rotation detection for both stars that can used as potential NIR activity proxies for M dwarfs. These lists can be used to further develop techniques to analyse and mitigate the effects of stellar activity in the NIR spectra of M dwarfs. Furthermore, these activity-sensitive lines can be used to reduce the effects of stellar activity on NIR RV by excluding them from CCF masks, LBL analysis, and template matching methods. Additionally, they can also be used to detect rotation periods and identify false-positive exoplanet detections.

Data availability

The tables described in Table E.1 are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/700/A177

Acknowledgements

This work has made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna. This publication makes use of The Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva (CH) dedicated to extrasolar planets data visualisation, exchange and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. The DACE platform is available at https://dace.unige.ch. JGdS, ED-M, NCS, SCB, ARCS & EC acknowledge the support from FCT – Fundação para a Ciência e a Tecnologia through national funds by these grants: UIDB/04434/2020, UIDP/04434/2020. JGdS and ED-M acknowledge the support from the Fundação para a Ciência e a Tecnologia (FCT) through national funds by the grant 2022.04416.PTDC. ED-M further acknowledges the support from FCT through Stimulus FCT contract 2021.01294.CEECIND. ED-M acknowledges the support by the Ramón y Cajal contract RyC2022-035854-I funded by MICIU/AEI/10.13039/501100011033 and by ESF+. Co-funded by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. ASM, NN, JIGH, FGT, VMP, JLR, RR & AKS acknowledge financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) projects PID2020-117493GB-I00 and PID2023-149982NB-I00. XDe, XB, ACar, TF & VY acknowledge funding from the French ANR under contract number ANR18-CE310019 (SPlaSH), and 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 Grenoble-Alpes University. ÉA, RA, FBa, LB, BB, AB, CC, NJC, AD-B, LD, RD, AL, PLam, OL, LMa, LMo, JSA, PV, TV & JPW acknowledge the financial support of the FRQ-NT through the Centre de recherche en astrophysique du Québec as well as the support from the Trottier Family Foundation and the Trottier Institute for Research on Exoplanets. ÉA, TA, FBa, RD, LMa, J-SM, MO, JS-A & PV acknowledges support from Canada Foundation for Innovation (CFI) program, the Université de Montréal and Université Laval, the Canada Economic Development (CED) program and the Ministere of Economy, Innovation and Energy (MEIE). NN acknowledges financial support by Light Bridges S.L, Las Palmas de Gran Canaria. NN acknowledges funding from Light Bridges for the Doctoral Thesis “Habitable Earth-like planets with ESPRESSO and NIRPS”, in cooperation with the Instituto de Astrofísica de Canarias, and the use of Indefeasible Computer Rights (ICR) being commissioned at the ASTRO POC project in the Island of Tenerife, Canary Islands (Spain). The ICR-ASTRONOMY used for his research was provided by Light Bridges in cooperation with Hewlett Packard Enterprise (HPE). The Board of Observational and Instrumental Astronomy (NAOS) at the Federal University of Rio Grande do Norte’s research activities are supported by continuous grants from the Brazilian funding agency CNPq. This study was partially funded by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior–Brasil (CAPES) – Finance Code 001 and the CAPES-Print program. KAM acknowledges support from the Swiss National Science Foundation (SNSF) under the Postdoc Mobility grant P500PT_230225. RA acknowledges the Swiss National Science Foundation (SNSF) support under the Post-Doc Mobility grant P500PT_222212 and the support of the Institut Trottier de Recherche sur les Exoplanètes (IREx). SCB acknowledges the support from Fundação para a Ciência e Tecnologia (FCT) in the form of a work contract through the Scientific Employment Incentive program with reference 2023.06687.CEECIND. LB acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project SPICE DUNE, grant agreement No 947634). This material reflects only the authors’ views and the Commission is not liable for any use that may be made of the information contained therein. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. BLCM & AMM acknowledge CAPES postdoctoral fellowships. BLCM acknowledges CNPq research fellowships (Grant No. 305804/2022-7). ARCS acknowledges the support from Fundaçao para a Ciência e a Tecnologia (FCT) through the fellowship 2021.07856.BD. NBC acknowledges support from an NSERC Discovery Grant, a Canada Research Chair, and an Arthur B. McDonald Fellowship, and thanks the Trottier Space Institute for its financial support and dynamic intellectual environment. LD acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) and from the Fonds de recherche du Québec (FRQ) – Secteur Nature et technologies. JRM acknowledges CNPq research fellowships (Grant No. 308928/2019-9). XDu acknowledges the support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement SCORE No 851555) and from the Swiss National Science Foundation under the grant SPECTRE (No 200021_215200). DE acknowledge support from the Swiss National Science Foundation for project 200021_200726. The authors acknowledge the financial support of the SNSF. FG acknowledges support from the Fonds de recherche du Québec (FRQ) – Secteur Nature et technologies under file #350366. H.J.H. acknowledges funding from eSSENCE (grant number eSSENCE@LU 9:3), the Swedish National Research Council (project number 2023-05307), The Crafoord foundation and the Royal Physiographic Society of Lund, through The Fund of the Walter Gyllenberg Foundation. AL acknowledges support from the Fonds de recherche du Québec (FRQ) – Secteur Nature et technologies under file #349961. ICL acknowledges CNPq research fellowships (Grant No. 313103/2022-4). LMo acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), [funding reference number 589653]. CMo acknowledges the funding from the Swiss National Science Foundation under grant 200021_204847 “PlanetsInTime”. CPi acknowledges support from the NSERC Vanier scholarship, and the Trottier Family Foundation. CPi also acknowledges support from the E. Margaret Burbidge Prize Postdoctoral Fellowship from the Brinson Foundation. AKS acknowledges financial support from La Caixa Foundation (ID 100010434) under the grant LCF/BQ/DI23/11990071. TV acknowledges support from the Fonds de recherche du Québec (FRQ) – Secteur Nature et technologies under file #320056.

Appendix A GLS plots of the known activity indicators

Fig. A.1 HARPS (top panels) and NIRPS (bottom panels) activity index GLS for Gl 581 (a) and Proxima (b). The vertical red band in the GLS marks the literature rotation period span between the minimum and the maximum rotation periods found in the literature. The period value in red represents the highest peak closest to the literature values. The horizontal dotted lines in GLS indicate the FAP = 0.1% level

Appendix B NIRPS master spectra

In this section we describe the computation of the NIRPS master spectra for Gl 581 and Proxima, which was used to identify the spectral lines to be analysed and to chose the best model, Gaussian, Lorentzian, or Voigt, to fit the individual observations line profiles.

B.1 Spectra quality selection

Before proceeding to the calculation of the master spectrum, we analysed the quality of the observations to include only high-S/N spectra and spectral orders. As can be observed in Table 2, the median S/N of the observations is significantly higher for NIRPS spectra than for HARPS. For HARPS, all observations with a median S/N below 3 were excluded from the study. This includes 12 observations for Proxima and one for Gl 581. Furthermore, the S/N is not constant for a spectrum along the wavelength range. Figure B.1 shows the S/N per extracted spectral order for all available spectra of Gl 581. In the regions between spectral orders 42 and 47 (approximately between 13500 Å and 14100 Å()) and orders higher than 67 (≳18000 Å) the S/N quickly drops to values close to zero. The 42–47 orders region corresponds to the unusable domain of the deep telluric water band between the J and H photometric bands. After a thorough visual inspection of the individual spectra, we decided to exclude these spectral orders from further analysis in this work.

Fig. B.1 S/N per order for all NIRPS observations of Gl 581. The superposed black line represents the median S/N of all spectra per order. The approximate positions of the NIR Y, J, and H photometric bandpasses are marked above the x-axis.

B.2 Master spectrum

To identify the spectral lines to be analysed, we first built a master spectrum, i.e. a high-S/N co-added spectrum comprising all observations, for each star. When the individual spectra of each observation were shifted to the stellar rest frame by subtracting the stellar systemic and Barycentric Earth RVs, their wavelength grids become slightly different, with different initial and final wavelength values. Therefore, we need to interpolate all spectra into a common wavelength grid to be able to co-add them. To construct the new grid, we determined each spectrum initial and final wavelength values together with the median wavelength step for each order and selected the maximum initial, minimum final, and minimum median step values as the values in the grid for each spectral order, a method also used by Dumusque (2018). We then interpolated each individual spectral order of each individual spectrum with cubic interpolation and co-added them using specutils⁷ (Earl et al. 2024), which conveniently takes care of flux errors during these processes.

During the co-adding process we truncated the spectrum extremities of each spectral order by 10 Å to avoid the drop in S/N at the edges due to the decrease in the blaze function flux, which will affect the normalisation process. After co-adding all spectra, the resulting master spectrum was normalised to unity. We used a two-step polynomial normalisation algorithm. First we performed a sequential σ-clip to approximate the spectral order continuum by eliminating all flux values below the median minus 2σ and fitted a fifth-order polynomial to the continuum. Then, we repeated the process but this time using a sequential σ-clip with 1.7σ and fitted the resulting continuum using a third-order polynomial. The resulting continuum polynomial was then interpolated into the original wavelength grid and divided from the original flux.

Appendix C Identification of spectral lines in the NIRPS spectra

C.1 Line identification using the Cretignier method

One of the important issues with line identification has to do with the definition of the line central bandpass. Some line identification methods use a fixed bandpass in wavelength; however, this is not a good practice since the wavelength step is not constant in wavelength for most unprocessed spectra. Fixing a bandpass value in RV space solves the wavelength step issue, but does not take into account the width of each line. We were thus interested in a method that could take into account each line width.

Cretignier et al. (2020) developed a method to identify and provide parameters such as line centre, depth and central bandpasses for symmetric spectral lines. The method is based on the first and second derivatives of the spectrum whose maxima and minima provide constraints for the parameters and for the symmetry assessment. Strong asymmetry in lines is assumed to be caused by line blending, and lines exhibiting asymmetry above a certain threshold are rejected by the selection criteria described in Cretignier et al. (2020, Appendix C). The central bandpasses we use to fit the line profiles and calculate the pEW are the ones returned by this method, which we denominate as Δλ_C20.

To identify the spectral lines, We applied this method to the master spectrum of both stars. Due to the overlapping of spectral orders our initial line list included duplicate lines at the orders extremities. We identified these cases by selecting lines whose centres were closer than the central bandpass of one of the lines, and removed the line closer to the spectral order edges. Since the S/N decreases towards the extremities due to the lower blaze function values in those regions, we selected lines from the duplicate list that were farther away from the wavelength initial and final values. Our initial line search resulted in a list containing 6210 lines for Gl 581 and 9396 for Proxima. After applying the selection criteria discussed in Cretignier et al. (2020, Appendix C)⁸ to this line list, we recovered a list of lines including 2391 for Gl 581 and 3125 for Proxima.

C.2 Master spectra fitting and resulting line list

To find the best profile model to fit to the observations and assess our fitting procedure, we fitted all lines we recovered with the Cretignier et al. method using the master spectrum as described in Sect. 5.2. The algorithm was incapable of properly fitting 12 of the Gl 581 lines and three lines of the Proxima list. Failure to fit identified lines can come from noise variability caused by low S/N for very shallow lines, residuals from telluric correction affecting the line profiles and/or blended lines that passed the selection criteria. This resulted in a list of 2379 lines for Gl 581 and 3122 lines for Proxima.

Some regions of the observations spectra have zero-valued flux counts that change from observation to observation. These zero flux values are due to intense telluric lines with almost zero flux. The pipeline forces them to zero to avoid spurious residuals. Before fitting the observations, we excluded all lines with zero fluxes in the central bandpasses in at least one observation, excluding 58 lines from the Gl 581 and 82 from the Proxima line lists. Consequently, the final set of lines to be fitted to the observations includes 2321 lines for Gl 581 and 3040 for Proxima.

Appendix D Accuracy of the fitted observations parameters

The accuracy of our observation fitting procedure was tested by comparing the median values of the time series from the observations fitting with those obtained from fitting the master spectrum to assess the accuracy of the time series. The results are provided in Fig. D.1. We are accurate to ~200 m/s in line centre position, ~750 m/s in FWHM, and ~2% in line depth, at the 1σ level. Lines with observation median parameters whose difference to the master values are above the 3σ levels from the median were removed. Most of these lines are located at the edges of the photometric bands Y, J and H, as can be observed in the line centre panels. This is due to the significant decrease in S/N at these regions as can also be observed in Fig. B.1. In terms of line depth, there is a tendency for the non-accurate lines to have lower depths as fitted from the observations than the master fit depths.

Fig. D.1 Difference between the median parameter values extracted from the fit to observations and the master spectra parameters. The horizontal lines mark the position of the 3σ levels around the median. The blue circles are lines with values below the 1σ levels; the red circles are the removed lines with values above the 3σ levels. Values of 1σ and 3σ levels are displayed for each parameter.

Appendix E Description of the CDS tables of the line lists of Gl 581 and Proxima

Appendix F Lines with significant rotation period detection for both stars

References

All Tables

All Figures

	Fig. 1 Cumulative distributions of the median uncertainties and standard deviation for extracted fit parameters and pEW for Gl 581 (top panels) and Proxima (bottom panels). The red lines are the distribution of uncertainties and the blue lines the distribution of standard deviations. The median value, standard deviation, and error are marked for each parameter.
In the text

	Fig. 2 Number of lines with significant correlations between line profile parameters, pEW, and the activity indices for Gl 581 (top panel) and Proxima (bottom panel).
In the text

Fig. 3 Distribution of lines with significant correlations (|ρ| ≥ 0.4, |ρ| > |ρBERV|, and p-value ≤0.01%) with activity in at least one parameter as a function of line profile parameters. The black histogram includes all the analysed lines, blue represents the lines with significant correlations between at least one parameter and one activity index, and red shows the lines with significant correlations common to both stars. The purple line indicates the fraction of lines with correlations in at least one profile parameter and at least one activity indicator to the full line list; the errorbars are given by Poisson statistics.

In the text

	Fig. 4 Distribution of lines with significant Prot detection (FAP ≤ 0.1%) as a function of line profile parameters. The colours and symbols are the same as in Fig. 3.
In the text

Fig. A.1 HARPS (top panels) and NIRPS (bottom panels) activity index GLS for Gl 581 (a) and Proxima (b). The vertical red band in the GLS marks the literature rotation period span between the minimum and the maximum rotation periods found in the literature. The period value in red represents the highest peak closest to the literature values. The horizontal dotted lines in GLS indicate the FAP = 0.1% level

In the text

	Fig. B.1 S/N per order for all NIRPS observations of Gl 581. The superposed black line represents the median S/N of all spectra per order. The approximate positions of the NIR Y, J, and H photometric bandpasses are marked above the x-axis.
In the text

	Fig. D.1 Difference between the median parameter values extracted from the fit to observations and the master spectra parameters. The horizontal lines mark the position of the 3σ levels around the median. The blue circles are lines with values below the 1σ levels; the red circles are the removed lines with values above the 3σ levels. Values of 1σ and 3σ levels are displayed for each parameter.
In the text