# Measuring Sub-Kelvin Variations in Stellar Temperature with High-resolution Spectroscopy

Étienne Artigau, Charles Cadieux, Neil J. Cook, René Doyon, Laurie Dauplaise, Luc Arnold, Maya Cadieux, Jean-François Donati, Paul Cristofari,

Xavier Delfosse, et al.

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

To cite this version:

Étienne Artigau, Charles Cadieux, Neil J. Cook, René Doyon, Laurie Dauplaise, et al. Measuring Sub-Kelvin Variations in Stellar Temperature with High-resolution Spectroscopy. The Astronomical Journal, 2024, 168, 10.3847/1538-3881/ad7b30. insu-04835922

HAL Id: insu-04835922

https://insu.hal.science/insu-04835922v1

Submitted on 13 Dec 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Measuring Sub-Kelvin Variations in Stellar Temperature with High-resolution Spectroscopy

Étienne Artigau¹,² ID, Charles Cadieux¹ ID, Neil J. Cook¹ ID, René Doyon¹,² ID, Laurie Dauplaise¹ ID, Luc Arnold³ ID, Maya Cadieux¹, Jean-François Donati⁴ ID, Paul Cristofari⁵ ID, Xavier Delfosse⁶ ID, Pascal Fouque⁴ ID, Claire Moutou⁴ ID, Pierre Larue⁷ ID, and Romain Allart¹ ID

¹ Institut Trottier de recherche sur les exoplanètes, Département de Physique, Université de Montréal, Montréal, Québec, Canada; etienne.artigau@umontreal.ca

² Observatoire du Mont-Mégantic, Québec, Canada

³ Canada-France-Hawaiʻi Telescope, CNRS, Kamuela, HI 96743, USA

⁴ CNRS, OMP, Université de Toulouse, 14 Avenue Belin, F-31400 Toulouse, France

⁵ Center for Astrophysics—Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

⁶ Observatoire de Genève, Département d'Astronomie, Université de Genève, Chemin Pegasi 51b, 1290 Versoix, Switzerland

⁷ Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

Received 2024 June 17; revised 2024 September 10; accepted 2024 September 11; published 2024 November 12

## Abstract

The detection of stellar variability often relies on the measurement of selected activity indicators, such as coronal emission lines and nonthermal emissions. On the flip side, the effective stellar temperature is normally seen as one of the key fundamental parameters (with mass and radius) to understanding the basic physical nature of a star and its relation with its environment (e.g., planetary instellation). We present a novel approach for measuring disk-averaged temperature variations to sub-Kelvin accuracy inspired by algorithms developed for precision radial velocity (pRV). This framework uses the entire content of the spectrum, not just preidentified lines, and can be applied to existing data obtained with high-resolution spectrographs. We demonstrate the framework by recovering the known rotation periods and temperature modulation of Barnard star and AU Mic in data sets obtained in the infrared with SPIRou at CHFT and at optical wavelengths on ϵ Eridani with HARPS at ESO 3.6 m telescope. We use observations of the transiting hot Jupiter HD189733b, obtained with SPIRou, to show that this method can unveil the minute temperature variation signature expected during the transit event, an effect analogous to the Rossiter-McLaughlin effect but in temperature space. This method is a powerful new tool for characterizing stellar activity, and in particular temperature and magnetic features at the surfaces of cool stars, affecting both pRV and transit spectroscopic observations. We demonstrate this method in the context of high-resolution spectroscopy but it could be used at lower resolution.

*Unified Astronomy Thesaurus concepts: Stellar effective temperatures (1597); High resolution spectroscopy (2096); Time series analysis (1916)*

*Materials only available in the online version of record: machine-readable table*

## 1. Introduction

The ability to detect and analyze exoplanets has ushered in an era where characterizing their atmospheres, surfaces, and potential habitability stands at the forefront of astronomical pursuits. The study of exoplanets is performed through several technical means. Radial velocity measurements, involving the detection of the kinematic effect on the host star by the Keplerian motion of its planets, have been the first successful discovery technique. Although this method is still used to identify a large number of planets, it has been overtaken in the sheer number of discoveries by photometric surveys where one looks for the minute flux change during a transit when a planet passes in front of its parent star. While precision radial velocity (pRV) and photometric surveys have contributed to the discovery of the majority of known exoplanets to date, both methods remain strongly biased toward the discovery of close-in planets (<0.1 au for transit, <5 au for pRV). Direct imaging surveys with adaptive optics, as well as very deep imaging in the infrared, complement the current portrait of exoplanets with

companions now known at separations up to thousands of astronomical units.

As exoplanet pRV searches gained momentum, they were confronted with a series of challenges, one of the most significant being the intricate interplay between the activity of host stars and the Keplerian motion of exoplanets. This complex relationship introduces variability and noise at all timescales, ranging from minutes for oscillation (B. J. Brewer et al. 2007; M. Bazot et al. 2012) and granulation (T. Kallinger et al. 2014; H. M. Cegla et al. 2018) to years for magnetic cycles (S. L. Baliunas et al. 1995; J. Gomes Da Silva et al. 2012). Activity filtering techniques have emerged as crucial tools to mitigate the confounding effects of stellar activity on exoplanet observations, ensuring that signals emerging from stellar variability are not falsely interpreted as kinematic signatures of exoplanets (e.g., Alpha Cen B b, V. Rajpaul et al. 2016; Gl699 b, J. Lubin et al. 2021; AD Leo b, A. Carmona et al. 2023; LHS 1140d, C. Cadieux et al. 2024b; TW Hya b, N. Huélamo et al. 2008).

High-accuracy photometry, ideally obtained contemporaneously with radial velocity (RV) measurements, informs on the distribution of features on the photosphere that depart from the mean brightness of the star. A passing stellar spot or a plage on a rotating star will lead to an RV jitter as it affects the mean

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
line shape while also imparting a photometric signature. One can show that at first order, for a flux $f$ and its time derivative $f'$, the RV signature of the spot and plages will scale as $f \cdot f'$ (S. Aigrain et al. 2012). An additional term is added to account for the convective blueshift inhibition (N. Meunier et al. 2010a, 2010b), leading to a correspondence between the flux and its time derivative and RV jitter. Convective blueshift inhibition refers to the reduction of the net blueshift effect caused by convective motions in a star's photosphere, often due to factors such as magnetic activity. This inhibition can introduce noise and systematic errors into RV measurements and becomes dominant against rotation-induced jitter for slowly rotating stars.

Several spectroscopic activity indicators have been used, each tracing a different facet of the activity and, ultimately, statistically tied to the perturbing RV jitter.

The indicators inform on the underlying physical common cause between an ancillary observable and the line barycenter shifts that mimic Keplerian reflex motions. Tracers of line profile distortions were first measured from the cross-correlated function (CCF) representing the averages of many spectral lines. These indicators include the bisector inverse slope (D. Queloz et al. 2001) measuring the skewness/asymmetry of the CCF (see also $V_{span}$; I. Boisse et al. 2011) or the full-width at half maximum (FWHM; D. Queloz et al. 2009) tracing the CCF width change. Inferred as a differential measurement from a stellar template spectrum, the chromatic index and the differential line width indicators also reflect variation in line shape from stellar activity (M. Zechmeister et al. 2018; É. Artigau et al. 2022).

R. D. Haywood et al. (2022) estimated the Sun's disk-averaged RV variations using the Fe I line observed by Solar Dynamics Observatory/HMI and modeled the magnetic activity, showing that a linear fit to the unsigned magnetic flux reduces RV variation by 62%. Their findings suggest that unsigned magnetic flux is an excellent proxy for activity-induced RV variations. Other tracers are based on emission in the core of specific lines, such as the S-index and $\log R'_{HK}$ indicators for the Ca II H and K lines on Sun-like stars (R. W. Noyes et al. 1984), or $H\alpha$ line for cooler stars (M. Kürster et al. 2003). These chromospheric indicators may, however, only weakly correlate with the RV jitter that arises from photospheric activity (P. Schöfer et al. 2019).

We propose a novel technique to measure stellar activity, namely a precise differential temperature measurement based on a high-resolution spectrum, measuring slight variations in the disk-averaged temperature due to the passage of active regions, either hotter than the photosphere (i.e., plages) and/or cooler (i.e., stellar spots). This technique has the advantage of applying to all existing pRV archival data. While the technique was devised with exoplanet searches in mind, it can be applied to data sets obtained with high-resolution spectrographs that are not stabilized at the m s⁻¹ level.

This paper is structured as follows. In Section 2, we present the mathematical framework of the method. We then describe the data sets used to build the spectral templates necessary for this method in Section 3. In Section 4, we apply this framework to data sets obtained with an infrared spectrograph (SPIRou; J.-F. Donati et al. 2020) and an optical spectrograph (HARPS; F. Pepe et al. 2002; M. Mayor et al. 2003) in the context of pRV monitoring and transit spectroscopy. We explore the change in the disk-averaged temperature during a transit in

Section 5. The proposed method being based on observations, we explore the caveats of using stellar models in our analysis in Section 6. We summarize our findings and highlight exciting future research avenues in Section 7.

## 2. Method

Our approach to temperature measurement rests on the line-by-line (LBL) technique described in É. Artigau et al. (2022) and expands its usage to the temperature domain rather than velocity (K. Al Moulla et al. 2022). The LBL technique involves measuring velocity, or equivalently small wavelength shifts, from the first derivative expansion of a nearly noise-free template of a star matched to a given observed spectrum. Crucially, this technique operates independently over individual spectral lines, allowing for the identification and removal of spectral outliers that would otherwise bias the velocity measurement. This technique was initially designed with stellar activity filtering in mind (X. Dumusque 2018; M. Cretignier et al. 2020) but also proved to be efficient at removing instrumental signatures (M. Cretignier et al. 2021; M. Ould-Elhkim et al. 2023). While a powerful tool at optical wavelengths, the LBL technique presents an even greater interest in the near-infrared where deep telluric absorption affects RV measurements (P. Figueira et al. 2016; S. X. Wang et al. 2022) and detector cosmetics are notably worse than for CCDs (E. Artigau et al. 2018a). Near-infrared LBL velocity measurements reach the <2 m s⁻¹ accuracy (É. Artigau et al. 2022, 2024) and are, to date, the only method consistently providing constraints on planetary signatures at near-infrared wavelengths at the m s⁻¹ level and better (e.g., C. Cadieux et al. 2022; T. Gan et al. 2022; C. Moutou et al. 2023).

We introduce an adaptation of the LBL framework by substituting the velocity derivative with the temperature derivative of the spectrum. This modification enables a differential evaluation of temperature relative to a template. Leveraging the full spectrum's information content, this adapted approach prioritizes precision while minimizing error.

Central to the LBL algorithm is a formalism introduced by F. Bouchy et al. (2001), in which one finds the shift $\delta\lambda(i)$ in wavelength space (and correspondingly, in velocity space) between a spectrum $A(i)$ and its (ideally) noiseless template $A_0(i)$ for pixel $i$ as:

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

where the $\frac{\partial A_0(i)}{\partial \lambda(i)}$ term corresponds to the velocity gradient of the template (i.e., change in flux for an infinitesimal velocity offset). This can be seen as the first term of a Taylor expansion of flux against wavelength. This can be expressed as a velocity gradient by using

$$\frac{\delta V}{c} = \frac{\delta \lambda}{\lambda}, \quad (2)$$

where $\delta V$ is the velocity shift, $\delta \lambda$ is the wavelength shift, and $c$ is the speed of light. Within this framework, one can take the derivative of the spectrum against any relevant quantity. Here, we explore the derivative in temperature, $\frac{\partial A_0(i)}{\partial T(i)}$. While the derivative against velocity can be determined from a high signal-to-noise ratio (SNR) template, the derivative in temperature requires a library of spectra of objects with known
temperatures. By analogy to the velocity-domain weighting of residuals presented in F. Bouchy et al. (2001), one can derive a per-pixel temperature change expressed as:

$$ \delta T(i) = \frac{A(i) - A_0(i)}{\partial A_0(i)/\partial T(i)}, \qquad (3) $$

and a corresponding error in the temperature as derived for that spectral element;

$$ \sigma_T(i) = \frac{\sigma(i)}{\partial A_0(i)/\partial T(i)}. \qquad (4) $$

Here, $\sigma(i)$ is defined as the running standard deviation of the difference between the spectrum of interest and the corresponding noiseless template.

One must note that the velocity and temperature derivatives may have a cross term: if the velocity of the template and that of the spectrum of interest differ, then there will be an error in the determination of $\delta T$. We, therefore, compute the $\delta T$ after convergence of the velocity iterations of the LBL. While a small ($\ll 1$ line width) velocity offset does not impact temperature measurements, larger velocity residuals will lead to a mismatch between the temperature gradient spectrum and the spectrum being analyzed.

Considering that an optimal weighted sum is performed with weights proportional to the inverse of the square of uncertainties, one therefore has a mean temperature change for the domain of interest. This change in the disk-averaged temperature is designated as $dT_{\text{Temp}}$ and determined in a weighted sum through the following relation:

$$ d\text{Temp} = \frac{\sum \delta T(i) \sigma_T(i)^{-2}}{\sum \sigma_T(i)^{-2}}. \qquad (5) $$

The uncertainty of a weighted sum becomes:

$$ \sigma_{dT_{\text{Temp}}} = \frac{1}{\sqrt{\sum \sigma_T(i)^{-2}}}. \qquad (6) $$

We keep the framework of the LBL and opt for a simple mixture model (see Appendix B in É. Artigau et al. 2022) to compile all per-line measurements into a final temperature estimate.

One key element in the above framework is to determine the derivative as a function of temperature, which requires a library of relatively high-SNR (typically >100) template spectra covering a wide range of effective temperatures. Here, we explore temperature changes for G, K, and M main-sequence dwarfs (6000 K and cooler), but the technique could be expanded to any temperature range and/or surface gravity provided that one has access to a similar library of template spectra.

Measurements of $dT_{\text{Temp}}$ for G, K, and M main-sequence stars can be performed within LBL by requesting a given combination of temperature gradient and rotational broadening. The $dT_{\text{Temp}}$ measurements, along with their associated uncertainties, are computed simultaneously with velocity measurements and returned in a single table.

### 3. Template Construction

In this demonstration of the differential temperature measurement, we focus on data sets obtained with the SPIRou (J.-F. Donati et al. 2020) and HARPS (F. Pepe et al. 2002;

M. Mayor et al. 2003) instruments to showcase the performance of this activity indicator at both optical and infrared wavelengths. The M dwarf templates from SPIRou that were used are observations from the SPIRou legacy survey (P. Fouqué et al. 2023; C. Moutou et al. 2023), while K and G stars templates are archival observations obtained in the context of exoplanet transit surveys. For SPIRou, these templates are byproducts of the APERO (N. J. Cook et al. 2022) data reduction software. For HARPS, we performed a bulk archival query from the ESO archive (N. Delmotte et al. 2006) for all objects within 50 pc having >100 visits. HARPS spectra were telluric-corrected with the correction routine available within the LBL package (É. Artigau et al. 2022), which nulls the correlation function of telluric absorbers by adjusting the optical depth of telluric absorption from models (J. L. Bertaux et al. 2014). This technique closely matches the mathematical framework by R. Allart et al. (2022) with the difference that one adjusts an optical density rather than an effective pressure and optical density.

Line profiles are, overall, affected by four main parameters: stellar temperature, metallicity, surface gravity, and stellar rotation. Since we are interested only in the temperature dependency of our templates, we selected only slow rotators from all HARPS and SPIRou templates. From the FWHM of the cross-correlation function of our templates, we estimate that all objects have a $v \sin i < 3$ km s$^{-1}$, but note that constraints on the $v \sin i$ significantly below the combined contribution of the instrument resolution (2.6 and 4.3 km s$^{-1}$ for HARPS and SPIRou) and intrinsic stellar line widths is challenging. Selecting objects with a uniform metallicity distribution is more challenging. Metallicity measurements are notoriously challenging for M dwarfs (e.g., F. Jahandar et al. 2024 and references therein) and all references overlapped with only a handful of targets within our template list. Rather than using a very heterogeneous set of references, some with contradictory values for the same targets, we opted to use the metallicity measurements from the GAIA DR3 catalog (Vizier catalog I/355/PARAMP; Gaia Collaboration 2022) that provides, among other parameters, an estimate of the temperature, surface gravity, and metallicity for most of our template stars. Stars for the current analysis were kept if they met the following criteria:

1. Have parameters defined in the Gaia DR3 catalog.

2. A temperature between 2500 and 6500 K.

3. A metallicity comprised between -0.5 and +0.5 dex.

4. A surface gravity consistent with the main sequence. From the GAIA catalog, log g decreases monotonically from 4.9 to 4.2 respectively for $T_{\text{eff}}$ of 3000–6000. We kept objects within 0.2 dex of this relation.

For all the analysis that follows, the temperatures used were drawn from this catalog. While not necessarily the most precise, it has the advantage of being homogeneous in its biases for the entire sample.

Table A1 compiles all stars used as templates, the instrument (s) concerned, the effective temperature, and the systemic velocity used for spectra registering. The temperature distribution for the SPIRou and HARPS templates is shown in Figure 1. The M dwarf domain (2800–4000 K) is well covered for SPIRou, while mid-K dwarfs to G stars are very well covered for HARPS with a sparser sampling of mid- to late-M dwarfs (<3200 K). Each template is constructed from the median combination of all observations of each target.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Temperature distribution of the stars used for the SPIRou and HARPS templates. A total of 44 stars were used for SPIRou, 105 with HARPS, among which eight are in common.</figcaption></figure>

All spectra were registered to a common systemic velocity and high-passed with a median filter on a scale of 100 km s⁻¹. As illustrated in Figure 2 for sample domains, one sees a coherent evolution of the depth of spectral features with the temperature of the templates. For a given spectral pixel (see Figure 3), the relative flux depth can be fitted with a low-order polynomial as a function of temperature.

For the 2500–6000 K domain, we empirically found that a fourth-order polynomial with temperature is sufficient to fit our data properly. The analytical derivative of this polynomial gives its derivative at any given temperature. This provides the ∂A₀/∂T needed in Equations (3) through (5).

Figure 3 illustrates the ∂A₀/∂T spectra over the 3000–6000 K domain for temperature steps of 500 K for the sample domains shown in Figures 2 and 3 for the SPIRou spectra. The spectrum is expressed in K⁻¹, consisting of a fractional variation in the depth of a feature expressed as a fraction of the local low-pass filtered pseudo-continuum. The amplitude of the temperature gradient spectra is generally very low, typically less than a few 10⁻⁴ K⁻¹ (i.e., at most a few tens of percent in fractional line depth change for a 1000 K change in temperature). Figure 4 illustrates the mean templates and the corresponding temperature gradients for a representative domain within the SPIRou domain.

The scatter to the fit shown in Figure 3 informs on the precision one needs on the stellar temperatures. Given the rms to the polynomial fit with temperature (see Figure 3), considering that typical spectral gradients are on the order of 10⁻⁴ K⁻¹, the required precision for the stellar temperatures such that is a minor contributor to the dispersion is 100–300 K (i.e., with the gradients considered, a 100 K error contributes, in quadrature, 0.01 to the scatter). Since the typical error in temperature for the GAIA DR3 values, compared to detailed spectroscopic analysis, is 110 K (R. Andrae et al. 2023), we conclude that this is a relatively minor contributor to the overall scatter in the relations shown earlier. Going into a detailed, target-per-target T<sub>eff</sub> determination through spectroscopic analysis would not, in the present context, change our results.

We note that the concept of measuring precise differential temperature using per-line depth measurements is not new. D. F. Gray & H. L. Johanson (1991) and D. F. Gray (1994) used line-depth ratios to measure temperatures of late-F to

early-K stars down to differential precision of ~10 K by selecting iron and vanadium, which are highly sensitive to temperature changes within the r photometric bandpass. They used an empirical line-depth ratio relation with color, as a proxy of temperature, to derive a polynomial relation that is similar conceptually to the generalized per-spectral-pixel polynomial relations presented here (see in particular Figures 1 and 2 in D. F. Gray & H. L. Johanson 1991). One key difference with this early work is that we exploit the full information content of the spectra and use a purely empirical determination of the temperature changes in spectra without prior knowledge of the dominant absorber at each wavelength. This is particularly important for M dwarf spectra that are dominated by molecular bands that are notoriously challenging to reproduce, at the spectral resolution considered here, with stellar models (e.g., E. Artigau et al. 2018b; F. Jahandar et al. 2024).

## 4. Application to pRV Data Sets

We used the framework described above and applied it to three representative stars with the goal of constraining their rotation period through the measured temperature variations induced by stellar heterogeneities:

1. AU Mic:<fnref n="8" /> a very active M dwarf at infrared wavelengths.

2. Barnard's star:<fnref n="9" /> an RV-quiet M dwarf at infrared wavelengths.

3. ε Eridani:<fnref n="10" /> an active K dwarf at optical wavelengths.

### 4.1. Active M Dwarf, Near-infrared: AU Mic

AU Mic is among the best examples of active M dwarfs in the solar vicinity, with a dust disk (P. Kalas et al. 2004). The star has two known transiting planets (P. Plavchan et al. 2020; E. Martioli et al. 2021; G. M. Szabó et al. 2021, 2022), with the detection of a third planet through transit-timing variations (J. M. Wittrock et al. 2023) and a tentative detection of a fourth one from radial velocity monitoring (J.-F. Donati et al. 2023; J. M. Wittrock et al. 2023). The data set analyzed here is the same as presented in J.-F. Donati et al. (2023), and has been reduced with APERO version 0.7.284. In our analysis we used the temperature gradient for 3500 K convolved with a rotational kernel of 9 km s⁻¹. The resulting dTemp time series is shown in Figures 5 and 6. The dTemp time series shows a strong peak-to-peak variation at the level of ~40 K.

The time series was fitted with a Gaussian process (GP) using a quasiperiodic covariance function (R. D. Haywood et al. 2014; V. Rajpaul et al. 2015; S. Stock et al. 2023):

$$k(\tau) = \sigma_{\text{GP}}^2 \exp\left(-\frac{\tau^2}{2l^2} - \Gamma \sin^2\left(\frac{\pi\tau}{P_{\text{tot}}}\right)\right) + \sigma_{\text{jit}}^2 \quad (7)$$

As defined in S. Stock et al. (2023), the parameters are described here accordingly:

1. σ<sub>GP</sub> is the amplitude of the GP component given in Kelvin.

<footnote n="8">⁸ CFHT run IDs: 19AD97, 19AF26, 19AH28, 19AP42, 19BD97, 20AH93, 20AP42, 20BP42, 21AP42, 21BP42, 22AP42, 22BP45, 23AP45, 23BP45.</footnote>

<footnote n="9">⁹ CFHT run IDs: 18BE97, 19AP40, 19BP40, 20AP40, 20BP40, 21AP40, 21BP40, 22AP40, 23AQ57, 23BQ57, 24AQ57.</footnote>

<footnote n="10">¹⁰ ESO program ID: 0104.C-0863(A), 072.C-0488(E), 073.C-0784(B), 074.C-0012(A), 074.C-0012(B), 076.C-0878(A), 077.C-0530(A), 078.C-0833(A), 079.C-0681(A), 192.C-0852(A), 60.A-9036(A).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Sample spectroscopic features in the HARPS (left-hand panel) and SPIRou (right-hand panel) for sample domains. The temperature dependency of the spectral element highlighted with a dashed line is shown in Figure 3. The line centered at 632.935 nm is a Ni I feature while the feature at 1638.50 nm is tentatively attributed to an Fe II line.</figcaption></figure>

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

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

<figure><img src="image_9.png" /><figcaption>**Figure 3.** Temperature dependency of the spectral element shown in Figure 2 with the fourth-order polynomial fit for that pixel. Residuals to the fit have an rms of 0.087 and 0.036 for SPIRou and HARPS, respectively.</figcaption></figure>

2. *l* is the correlation length scale, in units of days, of the GP squared-exponential component, and represents the coherence timescale of repeating structures in the periodic signal.

3. Γ defines the dimensionless relative weight between the GP sine-squared component and squared-exponential component.

4. *P*<sub>rot</sub> is the period of the GP sine-squared component in days.

5. τ is the time lag, in days, between observations for which we are determining the covariance.

<figure><img src="image_10.png" /><figcaption>**Figure 4.** (Top panel) Sample average spectra at 500 K steps for SPIRou. The corresponding gradient spectra are shown on the bottom row. One can see that lines that change in depth in the top plot correspond to strong temperature gradient features in the bottom plot.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 5.** Sample sequence within the AU Mic *dTemp* time series. A Gaussian-process fit to the data is shown as an overplot (top panel), with the residual to that fit in the bottom panel. One sees the envelope of the signal decaying in strength through the first half of the sequence as well as the notable presence of a strong second harmonic signal (secondary bump appearing around JD = 2459030). The small-scale magnetic field measurements form J.-F. Donati et al. (2023) are shown for comparison. *dTemp* and ⟨*B*⟩ show a remarkable agreement in their overall behavior, with a strong anticorrelation.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 6.** Same as Figure 5 over the entire time series.</figcaption></figure>

6. $\sigma_{jit}$ is the excess noise (jitter term; in Kelvin) to fully account for the variance in our data. A value of $\sigma_{jit}$ converging to zero implies that the variance in the data is fully described by the GP model and the propagated uncertainties in the data.

The recovered GP model parameters and the corresponding constraints on the physical parameters of the evolution of the *dTemp* time series of AU Mic are given in Table 1. The time series of *dTemp* and the small-scale magnetic field (⟨*B*⟩) of AU Mic taken from J.-F. Donati et al. (2023) are shown in Figure 5 (subset of 130 days for clarity), with the residual to the

corresponding GP fit. The full time series is shown in Figure 6. The derived rotation period of 4.8620$^{+0.0046}_{-0.0029}$ days is in agreement with the 4.8636$^{+0.0021}_{-0.0019}$ days for the posterior value obtained with a GP fit on ⟨*B*⟩ in J.-F. Donati et al. (2023; Table 2 therein).

As shown in Figure 7, the *dTemp* time series is remarkably anticorrelated with the small-scale magnetic field, with a Pearson-*r* correlation between *dTemp* and ⟨*B*⟩ of -0.924, indicating that the two draw their origin in the same underlying physical processes. Stronger magnetic fields are linked to an apparent decrease in the bulk temperature of AU Mic.
<table><caption>**Table 1**
Priors and Posteriors on the GP Model of AU Mic</caption><thead><tr><th>Parameter</th><th>Prior</th><th>Posterior</th></tr></thead><tbody><tr><td>σ<sub>GP</sub> (K)</td><td>LU(0.15, 1500)</td><td>15.9<sup>+1.7</sup><sub>-2.7</sub></td></tr><tr><td>l (days)</td><td>LU(10, 1000)</td><td>124<sup>+20</sup><sub>-14</sub></td></tr><tr><td>Γ</td><td>LU(0.1, 10)</td><td>1.38<sup>+0.36</sup><sub>-0.34</sub></td></tr><tr><td>P<sub>rot</sub> (days)</td><td>LU(4.8, 4.9)</td><td>4.8620<sup>+0.0046</sup><sub>-0.0029</sub></td></tr><tr><td>σ<sub>jitter</sub> (K)</td><td>LU(0.15, 1500)</td><td>1.64<sup>+0.10</sup><sub>-0.12</sub></td></tr></tbody></table>

Note. LU(a, b) refers to the the log-uniform (“Jeffreys”) distribution.

This anticorrelation provides an effective means of determining the small-scale magnetic field $\langle B \rangle$ of any M dwarfs through a dTemp measurement derived from a high-resolution spectrograph. From Figure 7, we derive the following empirical calibration for AU Mic between $\langle B \rangle$ and dTemp (respectively, in kGauss and Kelvin), accounting for the uncertainties in both $\langle B \rangle$ and dTemp (i.e., using an orthogonal distance regression; P. J. Brown et al. 1990):

$$ \langle B \rangle = (-0.0086 \pm 0.0003)d\text{Temp} + (2.733 \pm 0.004). \quad (8) $$

The above relation has a typical dispersion of 0.047 kG compared to a median error on $\langle B \rangle$ of 0.027 kG. While the anticorrelation is remarkable in the case of AU Mic, it remains to be seen if it is also present for stars that are much less active and at different temperatures; we defer this question for future work. Considering that $\langle B \rangle$ is an excellent proxy for RV jitter (R. D. Haywood et al. 2022) and that this quantity displays, at the very least in the case of AU Mic, a strong anticorrelation of dTemp, we expect that dTemp will be similarly good for RV activity filtering.

## 4.2. Inactive M Dwarf, Near-infrared: Barnard's Star

Barnard's star (E. E. Barnard 1918) is the second nearest stellar system to the Sun at a distance of 1.8282 ± 0.00013 pc (Gaia Collaboration et al. 2021), and has been the topic of much interest in exoplanet searches. Despite numerous attempts at finding planets, all discovery claims to date have been challenged, from early astrometric claims (P. van de Kamp 1963) to the more recent claim by I. Ribas et al. (2018) of a planet with a 233 days orbital period.

Using the data set presented in É. Artigau et al. (2022), we compute the temperature variations of Barnard's star, which is notable for its low level of activity, and correspondingly low level of RV jitter. We did not convolve the derivative spectrum because Barnard's star is a very slow rotator with an equatorial velocity of <100 m s⁻¹. Contrary to AU Mic, one expects a much more modest activity signal. Using the same GP kernel as in Section 4.1 we constrain the time covariance of the time series. Figure 8 shows the evolution of Barnard's star temperature over nearly 4 yr. Its temperature fluctuates by ±2 K peak-to-peak with a semiperiodic behavior. The GP fit parameters are provided in Table 2.

The rotation period recovered here (153.3<sup>+3.1</sup><sub>-2.4</sub> days) is close to the published value inferred through activity-rotation relations (148.6 days; A. Suárez Mascareño et al. 2015) and spectropolarimetry (137 ± 5 days; P. Fouqué et al. 2023). The value is mildly discrepant with the photometric measurements of 130.4 days obtained by G. F. Benedict et al. (1998) with Hubble Space Telescope fine-guider photometry; however, we note that this reference does not report uncertainties.

z

<figure><img src="image_13.png" /><figcaption>**Figure 7.** $\langle B \rangle$ as a function of dTemp for the time series shown in Figure 5. An increase in temperature is associated with a decrease in the small-scale magnetic field. This is qualitatively understood as a decrease in the number of active regions on the observed hemisphere, leading to a smaller level of activity, fewer observed stellar spots, and hence a hotter average temperature. The correlation between the quantities (–0.924) is excellent, suggesting that they trace the same underlying physical processes. The 1σ envelope of Equation (8) is highlighted in purple.</figcaption></figure>

Furthermore, M dwarfs may very well experience differential rotation; therefore, features inducing photometric modulations in G. F. Benedict et al. (1998) may not occur at the same latitude as structures that led to other rotation measurements. The uncertainties on GP parameters from dTemp are comparable to those derived in P. Fouqué et al. (2023) from longitudinal magnetic field measurements.

## 4.3. $\epsilon$ Eridani

$\epsilon$ Eridani is the second closest K dwarf to the Sun after α Centauri B and the closest isolated one. It is located 3.220 ± 0.004 pc (Gaia Collaboration et al. 2021) from the Sun and has had its radius and temperature, 0.74 ± 0.01 R<sub>☉</sub> and 5039 ± 126 K, determined through interferometric measurements (E. K. Baines & J. T. Armstrong 2012). RV measurements point to a planetary-mass companion of a ~7 yr orbit but pRV detection of planets around $\epsilon$ Eridani is challenging due to the relatively high level of activity of this 0.4–0.8 Gyr old star (E. E. Mamajek & L. A. Hillenbrand 2008) and orbital parameters (planet mass and period) are inconsistent between various recent studies of the system (A. P. Hatzes et al. 2000; J. Llop-Sayson et al. 2021; F. Feng et al. 2023). This star shows activity-induced jitter at the ~10 m s⁻¹ level (A. F. Lanza et al. 2014; M. J. Giguere et al. 2016; A. Hempelmann et al. 2016) which is comparable to the RV semiamplitude of $\epsilon$ Eridani b (F. Feng et al. 2023). While activity can be filtered to 2–3 m s⁻¹ (M. J. Giguere et al. 2016; R. M. Roettenbacher et al. 2022), this has only been demonstrated on timescales commensurate with the rotation
<figure><img src="image_14.png" /><figcaption>**Figure 8.** *dTemp* time series for Barnard's star (top panel). As expected from the modest level of activity of this thick disk star, a much milder variability is observed than for AU Mic, but nonetheless one notes the rotation-induced variability. The residual to the GP fit (bottom panel) shows sub-K rms, consistent with the $\sigma_{jit} = 0.493^{+0.037}_{-0.040}$ K reported in Table 2.</figcaption></figure>

<table><caption>**Table 2**
Barnard's Star Priors, Posteriors, and Constraints on the GP Model</caption><thead><tr><th>Parameter</th><th>Prior</th><th>Posterior</th></tr></thead><tbody><tr><td>σ<sub>GP</sub> (K)</td><td>LU(0.015, 150)</td><td>3.10<sup>+0.70</sup><sub>-1.11</sub></td></tr><tr><td>l (days)</td><td>LU(300, 1000)</td><td>309.0<sup>+6.8</sup><sub>-15.1</sub></td></tr><tr><td>Γ</td><td>LU(0.1, 10)</td><td>0.83<sup>+0.32</sup><sub>-0.61</sub></td></tr><tr><td>P<sub>rot</sub> (days)</td><td>LU(100, 170)</td><td>153.3<sup>+3.1</sup><sub>-2.4</sub></td></tr><tr><td>σ<sub>jit</sub> (K)</td><td>LU(0.015, 150)</td><td>0.493<sup>+0.037</sup><sub>-0.040</sub></td></tr></tbody></table>

period and not on a timescale similar to the activity cycle of K dwarfs.

ε Eridani has been observed with HARPS 741 times over the 2003–2019 period, but we focus our analysis on the observations taken between 2019 October 4 and 2019 December 27. Over this period of 83 days, it was observed on 62 nights with, on most occasions, three spectra per night (190 spectra total). This is the densest monitoring of this star and it is the best data set to constrain the evolution of its *dTemp* time series. The template used in the *dTemp* computation was constructed from all observations but only the 2019 sequence is sufficiently dense for proper retrieval of the rotation period. For our analysis, we used a temperature gradient spectrum corresponding to 5000 K convolved with a 3 km s⁻¹ rotation kernel (M. J. Giguere et al. 2016). The *dTemp* time series of ε Eridani shown in Figure 9 allow for an accurate GP retrieval of the rotation period; the properties of the fit are reported in Table 3.

Notable among main sequences stars in the solar neighborhood, ε Eridani is sufficiently close to directly resolve its active regions through interferometry, and the correspondence between spatially resolved starspot and RV jitter measured with an optical pRV spectrograph (EXPRES) is discussed in R. M. Roettenbacher et al. (2022). Performing a joint analysis of the *dTemp* with interferometric data would be even more compelling because it would lift degeneracies between filling factors and starspot temperature contrast.

The GP fit to the ε Eridani time series differs from those derived through optical and near-infrared high-resolution spectropolarimetry and TESS photometry (P. Petit et al. 2021). The rotation period derived is $P_{eq} = 10.77 \pm 0.06$ days at the equator (see Figure 9 therein) with a differential rotation

$d/\Omega = 0.11 \pm 0.01 \text{ rad day}^{-1}$. The equatorial velocity is discrepancy can be readily explained if the dominant features creating the modulation shown in Figure 9 are not equatorial. P. Petit et al. (2021) derive a polar rotation of 13.3 days, and using Equation (6) therein we determine that the temperature differences we have measured occur at a latitude of $53.0^{+8.6^{\circ}}_{-7.6^{\circ}}$. This is at higher latitudes than for the Sun (typically within latitudes of $\pm 30^{\circ}$) but reported for some Kepler quarters for the Sun-like star HD 173701 through asteroseismology (M. Bazot et al. 2018) and similar to values derived for the K1V star HD 106225 (A. P. Hatzes 1998).

## 5. Transit-induced Temperature Changes

While variations of the temperature over a full rotation were the focus of the above demonstration data sets, an interesting case of short-term apparent disk-averaged temperature change appears during a transit. As the stellar disk is occulted by a companion (an exoplanet in the case presented here), the averaged temperature will change slightly. Masking the area close to the limb will remove cooler-appearing regions (from the line-of-sight geometry at that position on the stellar surface) and the disk-averaged temperature will therefore increase. Conversely, masking the center of the star, the hottest-appearing area on the star (in the simplistic case without active regions), one will see an apparent cooling of the star. Since the transiting planet HD 189733b has been observed several times with SPIRou,<fnref n="11" /> we predicted the amplitude of this unreported transit effect in our data. The *dTemp* was measured with a gradient spectrum at 5000 K and assuming a 3 km s⁻¹ v sin i (E. Cristo et al. 2024). The analysis accounts for the Rossiter–McLaughlin (RM) effect (D. B. McLaughlin 1924; R. A. Rossiter 1924); as mentioned in Section 2, the *dTemp* is computed after determining the average velocity.

### 5.1. Constructing a Toy Model of dTemp

While the limb of a star has a temperature that is cooler than the disk-averaged stellar spectral energy distribution, its exact

<footnote n="11">¹¹ CFHT run IDs: 18BE99, 19AP40, 20AC01, 21BC16, 23AQ57, and 23BQ57.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 9.** *dTemp* time series for $\epsilon$ Eridani at optical wavelength with HARPS. Rotation-induced variability is detected in the second half of the data set and is readily described by a semiperiodic GP fit (top). The residuals to the GP fit are at the 0.1 K level (bottom). In both time series, individual points have been binned on a nightly basis.</figcaption></figure>

<table><caption>**Table 3**
$\epsilon$ Eridani Priors, Posteriors, and Constraints on the GP Model</caption><thead><tr><th>Parameter</th><th>Prior</th><th>Posterior</th></tr></thead><tbody><tr><td>σ<sub>GP</sub> (K)</td><td>L<sub>U</sub>(0.008, 80)</td><td>1.08<sup>+0.25</sup><sub>-0.44</sub></td></tr><tr><td>l (days)</td><td>L<sub>U</sub>(10, 1000)</td><td>18.6<sup>+3.6</sup><sub>-3.7</sub></td></tr><tr><td>Γ</td><td>L<sub>U</sub>(0.1, 10)</td><td>0.74<sup>+0.30</sup><sub>-0.43</sub></td></tr><tr><td>P<sub>rot</sub> (days)</td><td>L<sub>U</sub>(8, 20)</td><td>12.27<sup>+0.31</sup><sub>-0.36</sub></td></tr><tr><td>σ<sub>jit</sub> (K)</td><td>L<sub>U</sub>(0.008, 80)</td><td>0.023<sup>+0.012</sup><sub>-0.030</sub></td></tr></tbody></table>

spectrum is not identical to that of a cooler star. To construct a model prediction of the *dTemp* equivalent of the RM effect, we retrieved the PHOENIX Goettigen atmosphere models (T.-O. Husser et al. 2013) with their full radiation field. We constructed a 3D cube with wavelength and spatial position on the stellar disk and constructed a time series of the total spectrum through transit. We assumed no stellar rotation (i.e., stellar spectra are not RV-shifted at the corresponding limb), which is a valid assumption for a slowly rotating star. This leads to a gray light curve as the sum of the flux over the entire disk. It represents the classical transit light curve, characterized by the standard central dip centered at midtransit due to limb darkening of the star (see panel (1) of Figure 10). Furthermore, we allow for unocculted regions on the stellar disk, where a given fraction of the star, outside of the transit chord, is covered with regions at different temperatures from the rest of the photosphere. These unocculted regions could be either cooler (spots) or hotter (plages) than the bulk of the photosphere. For unocculted regions, we assumed the same limb darkening as for the bulk of the photosphere. Figure 10 illustrates the geometry of the toy model used here.

To test the accuracy of our methodology, we produced a model light curve within the TESS bandpasses and compared it with the TESS phase-folded light curve for HD 189733b integrated over the 3D cube, the agreement is excellent, with a <0.5 mmag accuracy on the light curve for a ~2.8% transit depth.

<figure><img src="image_16.png" /><figcaption>**Figure 10.** Cartoon representation of the toy model used to predict the amplitude of the *dTemp* effect in the HD 189733b transit data. The impact parameter of 0.69 corresponds to that of the system; the relative radii of the host star and planet correspond to literature values. Unocculted spots covering 2.8% of the disk are shown (filling factor from D. K. Sing et al. 2011).</figcaption></figure>

In the same framework, we derived a residual transit spectrum by taking the ratio of the out-of-transit spectrum to that at any transit position. This residual was projected onto the derivative of the spectrum with temperature, determined within the same framework as described in Section 3, but constructed with models at temperatures immediately above and below that of the temperature of interest (in 100 K steps here). The predicted *dTemp* measurement through transit is shown in Figure 11, panel (1), with impact parameters of 0.0, 0.5, 0.69 (nominal value), to 0.9. The impact of the light source effect is highlighted in Figure 11, panel (3). Unocculted regions occupying 2.8% (D. K. Sing et al. 2011) of the photosphere and a temperature difference ranging from 0 to 2000 K. The impact on the *dTemp* signal is significant, with an evolution of the morphology of the curve; the central dip strengthens with cooler unocculted regions. The impact on the transit depth, on the contrary, is minimal and would be identical to a change in planet-to-star radius ratio.
<figure><img src="image_17.png" /><figcaption>**Figure 11.** (1) Model *dT*Temp for different impact parameters for the HD 189733 system and the corresponding broadband light curves (2) integrated over the SPIRou bandpass. Low-impact-parameter transits would lead to first an apparent increase of temperature (planet occulting the limb) and a decrease of the mean temperature when the planet occults regions close to the disk center. Grazing transits lead to a positive excursion only. (3) For the nominal impact parameter (*b* = 0.69) and a 2.8% spot filling factor (D. K. Sing et al. 2011), we show the predicted *dT*Temp signal through transit (3) as well as the corresponding broadband light curve (4).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 12.** *dT*Temp time series of four HD189733b transit observed with SPIRou and phase-binned values. One notes the ~1 K increase in the apparent temperature of the host star close to ingress and egress as cooler regions of the stellar disk are masked by the planet. The beige domain corresponds to the window between first and fourth contact for an impact parameter of 0.69.</figcaption></figure>

As detailed in Section 6, the gradient spectrum derived from models only partially matches the empirically derived one. The toy model constructed here assumes that the atmosphere models are correct in a differential way. Even if all line positions are not accurate, changing the filling factor during transit (e.g., light source effects in Figure 11) should be correctly predicted, to a first order, provided that gradient spectra used to compute the *dT*Temp values are derived consistently with the spectra used to construct these effects.

5.1.1. HD 189733 Transit dTemp Time Series

Figure 12 overplots five transits of HD 189733b in the temperature space. Each time series is subtracted from a baseline between 1.8 and 1.0 half-transit durations on either

side of the event. Phase-binned points are overplotted and one sees a temperature excess just after ingress and just before egress. One sees an increase in *dT*Temp of ~1 K just after ingress and just before degree, with a ~0.5 K decrease below the baseline *dT*Temp value at midtransit. Qualitatively, this is understood as an increase in the disk-averaged temperature when the limb is occulted and a decrease in disk-averaged temperature, at midtransit, when central regions of the disk are occulted. This *dT*Temp variation through transit is analogous to the RM, but in temperature rather than velocity space.

While this qualitative behavior is reproduced by the model described in Section 5.1, no simple model reproduced quantitatively the observed *dT*Temp signal. Figure 13 illustrates the binned *dT*Temp signal in comparison with the predicted signature varying the impact parameter of HD 189733b. The
<figure><img src="image_19.png" /><figcaption>**Figure 13.** Binned *dTemp* measurements through HD 189733b transits compared with a change in the impact parameter. The ~0.5 K dip in *dTemp* compared to out-of-transit values as well as the maximum in *dTemp* at a position of ±0.8 stellar radii would suggest an impact parameter of *b* ~ 0.5, which is significantly smaller than the literature value (*b* = 0.69). Here, we assume that the stellar photosphere is uniform (i.e., no spots, no plages).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 14.** Binned *dTemp* measurements through HD 189733b transit compared with a change in the temperature contrast of spots covering 2.8% of the photosphere. We assume that spots are outside of the transit chord (as in Figure 10). The dip below the baseline value at midtransit suggests a starspot contrast of ~1000 K, but this model strongly underpredicts the rise at the start and end of the transit. We note that having very cold spots (−2000 K) or no spots leads to a very similar signature; very cold spots contribute little to the overall spectrum and have no effects on the *dTemp* signal.</figcaption></figure>

overall behavior is best reproduced with an impact parameter of *b* ~ 0.5, which leads to a lower midtransit minimum than the literature impact parameter. Considering that the impact parameter is very well constrained by space-based photometry and accurate pRV work (e.g., E. Cristo et al. 2024), we do not claim that literature values of the impact parameter are erroneous, but other effects must be at play. Figure 14 illustrates the impact of unocculted spots with varying temperatures. The central “dip” in the *dTemp* signature is well reproduced with a contrast of ~1000 K, which is also consistent with the temperature contrast derived by D. K. Sing et al. (2011). While this unocculted spot model matches the

central part of the transit, it fails to reproduce the amplitude of the *dTemp* increase just after ingress and just before egress. No unocculted spot models reproduce the ~1 K increase in *dTemp* close to ingress and egress and their maxima happens too close to midtransit compared to the SPIRou time series. No simple combination of parameters (spot contrast and distribution) that assumed a transit geometry consistent with literature values has been found to fit the observed *dTemp* sequence. Understanding the underlying reasons for this mismatch is beyond the scope of the current work. The HD 189733 system has been observed with several high-resolution spectrographs (e.g., SOPHIE, HARPS, HIRES, ESPRESSO; A. F. Lanza et al. 2011; P. W. Cauley et al. 2017; E. Cristo et al. 2024). The joint analysis of the *dTemp* effect in these data sets, along with the SPIRou observations presented here, may shed light on the origin of this discrepancy between simplistic models and the observed transit *dTemp* effect. It remains to be seen if a more complex geometry (e.g., combining spots and plages) or if inconsistencies in stellar atmosphere models should be invoked.

Spot and plage-crossing events should also manifest themselves as transit-induced temperature changes. These are readily detectable in photometric and spectro-photometric data but may be harder to identify in ground-based, high-resolution spectroscopic transit observations. While no such event is seen in the data set in hand, starspot crossing events should manifest themselves as a positive *dTemp* excursion. The peculiar TOI-3884 system with a persistent starspot transited by a super-Neptune (J. M. Almenara et al. 2022; J. E. Libby-Roberts et al. 2023) would be a particularly compelling case study. Obtaining simultaneous photometric and *dTemp* measurements would provide additional constraints regarding spot filling ratio and temperature contrast. Overall, the *dTemp* through transit provides a new tool to quantify the transit light source effect (B. V. Rackham et al. 2018), which is emerging as a significant issue in JWST transit spectroscopy (e.g., O. Lim et al. 2023).

## 6. Models versus Empirical *dTemp* Spectra

The usage of empirical temperature sequences, as described in Section 2, implies that one must assemble a sizable library of archival spectra of stars at known temperatures. Ideally, a *dTemp* spectrum could be derived with models and applied to the time series, alleviating this requirement. We explored this approach with spectra using PHOENIX (T.-O. Husser et al. 2013) spectral library. S. Shahaf & B. Zackay (2023) also considered the spectral variation of spectra with wavelength. They explore the representation of radial velocity shifts in the context of Taylor expansions relative to various quantities, including velocity (the RV shift) and temperature (equivalent to $\partial A_0/\partial T$ here).

A *dTemp* spectrum was obtained by applying the high-pass filter to two spectra at increasing temperatures and applying the same method as for spectral templates. Figure 15 shows the temperature gradient at 3500 K for both theoretical and empirical spectra. The overall behavior of the temperature gradient is reproduced, as well as strong features such as the K I doublet. Numerous small features are not reproduced correctly, which is consistent with the fact that the details of M dwarf spectra are, generally, poorly reproduced at high spectral resolution, and one-to-one identification of lines can be challenging (E. Artigau et al. 2018b; F. Jahandar et al. 2024).
<figure><img src="image_21.png" /><figcaption>**Figure 15.** Empirical and model temperature gradient spectra for the whole SPIRou domain (top panel), centered on the J-band K1 doublet (middle panel) and a sample domain within H band for $T_{\text{eff}} = 3500$ K. The overall behavior of the temperature gradient is well reproduced and most of the deep lines are reproduced by models. One notes that numerous small features, particularly in the H band, are not accurately predicted by models.</figcaption></figure>

Deriving a *dT*Temp time series with a template that only partially correlates with actual temperature gradients will lead to a lower apparent temperature change, as well as a loss in the significance of the detection. We verified this by using the theoretical temperature gradient shown in Figure 15 on the AU Mic time series, with the results shown in Figure 16. The overall behavior is preserved and one can retrieve the rotation period with the same accuracy as with the empirical temperature gradient. The GP fit on the time series derived from model gradients gives a rotation period of $4.8635_{-0.0029}^{+0.0024}$ days, a scale length of $160_{-32}^{+27}$ days, and shape factor $\Gamma = 1.67_{-0.44}^{+0.43}$, all being consistent with values reported in Table 1. The derived amplitude $\sigma_{\text{GP}} = 7.81_{-1.32}^{+0.96}$ K is a factor ~2 smaller. This difference in amplitude arises from the fact that the model *dT*Temp spectrum only partially traces the empirical effect of a temperature change.

We determined the *dT*Temp with model spectra for HD189733; Figure 17 shows the *dT*Temp signal with both empirical and model temperature gradient spectral. The effect described in Section 5.1.1 disappears completely and one gets a much flatter time series. The median uncertainties are reduced from 0.97 to 0.39 K as models overpredict temperature gradients at 5000 K, but this comes with a true loss in the accuracy of the measurements. Overall, empirical gradients are strongly preferred for this type of analysis and *dT*Temp values retrieved in this framework with model-based gradients should be viewed, at best, with caution.

## 7. Conclusions and Future Prospects

We highlighted with sample data sets the interest of a new stellar activity indicator, the *dT*Temp value, in time series of high-resolution spectroscopy, in particular in the context of pRV exoplanet searches. The indicators provide a measurement of the differential temperature of the star in a time series and, for bright targets, can reach uncertainties well below 1 K per measurement. The activity indicator relies on a library of

spectral templates of known temperature, but this library only needs to be determined once for a given combination of spectral type and wavelength domain (see the Appendix for access to the templates). Therefore, new instruments do not have to reobserve a large set of templates to measure this indicator.

The framework opens a new window in stellar astrophysics, allowing for very subtle changes in disk-averaged temperatures to be measured, at the sub-K level. If one expresses this as the equivalent of a photometric change, a 1 K temperature change on a 4000 K blackbody (roughly the temperature of an M0) corresponds to a photometric variation of 1 mmag (assuming a $T_{\text{eff}}^4$ dependency of flux). Obtaining sub-mmag photometry simultaneously with pRV observations over baselines of months to years is virtually impossible from the ground; sub-mmag ground-based photometry from the ground has only been demonstrated within nights (e.g., G. Stefansson et al. 2017). With *dT*Temp measurements, the time sampling is contemporaneous with pRV measurement because both are drawn from the same data set. Future work will explore the equivalent of the $f \cdot f'$ (S. Aigrain et al. 2012) framework but using temperature rather than the flux to correct for activity jitter in pRV sequences.

It may be nontrivial to determine the underlying physical processes traced by a given measured *dT*Temp signal. In the case of AU Mic, the notable correspondence between small scale $\langle B \rangle$ and *dT*Temp points toward the magnetic origin of the temperature variations, but these may be purely geometric in the case of HD 189733. Stellar pulsations also probably lead to detectable *dT*Temp for G and K dwarfs, a case that has not been investigated here. One more subtle effect may be unrelated physical causes, such as variations in the Zeeman splitting from the magnetic field, that lead to a modification of line profiles with a nonzero cross term with the temperature derivative spectrum. Changes in chemical abundances, leading to a change in line depth, may also lead to an erroneous *dT*Temp interpretation. Another plausible mechanism for *dT*Temp
<figure><img src="image_22.png" /><figcaption>**Figure 16.** *dTemp* time series with empirical and model temperature gradient spectra for the same timespan as in Figure 5. The signature of the activity is recovered but with an amplitude reduced by a factor ~2 while the periodicity and correlation length of the GP fits are consistent.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 17.** Same as for Figure 12, but comparing the effect of using empirical and model-based gradient spectra. The overprediction of temperature gradients in models leads to artificially small uncertainties but the signature of the effect described in Section 5.1.1 disappears.</figcaption></figure>

signatures is through planet-star interaction that has yet to be explored. Close-in planets (e.g., A. Strugarek et al. 2014, 2019) interact with the large-scale magnetic field of the host star, as the planet rotates inside the magnetosphere. This second-order effect in $\Delta T_{\text{eff}}$ may be detectable as a modulation at the orbital period of the planet. Detecting a periodicity peak in *dTemp* at the orbital period of a tentative radial-velocity detection or its higher harmonics should therefore be seen with caution. Although it may point to a false positive of stellar activity masquerading as a planetary signal, it may also be indicative of planet-star interaction.

The framework presented here highlights the usage of the *dTemp* framework in the context of high-resolution spectroscopy at optical and infrared wavelengths. Conceptually, the framework should work at any spectral resolution, provided that the SNR is sufficiently high. Transit-spectroscopy time

series with JWST are particularly appealing for a low-resolution application of the *dTemp* framework. As pointed out earlier, the transit light source effect is a significant limitation to transit spectroscopy and must be accounted for in the interpretation of transit-spectroscopy data, in particular at the bluer end of the JWST domain. O. Lim et al. (2023), R. Doyon (2024), and C. Cadieux et al. (2024a) illustrate that for mid- to late-M dwarfs, the domain shortward of ~1.0 µm becomes dominated by these effects. Out-of-transit *dTemp* effects in JWST data could point to rotation-induced temperature modulation and, when combined with raw photometry, inform on the filling factor of active regions. In-transit *dTemp* effects, analogous to that of HD 189733b shown in Figure 12, would inform on light source effects and provide a strong prior for the models used in the interpretation of transit-spectroscopy data. One limiting factor in constructing a *dTemp* framework for JWST transit data is the lack of a large spectral library, uniformly stepped in effective stellar temperatures, obtained with each observation mode. One can expect such a library to become gradually available as more JWST transit data becomes public in the next few years.

## Acknowledgments

This work was based on observations obtained at the Canada-France-Hawai'i Telescope which is operated from the summit of Maunakea by the National Research Council of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawai'i. The observations at the Canada-France-Hawai'i Telescope were performed with care and respect from the summit of Maunakea which is a significant cultural and historic site. This work is partly supported by the Natural Science and Engineering Research Council of Canada, the Fonds Québécois de Recherche (Nature et Technologie) the Trottier Family Foundation through the Institute for Research on Exoplanets.
E.A., C.C., N.J.C., R.D., M.C., and L.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

J.F.D., C.M., and P.C. acknowledge funding from the European Research Council (ERC) under the H2020 Research & Innovation Program (grant agreement #740651 New-Worlds).

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).

Facilities: CFHT (SPIRou), ESO:3.6m (HARPS), ESO Science Archive.

### Software:

1. Astropy; a community-developed core Python package for Astronomy; The Astropy Collaboration et al. (2013, 2018).

2. george; S. Ambikasaran et al. (2016).

3. IPython; F. Perez & B. E. Granger (2007).

4. matplotlib, a Python library for publication quality graphics; J. D. Hunter (2007).

5. NumPy; C. R. Harris et al. (2020).

6. radvel; B. J. Fulton et al. (2018).

7. SciPy; SciPy 1.0 Contributors et al. (2020).

8. tqdm; C. Da Costa-Luis et al. (2021).

9. During the analysis work, E.A. and N.C. made extensive use of the Dace platform (N. Buchschacher et al. 2015).

## Appendix
### Data Availability

The template derivatives used here are provided for the HARPS and SPIRou domains in the LBL webpage.<fnref n="12" /> Table A1 provides the list of objects used as HARPS and/or SPIRou

<table><caption>**Table A1**
First Entries in the Full Table, Available in Electronic Form, of the Templates used for the Current Work</caption><thead><tr><th>SIMBAD<sup>a</sup><br>Identifier</th><th>Instrument(s)</th><th>T<sub>eff</sub><br>(K)</th><th>Fe/H<br>(dex)</th><th>log g<br>(dex)</th><th>V<sub>sys</sub><br>(m s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>G 157-77</td><td>SPIRou</td><td>2863</td><td>-0.47</td><td>5.10</td><td>-41066</td></tr><tr><td>PM J22114+4059</td><td>SPIRou</td><td>2919</td><td>-0.04</td><td>5.07</td><td>-16835</td></tr><tr><td>G 141-36</td><td>SPIRou</td><td>2929</td><td>-0.21</td><td>5.11</td><td>-34544</td></tr><tr><td>LP 791-18</td><td>SPIRou</td><td>2930</td><td>-0.23</td><td>5.11</td><td>14545</td></tr><tr><td>L 788-37</td><td>HARPS</td><td>3061</td><td>-0.13</td><td>5.05</td><td>-24292</td></tr><tr><td>G 112-50</td><td>SPIRou</td><td>3075</td><td>-0.05</td><td>5.02</td><td>38704</td></tr><tr><td>G 36-24</td><td>SPIRou</td><td>3089</td><td>-0.31</td><td>4.86</td><td>-6489</td></tr><tr><td>Ross 128</td><td>HARPS, SPIRou</td><td>3093</td><td>-0.20</td><td>5.06</td><td>-31111</td></tr><tr><td>TOI-1452</td><td>SPIRou</td><td>3170</td><td>-0.03</td><td>4.98</td><td>-34471</td></tr><tr><td>BD-15 6290</td><td>HARPS, SPIRou</td><td>3185</td><td>-0.00</td><td>4.78</td><td>-1490</td></tr></tbody></table>

**Notes.** A total of 105 have been used with HARPS and 44 with SPIRou, with eight in common between the two instruments.

templates, the $T_{\text{eff}}$ used for the polynomial fit in Section 3, the surface gravity, and the systemic velocity for template allignment. We plan to upgrade these templates with full coverage of the 400–2450 nm domain with archival data from instruments covering the far-red domain. The dTemp indicator is computed for temperatures ranging from 3000 to 6000 K in steps of 500 K and for rotational broadenings ranging from 0 to 19 km s⁻¹ in steps of 1 km s⁻¹. The LBL framework now accepts data from several instruments (i.e., SPIRou, NIRPS, SOPHIE, HARPS, HARPS-N, CARMENES-VIS, ESPRESSO, MAROON-X), and other instruments may be added upon reasonable request and pending resource availability.

## ORCID iDs

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

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

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

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

Laurie Dauplaise https://orcid.org/0009-0004-2993-7849

Luc Arnold https://orcid.org/0000-0002-0111-1234

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

Paul Cristofari https://orcid.org/0000-0003-4019-0630

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

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

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

Pierre Larue https://orcid.org/0009-0005-1139-3502

Romain Allart https://orcid.org/0000-0002-1199-9759

## References

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

Allart, R., Lovis, C., Faria, J., et al. 2022, A&A, 666, A196

Almenara, J. M., Bonfils, X., Forveille, T., et al. 2022, A&A, 667, L11

Al Moulla, K., Dumusque, X., Cretignier, M., Zhao, Y., & Valenti, J. A. 2022, A&A, 664, A34

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

Andrae, R., Fouesneau, M., Sordo, R., et al. 2023, A&A, 674, A27

Artigau, É., Bouchy, F., Doyon, R., et al. 2024, arXiv:2406.08304

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

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

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

Baines, E. K., & Armstrong, J. T. 2012, ApJ, 744, 138

Baliunas, S. L., Donahue, R. A., Soon, W. H., et al. 1995, ApJ, 438, 269

Barnard, E. E. 1918, AAS, 3, 191

Bazot, M., Campante, T. L., Chaplin, W. J., et al. 2012, A&A, 544, A106

Bazot, M., Nielsen, M. B., Mary, D., et al. 2018, A&A, 619, L9

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

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

Boisse, I., Bouchy, F., Hébrard, G., et al. 2011, A&A, 528, A4

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

Brewer, B. J., Bedding, T. R., Kjeldsen, H., & Stello, D. 2007, ApJ, 654, 551

Brown, P. J., Fuller, W. A. & American Mathematical Society, Institute of Mathematical Statistics, & Society for Industrial and Applied Mathematics 1990, in Statistical Analysis of Measurement Error Models and Applications: Proc. AMS-IMS-SIAM (Providence, RI: American Mathematical Society)

Buchschacher, N., Ségransan, D., Udry, S., & Díaz, R. 2015, in ASP Conf. Proc. 495, Astronomical Data Analysis Software an Systems XXIV, ed. A. R. Taylor & E. Rosolowsky (San Francisco, CA: ASP), 7

Cadieux, C., Doyon, R., MacDonald, R. J., et al. 2024a, ApJL, 970, L2

Cadieux, C., Doyon, R., Plotnykov, M., et al. 2022, AJ, 164, 96

Cadieux, C., Plotnykov, M., Doyon, R., et al. 2024b, ApJL, 960, L3

Carmona, A., Delfosse, X., Bellotti, S., et al. 2023, A&A, 674, A110

Cauley, P. W., Redfield, S., & Jensen, A. G. 2017, AJ, 153, 217

<footnote>a M. Wenger et al. (2000).</footnote>

<footnote>(This table is available in its entirety in machine-readable form in the online article.)</footnote>

<footnote n="12">¹² https://lbl.exoplanets.ca</footnote>
Cegla, H. M., Watson, C. A., Shelyag, S., et al. 2018, ApJ, 866, 55

Cook, N. J., Artigau, É., Doyon, R., et al. 2022, PASP, 134, 114509

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

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

Cristo, E., Esparza Borges, E., Santos, N. C., et al. 2024, A&A, 682, A28

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

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

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

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

Doyon, R. 2024, arXiv:2403.12617

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

Feng, F., Butler, R. P., Vogt, S. S., Holden, B., & Rui, Y. 2023, MNRAS, 525, 607

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

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

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

Gaia Collaboration 2022, yCat, I/355

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

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

Giguere, M. J., Fischer, D. A., Zhang, C. X. Y., et al. 2016, ApJ, 824, 150

Gomes Da Silva, J., Santos, N. C., Bonfils, X., et al. 2012, A&A, 541, A9

Gray, D. F. 1994, PASP, 106, 1248

Gray, D. F., & Johanson, H. L. 1991, PASP, 103, 439

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

Hatzes, A. P. 1998, A&A, 330, 541

Hatzes, A. P., Cochran, W. D., McArthur, B., et al. 2000, ApJL, 544, L145

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

Haywood, R. D., Milbourne, T. W., Saar, S. H., et al. 2022, ApJ, 935, 6

Hempelmann, A., Mittag, M., Gonzalez-Perez, J. N., et al. 2016, A&A, 586, A14

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

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

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

Jahandar, F., Doyon, R., Artigau, E., et al. 2024, ApJ, 966, 56

Kalas, P., Liu, M. C., & Matthews, B. C. 2004, Sci, 303, 1990

Kallinger, T., De Ridder, J., Hekker, S., et al. 2014, A&A, 570, A41

Kürster, M., Endl, M., Rouesnel, F., et al. 2003, A&A, 403, 1077

Lanza, A. F., Boisse, I., Bouchy, F., Bonomo, A. S., & Moutou, C. 2011, A&A, 533, A44

Lanza, A. F., Das Chagas, M. L., & De Medeiros, J. R. 2014, A&A, 564, A50

Libby-Roberts, J. E., Schutte, M., Hebb, L., et al. 2023, AJ, 165, 249

Lim, O., Benneke, B., Doyon, R., et al. 2023, ApJL, 955, L22

Llop-Sayson, J., Wang, J. J., Ruffio, J.-B., et al. 2021, AJ, 162, 181

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

Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264

Martioli, E., Hébrard, G., Correia, A. C. M., Laskar, J., & Lecavelier Des Etangs, A. 2021, A&A, 649, A177

Mayor, M., Pepe, F., Queloz, D., et al. 2003, Msgr, 114, 20

McLaughlin, D. B. 1924, ApJ, 60, 22

Meunier, N., Desort, M., & Lagrange, A.-M. 2010a, A&A, 512, A39

Meunier, N., Lagrange, A.-M., & Desort, M. 2010b, A&A, 519, A66

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

Noyes, R. W., Hartmann, L. W., Baliunas, S. L., Duncan, D. K., & Vaughan, A. H. 1984, ApJ, 279, 763

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

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

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

Petit, P., Folsom, C. P., Donati, J.-F., et al. 2021, A&A, 648, A55

Plavchan, P., Barclay, T., Gagné, J., et al. 2020, Natur, 582, 497

Queloz, D., Bouchy, F., Moutou, C., et al. 2009, A&A, 506, 303

Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279

Rackham, B. V., Apai, D., & Giampapa, M. S. 2018, ApJ, 853, 122

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

Rajpaul, V., Aigrain, S., & Roberts, S. 2016, MNRAS, 456, L6

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

Roettenbacher, R. M., Cabot, S. H. C., Fischer, D. A., et al. 2022, AJ, 163, 19

Rossiter, R. A. 1924, ApJ, 60, 15

Schöfer, P., Jeffers, S. V., Reiners, A., et al. 2019, A&A, 623, A44

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

Shahaf, S., & Zackay, B. 2023, MNRAS, 525, 6223

Sing, D. K., Pont, F., Aigrain, S., et al. 2011, MNRAS, 416, 1443

Stefansson, G., Mahadevan, S., Hebb, L., et al. 2017, ApJ, 848, 9

Stock, S., Kemmer, J., Kossakowski, D., et al. 2023, A&A, 674, A108

Strugarek, A., Brun, A. S., Donati, J.-F., Moutou, C., & Réville, V. 2019, ApJ, 881, 136

Strugarek, A., Brun, A. S., Matt, S. P., & Réville, V. 2014, ApJ, 795, 86

Suárez Mascareño, A., Rebolo, R., González Hernández, J. I., & Esposito, M. 2015, MNRAS, 452, 2745

Szabó, G. M., Gandolfi, D., Brandeker, A., et al. 2021, A&A, 654, A159

Szabó, G. M., Garai, Z., Brandeker, A., et al. 2022, A&A, 659, L7

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

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

van de Kamp, P. 1963, AJ, 68, 515

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

Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9

Wittrock, J. M., Plavchan, P. P., Cale, B. L., et al. 2023, AJ, 166, 232

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