© The Authors 2025

1 Introduction

Since the discovery of the first exoplanet around a solar-like star, 51 Peg b (Mayor & Queloz 1995), the development of increasingly precise and optimized instruments to measure planetary mass (e.g., HARPS¹, Mayor et al. 2003; HARPS-N², Cosentino et al. 2012; ESPRESSO³, Pepe et al. 2021; EXPRES⁴, Jurgenson et al. 2016; NEID⁵, Schwab et al. 2016) and radius (e.g., CoRoT⁶, Auvergne et al. 2009; Kepler, Borucki et al. 2010; TESS⁷, Ricker et al. 2015) has facilitated the detection and characterization of a growing number of exoplanets. In particular, studying their atmospheres provides insight into the physical properties of these planets, namely their chemical composition and surface temperature. Hot exoplanets, such as hot Jupiters, are ideal candidates for atmospheric studies given that their H and He-dominated atmospheres are significantly extended and have large scale heights. The proximity of these planets to their host stars results in intense heating of their thermobase by stellar radiation, causing hydrodynamic expansion and, in some cases, loss of gas from the upper atmospheric layers (Lammer et al. 2003; Vidal-Madjar et al. 2003). Atmospheric mass loss not only provides insights into the physical properties of exoplanets but can also help explain some of the puzzling observations seen in population studies. For instance, the lack of short-period Neptune-mass planets, also known as the hot Neptune desert, could be due to the fact that some planets are not massive enough to retain their expanding gaseous atmosphere (e.g., Lecavelier Des Etangs 2007; Davis & Wheatley 2009; Beaugé & Nesvorný 2013; Lundkvist et al. 2016; Mazeh et al. 2016; Owen 2019). Additionally, the clear gap in the radius distribution between 1.5 R_⊕ and 2 R_⊕, otherwise known as the radius valley, may be caused by some exoplanets retaining a substantial atmosphere while others losing it through photoevaporation (e.g., Lopez et al. 2012; Lopez & Fortney 2013; Owen & Wu 2013, 2017; Fulton et al. 2017; Rogers & Owen 2021) or core-powered mass loss (e.g., Ginzburg et al. 2016; Gupta & Schlichting 2019, 2020, 2021). Therefore, developing techniques to study the impact of various atmospheric escape mechanisms is crucial to improving our understanding of the evolution and demographics of exoplanets.

Historically, studying the H Ly-α line through transmission spectroscopy has been a prominent method to explore atmospheric escape. It has yielded significant results, especially for hot Jupiters orbiting bright stars like HD 209458b (Vidal-Madjar et al. 2003, 2008; Ehrenreich et al. 2008) and HD 189733b (Lecavelier Des Etangs et al. 2010; Lecavelier des Etangs et al. 2012; Bourrier et al. 2013, 2020). However, as H Ly-α is in the ultraviolet range, it is not observable from the ground and space-based observations are plagued with interstellar medium (ISM) absorption and geocoronal emission (Landsman & Simon 1993; Wood et al. 2002, 2005). Although these difficulties can be overcome by restricting the spectral analysis to the wings of H Ly-α, ISM absorption will in any case limit observations to extrasolar systems in close proximity to the Sun.

Another possible tracer of escaping atmospheres, initially theorized by Seager & Sasselov (2000), is the metastable HeI (2³S) triplet. This transition presents features at 10 832.1, 10 833.2, and 10 833.3 Å, a wavelength region devoid of ISM and geocoronal absorption (Indriolo et al. 2009). Furthermore, as these lines are in the near-infrared (IR) they can be studied with dedicated high-resolution ground-based instruments, such as HPF⁸ (Mahadevan et al. 2012, 2014), SPIRou⁹ (Donati et al. 2020), CARMENES¹⁰ (Quirrenbach et al. 2014), NIRPS (Bouchy et al. 2017), or NIGHT¹¹, which is currently in development (Farret Jentink et al. 2024). The He triplet has successfully been detected in several exoplanets, including, but not limited to, WASP-107b (Spake et al. 2018; Allart et al. 2019; Kirk et al. 2020; Spake et al. 2021), HAT-P-11b (Allart et al. 2018; Mansfield et al. 2018), HD 189733b (Salz et al. 2018; Guilluy et al. 2020), and most recently WASP-69b (Nortmann et al. 2018; Vissapragada et al. 2020, 2022). It is known that the He triplet’s core is formed at the chromosphere level and therefore is sensitive to magnetic fields (Guilluy et al. 2020). Magnetic fields are not supposed to evolve significantly over the timescale of a transit and it is therefore assumed that the He triplet is stable over such a timescale. However, the He triplet could also be sensitive to other stellar noises such as oscillations, granulation, and super-granulation that evolve on timescales of a few minutes to a few days, or be contaminated by telluric lines that can also evolve significantly over a few hours. In this paper, we investigate the variability of the solar He triplet over a few hours for several days and how the observed variability impacts the signature of an exoplanet atmosphere. Additionally, we search for links between this variability and stellar and telluric contamination tracers to better understand the origin of this helium transition’s variations. To this end, the He triplet’s wavelength range offers a unique advantage as it includes both a telluric water line and a stable photospheric SiI line. This positioning allows us to effectively separate telluric and stellar contributions and study them independently.

We provide a brief outline of the NIRPS instrument and the observations used in Section 2. We describe the methods used to study variations in the He triplet properties, their possible causes, and their effects on retrieved planetary parameters in Sect. 3. Finally, we present and discuss the implications of our results in Sect. 4.

2 NIRPS observations

NIRPS is a high-resolution, near-IR spectrograph installed on the 3.6 m telescope at the European Southern Observatory’s (ESO) La Silla site in Chile, along with the HARPS spectrograph. One of the goals of this instrument is to study the atmospheres and orbital architectures of a variety of exoplanets using high-resolution time-series spectroscopy, in both transmission and emission (Artigau et al. 2024; Bouchy et al. in prep.). Such studies could improve our understanding of the mechanisms driving exoplanet formation and evolution, as well as observed trends in exoplanet demographics. NIRPS complements the spectral range of HARPS, by covering the Y, J, and H bands (0.97–1.81 μm), and offers two modes of observation depending on the magnitude of the observed star: a high-accuracy mode (HA) and a high-efficiency mode (HE) that can reach spectral resolutions of 90 000 and 75 000, respectively. The HE mode, which uses a larger science fiber, is less affected by modal noise and thus provides more precise retrievals of atomic and molecular features from atmospheric spectra. We therefore only analyze the HE mode solar data in this paper.

The HARPS experiment for light integrated over the Sun (HELIOS) is installed outside the 3.6 m telescope building. It consists of a small robotic telescope that focuses the disk-integrated Sun’s light onto an integrating sphere and feeds it to HARPS and NIRPS via optical fibers. This telescope therefore allows for the acquisition of high-precision solar spectra for several hours per day, on any possible day. As we aim to study the variability of the He triplet over different timescales, we extracted the NIRPS spectra for multiple days of HELIOS observations. Details concerning these observations can be found in Table A.1.

Further insights into the design and performance of NIRPS and HELIOS are provided by Wildi et al. (2017); Blind et al. (2022); Cointepas et al. (2022); Frensch et al. (2022); Thibault et al. (2022); Wildi et al. (2022); Al Moulla et al. (2023); Bouchy et al. (2024) and the dedicated ESO webpages¹²,¹³.

3 Methodology

The following sections describe the methods used to characterize the He triplet variability and study its impact on retrieved planetary parameters. As multiple analyses are done on different portions of the data, we clarify what is done on what in the following points:

1. The long-term variability analysis (Sects. 3.2.5 and 4.2, and Fig. 5) was carried out with all available observation days.

2. The short-term variability analysis (Sects. 3.2.4 and 4.1, Appendices B and E, and Figs. 2, 3, and 4) was carried out with all days of observation. All days presented similar results, but we opted to present only the results for January 26, 2023 and January 27, 2023 as these two provided the most adequate representation of the trends witnessed across all days.

3. The injection-recovery tests (Sects. 3.3 and 4.3, Appendix C, and Figs. 6, 7, and 8) were carried out with all days of observation. We present the results for January 26, 2023 and report the results of all other analyses in Table C.1.

3.1 Data extraction and reduction

The raw NIRPS data, stored in the FITS format, were first processed with the latest version of the standard NIRPS pipeline¹⁴. Primary data products were then extracted, which included the wavelength solution in vacuum, raw and telluric-corrected spectra, and properties of the cross-correlation function (CCF) such as its radial velocity (RV), full-width at half-maximum (FWHM), and bisector span. We note that the wavelength solution for the raw spectra and model telluric spectra, used to perform the telluric correction, were the same and were both expressed in the stellar rest frame. Subsequent steps, outlined in the sections below, focused on reducing the data to remove effects that degrade the spectra and have the potential to influence our analyses. The following reduction and analyses focused solely on data from Échelle order 14 (physical order 135) that contains the He triplet.

3.1.1 Telluric correction

Ground-based observations are prone to contamination from absorption and emission lines of various atmospheric species in Earth’s atmosphere. These features that can bias stellar spectra and subsequent RV measurements must be removed. While at low resolution individual telluric lines cannot be resolved and removed, this becomes possible when moving to high resolution.

In the NIRPS YJH bands, molecules such as H₂O, O₂, CO₂, and CH₄ have the most significant impact. The absorption linked to these molecules can be modeled in first order from the air-mass, except for water vapor, which requires knowledge of the integrated water column density and its variations linked to the weather conditions. The telluric correction used is included in the NIRPS processing pipeline and follows the methodology presented by Allart et al. (2022). It involves building a model transmission spectrum of Earth’s atmosphere from a line-by-line radiative transfer code assuming a single atmospheric layer. The spectral lines of modeled species are retrieved from the HITRAN database¹⁵. For each telluric species, a Levenberg-Marquardt algorithm minimizes the χ² between the CCF derived from the model using a few strong telluric lines of this species and the CCF derived from the observation using the same line selection. The resulting best-fit is used to remove the telluric lines. Allart et al. (2022) demonstrate a correction better than 2% peak to valley across spectral types F to K dwarfs with this pipeline. Although their tests are performed on visible spectra, the near-IR region around the He triplet is not significantly contaminated by tellurics and thus we expected a similar correction level. An example best-fit telluric spectrum is shown in Fig. 1.

3.1.2 Blaze, dispersion, and black body correction

We corrected each raw solar spectrum for the Blaze function using the corresponding Tungsten-lamp calibration spectrum. Following this correction, the flux of each spectrum presented a slow decreasing monotonous trend. Although the trend could have been ignored as it was negligible over the He triplet’s wavelength range, we removed it to avoid any bias that could influence later diagnostics or analyses. The trend arose from two sources, thus we removed both their contributions to ensure the robustness of our correction. We first accounted for the wavelength-dependent dispersion of the spectrograph (i.e., the fact that the width in wavelength of a pixel, and therefore the amount of light it collects, varies as a function of wavelength) as it had the strongest contribution to the decreasing monotonous trend. We used the DLLDX product, computed by the standard NIRPS pipeline and recorded in the FITS files, to detrend our spectra from this effect. This product records the size of each pixel in the wavelength domain. Secondly, because we are probing the Sun’s spectral energy distribution (SED) in the IR, we had to account for the shape of the Planck function over the wavelength range considered. We modeled the Sun’s SED at an effective temperature T_eff = 5778 K with Astropy’s BlackBody function and detrended our spectra with this model.

3.1.3 S/N, RV, and airmass clipping

We excluded spectra affected by poor weather, telluric-, instrumental-, or stellar-induced anomalies by applying sigma-clipping to the signal-to-noise ratio (S/N) and RV values. Spectra deviating from the daily average by more than 20σ were discarded. This threshold was selected because it was the lowest value that effectively removed outliers from all datasets. Although S/N and RV outliers often overlap, we applied both steps for thoroughness. Additionally, we omitted data from observations taken at airmass greater than two that can present additional telluric contamination.

3.1.4 Normalization and interpolation

We performed normalization to remove variations in the continuum level across different spectra, thus simplifying the subsequent steps of our analysis. To do so, we applied the RASSINE framework, outlined by Cretignier et al. (2020). Notably, we forced RASSINE’s alpha-shape algorithm to avoid the narrow region around the He triplet, thus guaranteeing that none of the variations intrinsic to the He triplet were absorbed by the normalization.

To accurately study the properties of the He triplet across different timescales, all spectra had to be interpolated on a common wavelength grid. RASSINE automatically performs this interpolation. Therefore, by normalizing the spectral time series for each day with RASSINE, we simultaneously interpolated all spectra onto a common wavelength solution.

3.2 He triplet variability analysis

Our analysis of variability aims at characterizing temporal changes in the helium triplet and elucidating the underlying causes. In the following sections, we highlight the different strategies employed to pursue this goal.

In particular, we conducted a preliminary assessment of the He triplet variability by comparing individual spectra to a master spectrum and visually inspecting the residuals in search of variations across the day of observation. Subsequently, we performed a rigorous analysis of the variability by fitting the He triplet line with a custom model and studying the variations in the fitted properties. Finally, we describe the different proxies that were used to investigate the origin of the He triplet variability.

3.2.1 Master spectrum

The master spectrum is a direct output of the RASSINE pipeline, built by summing all the raw spectra from the provided spectral time series and subsequently normalizing this sum. This normalization uses the same anchor point locations as those used to normalize the spectral time series (see Sect. 3.1.4). By comparing each normalized spectrum to the master we built a residual map, as shown in Fig. 2. This map presents structures around the He triplet and SiI lines that we analyzed to determine their origin. Shifting the maps from the stellar reference frame to the Earth’s barycentric frame showed no transformation of the diagonal tracks, seen around 10 830 and 10 832 Å, into vertical features, ruling out telluric or instrumental contamination. This indicated that the features had to be of solar origin. Additionally, residual maps from other observation days exhibited radically different structures, further supporting the conclusion that these features were driven by solar processes.

3.2.2 Radiative transfer

We obtained a more detailed characterization of the He triplet properties and their variations by modeling the He absorption in the solar chromosphere with a plane-parallel single-layer approximation. The equations used in this routine, which we describe below, are based on previous works aimed at studying the impact of radiation pressure on interstellar gas (e.g., Lagrange et al. 1998).

To model the observed line profiles of the solar spectrum (F_obs), we used the solution of the radiative transfer equation considering only extinction and ignoring the three-dimensional structure of the solar atmosphere. $$F_{\mathrm{o}\mathrm{b}\mathrm{s}}=F_0\times e^{-{\tau }_0}.$$(1)

Here, τ₀ is the solar atmosphere’s total optical depth, and F₀ is the base profile, taken as the spectrum at the base of the solar atmosphere. Given the scope of this project, we considered this approach to be a reasonable compromise between accuracy and computational efficiency. We considered a limited number of atomic transitions in the absorption process, namely, the He triplet transitions and the neighboring SiI transition. Including the SiI line ensured an accurate retrieval of the He triplet properties given that its red tail blends with the He triplet feature around 10 832.1 Å. As a result, the absorption, captured in the e^−τ₀ component and that is applied to the base profile F₀ provided as input to our model, is calculated as follows:

Since all transitions are considered simultaneously τ₀ = n_He × (σ_1,He + σ_2,He + σ_3,He) + n_Si × σ_Si, with n_He, n_Si the column densities of He and Si, and σ_1–3,He, σ_Si the cross-sections of the transitions considered. Each cross-section, exemplified by σ₀, is determined as follows: $${\sigma }_0(\lambda )=\frac{e^2}{4{ϵ}_0\times m_e\times c\times \sqrt{\pi }}\times \frac{f_{\mathrm{o}\mathrm{s}\mathrm{c}}\times {\lambda }_0\times L(\lambda )}{w}.$$(2)

The left-hand side contains fundamental constants with e the elementary charge, m_e the electron mass, ϵ₀ the vacuum permittivity, and c the speed of light. The right-hand side consists of terms depending on intrinsic properties of the electronic transition considered with f_osc the oscillator strength, λ₀ the reference wavelength in vacuum, L the line profile, and w the line width. While the $1/w\sqrt{\pi }$ factor is due to the normalization of the line profile, the λ₀ factor arises from approximating the cross-section as a Dirac delta function centered on the reference wavelength. Two line profiles are implemented in our model. A Gaussian profile, $$L(\lambda )=\exp (-\frac{1}{w^2}{(\frac{c}{{\lambda }_0}(\lambda -{\lambda }_0)+rv)}^2),$$(3)

and a Voigt profile, $$z=\frac{1}{w}(\frac{c}{{\lambda }_0}(\lambda -{\lambda }_0)+rv)+i\gamma$$(4)

L(λ) = Re[W(z)].

In Eq. (4), γ is the Voigt damping parameter and W is the Faddeeva function. In Eqs. (3) and (4), λ is the wavelength grid assigned to the input profile F(τ₀) and rv is a radial velocity offset. The line width, w, is defined as follows: $$w^2=\underset{\begin{array}{c}\text{\hspace{0.17em}Thermal\hspace{0.17em}}\\ \text{\hspace{0.17em}broadening\hspace{0.17em}}\end{array}}{\underbrace{\frac{2\times k_B\times T}{m}}}+\underset{\begin{array}{c}\text{\hspace{0.17em}Instrumental\hspace{0.17em}}\\ \text{\hspace{0.17em}broadening\hspace{0.17em}}\end{array}}{\underbrace{{\left(\frac{c}{R}\right)}^2}}+\underset{\begin{array}{c}\text{\hspace{0.17em}Rotational\hspace{0.17em}}\\ \text{\hspace{0.17em}broadening\hspace{0.17em}}\end{array}}{\underbrace{(v\cdot \sin (i){)}^2}}.$$(5)

Here, k_B is the Stefan-Boltzmann constant, T and m are the temperature and mass of the species considered, i and v are the stellar inclination and rotational velocity, and R is the instrumental resolving power. All variables are expressed in the International System of Units (SI). The equations for the helium and silicon transitions are identical to Eqs. (3), (4), and (5), except for the values of f_osc, λ₀, and m that were replaced by those specific to each species and transition. All the physical constants were fixed to values according to the SI. vsin(i) was set to the Sun’s projected rotational velocity (1.97 km/s, Lang 2013), and R was set to the average resolving power of the instrumental mode used (R_HE = 74 000). The values for f_osc and λ₀ were extracted from the NIST Atomic Spectral Database¹⁶.

By selecting an appropriate base profile for the solar spectra, such as the continuum, we fitted for these variables and retrieved the He triplet properties. Since we were fitting normalized solar spectra, the continuum profile was described with a constant flux level whose value was adjusted to each fitted spectrum.

3.2.3 Model construction and fitting

Using the model detailed in Sect. 3.2.2, we fitted the solar spectra to retrieve the He triplet properties. To ensure an accurate retrieval, we included the SiI line at 10 830 Å in the model. Indeed, this line’s red tail blends with the He triplet and could otherwise have biased the fit.

The Sil line exhibited a complex line shape that neither a Gaussian nor Voigt profile can capture individually. Thus, we adopted a hybrid approach, fitting its wings with a Voigt profile and its core with a Gaussian profile. In contrast, each He triplet line was well-described with a Gaussian profile. This setup led us to choose a fitting range, fixed for all spectra, that included the He triplet and SiI lines as much as possible while excluding other neighboring lines. The selected fit range is illustrated in Fig. 1.

The final model included eight free parameters: shared temperature, column density, and radial velocity offset for all three He triplet transitions, separate values for the Sil transition, a damping parameter for the SiI wings’ Voigt profile, and a multiplicative factor controlling the continuum level. We fitted each telluric-corrected, normalized solar spectrum using Scipy’s curve_fit function. As bounds were set for all free parameters, the Trust Region Reflective algorithm was needed to perform the minimization step of the nonlinear least squares method.

Given our hybrid approach to modeling the Sil line, we needed to determine the wavelength borders that defined the core and wings of the line. To do so, we first computed the Gaussian and Voigt profiles for the SiI line with the set of parameters at the minimization step considered. Using a least-square approach, we then found the wavelengths on either side of the SiI line core that minimized the overlap between the Gaussian and Voigt profiles. These wavelengths were used to compute a final model describing the full profile of the Sil line. We note that minimizing the overlap between the core and wings was crucial to creating a model that was as continuous as possible across the full SiI line.

After fitting, He triplet and SiI EWs were calculated for each normalized solar spectrum by integrating the flux within the wavelength ranges of the respective lines. The wavelength range was computed from the absolute transition wavelength (shown in Fig. 1) corrected from the fitted radial velocity offset of the He triplet and Sil lines. The continuum level used to perform this integration was set to the best-fit value of the multiplicative factor controlling the continuum profile.

Additionally, we improved our estimate of the best-fit parameter errors using a bootstrap approach. To do so, we generated k=500 realizations of each solar spectrum by sampling a Gaussian distribution with mean and standard deviation corresponding to the value and error on the flux of each point in the spectrum. The uncertainty on each original best-fit parameter was approximated by its corresponding standard deviation across the best-fits of the k realizations. Finally, for a more tangible measure of the He triplet’s width and depth, we converted the best-fit He triplet temperatures and column densities to line widths in km/s and contrasts in %. The line widths were obtained by evaluating Eq. (5) with the best-fit temperatures and retrieving w. To determine the contrast for each solar spectrum, we subtracted the flux at the center of the model He triplet from the continuum flux level, retrieved from the fitted multiplicative factor. We studied the evolution of the best-fit depth, position, width, and EW for the He triplet by plotting their time series as illustrated in Figs. 3 and 4. These measurements are discussed in detail in Sect. 4.1.

3.2.4 Variability diagnostics

Changes in the He triplet properties measured in our spectra could potentially be of telluric and/or stellar origin. In an effort to identify the true underlying cause, we examined proxies associated with each of these mechanisms.

Telluric contamination due to imperfect correction could arise from temporal variations in atmospheric components. In particular, the presence of a telluric water line overlapping with the He triplet prompted us to explore this possibility (see Fig. 1). These variations can be assessed by studying the total amount of a particular species present in the vertical column of air above the telescope, also known as the integrated column density (ICD). We retrieved the ICD for H₂O, O₂, and CO₂ from the NIRPS pipeline FITS files. As expected, the time series of these quantities were strongly correlated with the airmass. Thus, to study the unbiased effects of ICDs we detrended them using the airmass measurements provided in the NIRPS pipeline FITS files. These same measurements were used to explore the impact of airmass on the He triplet. In addition, we also fitted the water line overlapping with the He triplet with the following Gaussian profile to extract and subsequently study its properties.

$$G=(1-a\times \exp (-\frac{1}{2}{\left(\frac{x-\mu }{\sigma }\right)}^2))\times (b\times x+c),$$(6)

where a, μ, and σ are the depth, position, and width of the line, while b and c are the first-order polynomial coefficients describing the local continuum. The fitting procedure was performed with the curve_fit function. Due to the weakness of this water line with respect to the He triplet, we fitted its profile from the telluric model instead of the raw spectra.

In a second phase, we considered the effects of stellar activity on the He triplet. We first studied the SiI line at 10 830 Å A, whose properties were extracted with our custom model. Studying this line is crucial as its known stability and lack of nearby telluric lines suggest that variability in its best-fit parameters is associated with stellar activity.

The CCF, commonly used to study the overall spectrum or evaluate the quality of the normalization procedure, is also a proxy for stellar activity. The CCF is a mathematical technique that cross-correlates the observed spectrum to a template composed of the most prominent spectral lines, giving as output an average profile of all considered spectral lines (Baranne et al. 1996; Pepe et al. 2002). CCF properties, such as the FWHM and bisector span, can be used to evaluate the impact of magnetic stellar activity induced by the presence of active regions on the solar surface (e.g., Queloz et al. 2001; Desort et al. 2007; Dumusque et al. 2014) or granulation and supergranulation (e.g., Dumusque 2016; Cegla et al. 2019). However, it is important to acknowledge that factors other than stellar activity can impact the CCF. Instrumental effects, such as variations in the point-spread function or detector nonlinearity, alongside variations in atmospheric conditions, can induce changes in CCF properties like the FWHM and contrast.

We evaluate the impact of each mechanism by performing a correlation study. This is done by plotting the time series of each helium triplet best-fit parameter against all other quantities mentioned above. We note that our correlation studies include the HeI line’s depth, width, and EW, rather than relying solely on the EW. The EW is the spectral line’s area, and, therefore, variations in this quantity reflect variations in the line depth and width. These can vary in compensating ways, such as the line getting both deeper and narrower. Thus, relying only on the EW can lead to misleading correlation analyses, as a lack of correlation with the EW might mask significant variations in the line depth or width. For example, in Table E.1, correlations with the HeI depth, but not with the EW, indicate compensating variation in depth and width.

We computed Pearson and Spearman correlation coefficients for each comparison, which can capture linear and nonlinear monotonic correlations, respectively. We did so using the spearmanr and corrcoeff functions from Scipy and Numpy, respectively. Example correlation plots are shown in Appendix E. The implications of these correlations are discussed in Sect. 4.1.

3.2.5 Long timescales

After measuring short-term variations, our focus shifted to investigating the He triplet’s behavior over longer timescales. We consolidated individual time series for each He triplet property over the 12 days of observations into a unified time series. We performed correlation studies between these unified time series and the $R_{HK}^{\mathrm{\prime }}$ activity index. This index is associated with variations in the cores of the Ca II H and K lines that are closely linked to changes in the stellar magnetic field (e.g., Wilson 1968; Borgniet et al. 2015; Milbourne et al. 2021; Cretignier et al. 2024). As these lines are in the visible wavelength range, the index was derived from HARPS solar spectra. However, since HARPS and NIRPS observations were conducted at different times (from 1.8 to 7.3 hours apart), the $R_{HK}^{\mathrm{\prime }}$ values were first extrapolated onto the corresponding NIRPS observation windows (see Table A.1) before performing our analysis. To achieve this, we performed a Fourier expansion of the activity index, keeping up to the third term. The fundamental frequency used in the expansion was obtained from a periodogram analysis and corresponded to the synodic solar rotation period. The results of our correlation studies are shown in Fig. 5.

3.3 Injection-recovery tests

In the upcoming sections, we define the methods used to generate a planetary signal, introduce it into observed solar spectra, and subsequently retrieve it. We aimed to evaluate how variations in the solar He triplet over various timescales impacted a simulated planetary signal.

3.3.1 Signal generation

By making similar simplifying assumptions (i.e., ignoring the exoplanet orbital geometry, and the atmosphere’s emission and three-dimensional structure), we adapted the model detailed in Sect. 3.2.2 to simulate planetary signals. As mentioned previously, this model applies an absorption signature to a base profile. In contrast to our previous implementation, where the stellar continuum served as the base profile, we used the processed solar spectra as the base profile. In doing so, the introduced absorption signature acted as a mock planetary signal. While the planetary signal should ideally be introduced in the “true” solar spectra before performing instrumental convolution, we believe introducing it afterward yields comparable results especially since our model is a first-order approximation. Broadening induced by the exoplanet’s rotation was ignored since it was negligible compared to the other sources of broadening.

Unlike the model defined in Sect. 3.2.2, the planetary absorption signature applied to the processed solar spectra only included contributions from the planetary He triplet. Gaussian profiles were used to model the three He triplet transitions. These profiles shared the parameters controlling their temperature (T) and column density (n). The T and n values of the simulated planetary signal were fixed to typical exoplanet thermosphere values, as discussed in Sect. 4.3. The radial velocity shift, rv, was set to zero for simplicity.

3.3.2 Signal introduction and retrieval

Having established a method to generate a mock planetary signal, we now describe how the signal is incorporated into the solar spectra, and its properties are subsequently retrieved. We split the array of processed solar spectra into two sets, in-transit and out-of-transit, by defining a transit duration during which the simulated planet passes in front of the Sun. For observations taken over the time interval [t, t + Δt], we defined $[t+\frac{1}{4}\mathrm{\Delta }t,t+\frac{3}{4}\mathrm{\Delta }t]$ as in-transit and $[t,t+\frac{1}{4}\mathrm{\Delta }t]\cup [t+\frac{3}{4}\mathrm{\Delta }t,t+\mathrm{\Delta }t]$ as out-oftransit. In our case, t denotes the timestamp at the start of a given day, and Δt is the total duration of observations for that day. To further simplify our approach, we ignored the drop in flux during transit as it is locally achromatic, and assumed the planetary signal to be static in time. Additionally, we ignored planet-occulted line distortions as their impact on derived planetary atmosphere parameters is small, thanks to the Sun’s low vsin(i) (Dethier & Bourrier 2023).

After introducing the planetary signature into the in-transit set, we built master in-transit and out-of-transit spectra. Each master spectrum was created by taking the weighted average of the processed solar spectra in the associated set, with weights corresponding to the inverse squared flux errors. A planetary absorption spectrum, A, was extracted for each in-transit spectrum, F_in, with the following formula: $$A=\frac{{\overline{F}}_{out}-F_{in}}{{\overline{F}}_{out}},$$(7)

with ${\overline{F}}_{out}$ the master out-of-transit spectrum. To retrieve the simulated planet’s He triplet properties, which may have been affected by the solar He triplet, we fitted each planetary absorption spectrum with the model described in Sect. 3.3.1, keeping T and n free. Since n and T can reach large values, we fitted their logarithms instead. An example simulated absorption spectrum is shown in Fig. 6. We plot the time series of the resulting best-fit parameters in Figs. 7 and 8. To validate the results, we also retrieved the best-fit parameters for the master planetary absorption spectrum, obtained by applying Eq. (7) to the master in-transit spectrum. Finding that both master best-fit parameters were within <0.01σ of the expected results, we did not include them in our figures. By comparing the weighted average of the retrieved parameters to the true values used as input for the simulation, we evaluated the impact of the solar He triplet variability on retrieved planetary parameters.

We performed these injection-recovery tests both on individual days and across the entire time series of observations. This approach allowed us to study the influence of short- and long-term He triplet variations on the retrieved exoplanet properties. For the full time series’ injection-recovery test, the master spectra were constructed on a per-day basis. Specifically, we treated each day as an independent transit event, using its corresponding out-of-transit and in-transit spectra to generate the associated absorption spectra. We present and discuss the results of our various injection-recovery tests in Sect. 4.3.

4 Results and discussion

We showcase our results with observations from January 26, 2023 and January 27, 2023. We select these days based on the abundance and quality of measurements, as well as the significant differences that we see in their correlation analyses that highlight the variable nature of solar He triplet measurements. We present images of the Sun’s surface for those days, captured at 1700 Å with the Atmospheric Imaging Assembly (AIA) of the SDO (Lemen et al. 2012), in Figs. D.1 and D.2.

4.1 Short-term variability

As shown with the residual map presented in Fig. 2, we observe notable variability in the solar He triplet. Further support for this variability is obtained by fitting the solar He triplet and examining the time series of the fit parameters, as shown in Fig. 3. In particular, we observe variations on two distinct timescales: rapid oscillations spanning several minutes and a longer trend over the entire day of observation. We witness similar variations for all days of observations.

In our analysis of the variability (see Sect. 3.2.4), we explore two potential sources: stellar activity and telluric contamination. We obtain polarizing results from correlation studies between the solar He triplet and the selected stellar activity probes. While certain days exhibit strong correlations (r_S > 0.5) with multiple CCF properties, such as the FWHM, bisector span, and RV, we consistently find no significant correlations between the properties of the SiI line and the He triplet. The latter is expected, as these lines form in different regions of the solar atmosphere and the SiI line is generally stable. Interestingly, preliminary analyses revealed strong He-Si correlations, particularly in their EWs, that were highly sensitive to the normalization method used. In Appendix B we explore this dependence and its possible origins, and discuss our motivation for adopting RASSINE normalization in our fiducial analysis. The differences in correlations with stellar activity probes across different days, or lack thereof, are exemplified in Figs. E.1, E.2, E.3, and E.4 (CCF), and Figs. E.5, E.6, E.7, and E. 8 (Si).

Additionally, some days show strong correlations between He triplet properties and multiple telluric contamination tracers (e.g., airmass, telluric water line properties, O₂, and CO₂ ICD). Remarkably, these are the same days that exhibit strong correlations with the CCF properties (see Table E.1). This suggests that improper telluric correction may contribute to the observed He triplet variability. Indeed, telluric contamination, traced with the airmass and ICDs, can impact CCF properties such as position and shape. Furthermore, low-amplitude telluric absorption lines, such as the water line near the He triplet, are challenging to account for with the NIRPS telluric correction pipeline. This limitation can result in the He triplet varying similarly to the telluric water line properties. These results are exemplified in Fig. E.9 (airmass), Figs. E.10, E.11, and E.12 (ICDs), and Figs. E.13, E.14, E.15, and E.16 (water line). We further emphasize the link between CCF properties and telluric contamination in Table E.2. Notably, the CCF properties most strongly correlated with the He triplet (i.e., CCF RV and FWHM) also show the strongest correlations with the telluric contamination proxy. Moreover, the lack of significant correlations on January 26, 2023 contrasted with their presence on January 27, 2023 mirrors the pattern seen in Table E.1. This provides additional evidence that the observed correlations (or lack thereof) between He triplet properties and CCF properties arise from telluric contamination.

Despite identifying telluric contamination as a possible source for the solar He triplet variability, our analysis examined each mechanism independently. Investigating the interactions between these processes through comprehensive modeling would provide a more reliable explanation for the observed variations.

Although our analysis focuses on the solar He triplet variations across an entire day of observation, short-term oscillations, lasting from several minutes to an hour, are also present (see Figure 3). These rapid oscillations were not studied in detail due to their negligible impact on atmospheric retrievals. Indeed, they are typically smoothed out when combining multiple transit observations. They may be of stellar or instrumental origin, warranting further investigation.

4.2 Long-term variability

As illustrated in Fig. 5, the solar He triplet properties also exhibit variations over the entire duration of our study. Both the depth and EW of the He triplet present significant correlations with the $R_{HK}^{\mathrm{\prime }}$ index. Due to the pertinence of these properties in atmospheric studies (e.g., depth provides measurements of the exospheric helium content), we provide linear relationships linking the $R_{HK}^{\mathrm{\prime }}$ index to stellar He triplet depth and EW so that future studies can evaluate or predict the impact of stellar activity on observed stellar He triplet variations. $$\begin{array}{rl}\mathrm{E}\mathrm{W}(\text{Å})& =(3.8\pm 0.8\text{Å})\cdot R_{HK}^{\mathrm{\prime }}+(34\pm 4\text{Å})\\ {\log }_{10}n& =(0.5\pm 0.1)\cdot R_{HK}^{\mathrm{\prime }}+(2.5\pm 0.6).\end{array}$$(8)

Since $R_{HK}^{\mathrm{\prime }}$ is sensitive to variations in magnetic field, this correlation could be attributed to rotational modulation due to magnetic regions on the Sun’s surface. These active regions affect the flux at the photosphere and the shape of spectral lines (e.g., Basri et al. 1989; Shcherbakov et al. 1996; Cretignier et al. 2021, 2023), thus leading to changes in the line depth and EW. These findings align with Guilluy et al. (2020), confirming the He triplet’s sensitivity to variations in magnetic field. As illustrated in Figs. D.1 and D.2, the witnessed rotational modulation can be attributed to the presence of faculae on the solar surface.

Our long-timescale results diverge from those presented by Krolikowski et al. (2024). Although our EW’s mean absolute deviation (MAD_EW = 2.35 mÅ) agrees with their power-law fit in Fig. 8, we find opposing results when investigating the origin of the EW variations. Indeed, while we find a strong correlation between the He triplet EW and the solar rotation (as traced by the $R_{HK}^{\mathrm{\prime }}$ index), a majority of stars in their sample do not present a similar correlation. Their sole outlier, V1298 Tau, shows a strong correlation with stellar rotation over a brief time range that they attribute to an intense flaring or mass loss event. Such discrepancies in our results may arise from differences in the stellar populations considered. The Sun, an old and relatively quiet star, could have its long-timescale He triplet variability attributed to star spots and faculae, whereas He triplet variability in young stars may be dominated by other stellar activity processes unaffected by rotational modulation. In light of these results, we recommend that future studies complement near-IR spectroscopic observations with concurrent photometric monitoring and optical spectroscopic observations of known activity tracers. This approach would help identify and characterize potential stellar He triplet variability. As emphasized by Krolikowski et al. (2024), this is particularly important for observations of young stars.

4.3 Impact on planetary atmosphere retrieval

We investigate the impact of stellar He triplet variability on retrieved planetary parameters through injection-recovery tests. The results, presented in Table C.1, demonstrate that the expected parameter values consistently fall within the 1σ region for all the days of observation, where σ is the weighted standard deviation of the retrieved parameters. An example showcasing these results is shown in Fig. 7. The baseline case used to study the response over all the days of observation is T = 10 000 K and n_He = 5 × 10¹⁴ cm⁻², inducing an absorption signature of ~5%. These values reflect conditions in the upper atmospheric layers of hot giant planets, which are the main targets of atmospheric escape studies due to their significant transmission spectroscopy signatures. These layers can reach temperatures ranging from 5000 to 20 000 K (e.g., Guilluy et al. 2020; Lampón et al. 2023) and, being primarily composed of light elements such as H and He, can induce absorption signatures of several percent (e.g., Allart et al. 2018, 2019), making the choice of input conditions realistic.

From these observations, and the values presented in Table C.1, we conclude that the solar He triplet variability, present on both ultra-short timescales (~10 min) and over the entire day of observation (as highlighted by the binning in Fig. 7), does not have a significant impact on the retrieved planetary parameters. Thus, assuming the stellar He triplet is stable throughout the transit duration remains a valid approximation. Furthermore, our model, although a first-order approximation reliant on basic assumptions, successfully demonstrates the minimal impact of solar He triplet variability on the retrieved planetary parameters. This implies that introducing additional complexity to the model would lead to the same conclusion.

Our injection-recovery analysis of long-term He triplet variability, presented in Fig. 8 and the last row of Table C.1, yields similar results: the expected parameter values remain within the corresponding 1σ region. This result suggests that variations in the He triplet over days to months, due to stellar activity, do not significantly affect the retrieval of exoplanet properties.

We stress that our results are only valid for our analysis of a 5% planetary signature in the presence of solar variability. Further analyses would have to be conducted to evaluate the impact of stellar He triplet variability when more or less active stars and stronger or weaker planetary signatures are considered.

5 Conclusion

Our study revealed the presence of variability in the solar He triplet, a spectral line recently used in the study of atmospheric escape. We demonstrate that the He triplet exhibits variations across different timescales under the influence of various processes. We identify improper telluric correction as the potential source of short-term He triplet variability and stellar activity as the definitive source of long-term He triplet variability. We elected not to investigate the impact of instrumental noise (e.g., modal noise) on our results as it is hard to quantify, and its impact is likely negligible. Despite having identified sources for the He triplet variability, we recommend that additional studies be performed accounting for correlated noise and the combined effects of stellar activity and telluric contamination, for example using the ANTARESS workflow (Bourrier et al. 2024) to reduce and analyze the helium spectra. Further studies could also refine the long-timescale relations that we established between the $R_{HK}^{\mathrm{\prime }}$ index and the He triplet depth and EW (see Eq. (8)).

Importantly, our findings indicate that variations in the solar He triplet, over the timescale of a single transit or across multiple transits, have minimal impact on the retrieved planetary parameters crucial to the study of atmospheric escape. Moreover, given our use of a first-order model, we expect that more complex modeling approaches would reach similar conclusions. While we definitively observe strong variations induced by stellar activity on the timescale of solar rotation, as seen in Guilluy et al. (2020), our analysis demonstrates that these variations are well captured in the master-out spectrum and would therefore not impact studies focusing on the variability of exoplanet atmospheres through time.

These results are significant as the He triplet is gaining in popularity and availability for atmospheric studies, especially with the arrival of new instruments and telescopes, like NIRPS and NIGHT, that specifically target the IR range to measure planetary He signals and characterize stellar He variability. However, we note that the conclusions drawn here are only valid for old, quiet Sun-like stars. It is well-known that stellar activity varies significantly across different spectral types and stellar ages. As a result, the He triplet variability would need to be reevaluated when considering M-dwarfs, for example.

Data availability

The plots and tables illustrating the correlation studies between the He triplet properties, stellar activity tracers, and telluric contamination tracers, referred to as Appendix E, can be found online at https://zenodo.org/records/14674319.

Acknowledgements

S.J.M. acknowledges support from the Massachusetts Institute of Technology through the Praecis Presidential Fellowship. X.D. 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). 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. 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. M.Cr. acknowledges the Swiss National Science Foundation funding under grant P500PT 211024. W.D. and N.C.S. 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. 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. R.A. 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). R.A., F.Ba. and R.D. 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. F.Ba. and R.D. acknowledge 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). A.C. acknowledges funding from the French ANR under contract number ANR18CE310019 (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. Our solar data analysis made use of the following packages: Astropy¹⁷ (v5.3.4, Astropy Collaboration 2013, 2018, 2022), Matplotlib¹⁸ (v3.7.2, Hunter 2007), Numpy¹⁹ (v1.26.1, Harris et al. 2020), Pandas²⁰ (v1.5.3, Wes McKinney 2010), Scipy²¹ (v1.11.3, Virtanen et al. 2020).

Appendix A HELIOS observations

Appendix B He-Si EW correlation

The properties of the solar He triplet, formed in the upper chromosphere (Goldberg 1939; Zirin 1975), and the SiI line, formed in the upper photosphere (Bard & Carlsson 2008), are not expected to exhibit strong correlations. Nevertheless, our preliminary analyses reveal a significant correlation in their equivalent widths. Moreover, we find that this correlation depends on the normalization method used. This puzzling result prompted us to explore various normalization methods to determine the optimal approach for normalizing spectra and understand the cause of the observed correlation’s dependence on the chosen method.

In our initial approach, the stellar continuum around the He triplet of each un-normalized solar spectrum is described with an n-th order polynomial. This polynomial acts as the base profile onto which the absorption signature is applied in our model (outlined in Sect. 3.2.2). As reported in Table B.1, we retrieve significant, linear or nonlinear, He-Si EW correlations when using a zeroth-order polynomial. Despite obtaining low correlations with the first- and second-order polynomials, both correlation plots display noticeable structure (see Fig. B.1). We conclude that due to the complexity, particularly the symmetry, of this structure, the Pearson and Spearman coefficients fail to accurately capture the strength of the correlation. The variation in correlation coefficients across different polynomial orders suggests that the observed correlation is not physical but an artifact of our normalization procedure. This prompts us to adopt a more refined approach: each solar spectrum is first normalized with the RASSINE pipeline (Cretignier et al. 2020) before fitting it with our model.

As mentioned in Sect. 3.1.4, RASSINE utilizes an alpha-shape algorithm to scan spectra and identify local maxima, also known as anchor points. A key parameter in the alpha-shape algorithm is the rolling pin radius. This user-selected parameter can have important repercussions, as too small a rolling pin can lead to anchor points being placed inside stellar lines. This can result in the normalization procedure absorbing variations intrinsic to the stellar lines, thus biasing our results. Selecting an optimal rolling pin radius is therefore a crucial step in the normalization procedure. Thankfully, the rolling pin radius can be adjusted, and its impact on the anchor points can be observed, through the RASSINE graphical user interface (GUI). Using RASSINE’s GUI to determine the optimal pin radius thus ensures that anchor points are not placed inside the He triplet and SiI lines. Additionally, we perform an extra cautionary step to verify that RASSINE’s normalization does not absorb any of the variations in the He triplet and SiI lines. RASSINE allows users to manually exclude wavelength ranges, preventing anchor points found in this range from contributing to the continuum profile. Combining this option with the optimal rolling pin radius we find no significant He-Si EW correlation (see Table B.1). In doing so, we prove that RASSINE’s normalization method, using the optimal rolling pin radius, does not significantly absorb variations intrinsic to the He triplet and Si lines. The results presented in Table B.1, and Fig. B.1 are for observation day January 26, 2023.

Appendix C Planetary signal simulation

Appendix D SDO images

References

---

All Tables

All Figures