Measuring precise radial velocities on individual spectral lines. IV. Stellar activity correlation with line formation temperature

Context. Radial velocities (RVs) of stars contain both the Doppler reflex motion of potential planetary companions and the drowning and sometimes imitating effect of stellar activity. To separate the two, previous efforts have sought for proxys which only trace the activity signals, yet the sub-meter-per-second floor required for the detection of Earth-like planets remains difficult to break. Aims. In this work, we analyze a sample of 12 G- to early M-type stars in order to investigate the feasibility of detecting a differential effect of stellar activity with photospheric depth, as traced by the spectral line-forming temperature, for observations with different sampling and noise levels. Methods. We computed the average line formation temperature for each point in the observed wavelength grids using the spectral synthesis code PySME. The final line selection was curated to exclude blended and poorly synthesized lines. We thereafter computed the convective blueshift (CB) of the line cores of our master spectra (composed of the stacked individual spectra for each star). Finally, we extract RV time series for certain intervals of formation temperature using a template-matching approach. Results. We find the CB to follow a linear relation with the formation temperature of the line cores, and the CB slope to be steeper with increasing effective temperature. For the RV time series derived for different intervals of formation temperature, we find the RVs of line parts formed at higher temperatures, close to the spectral continuum, to be generally correlated with the S index, and RVs of line parts formed at cooler temperatures, close to the spectral line cores, to be generally anti-correlated, especially for stars with low noise levels and significant variations over their magnetic cycles.

⋆ Based on observations made with the HARPS instrument on the ESO 3.6m Telescope at the La Silla Observatory, Chile.
Previous efforts using stellar activity indicators, such as the full-width at half-maximum (FWHM) and the bisector inverse slope (BIS) of the cross-correlation function (CCF) (CCF; Queloz et al. 2001;Pepe et al. 2002) or the Mount Wilson S index extracted from the Ca II H and K lines (Baliunas et al. 1983), have not been able to sufficiently mitigate the activity signal since these indicators are seldom fully correlated with the RVs (Haywood et al. 2016).Recent studies (e.g., Cretignier et al. 2021Cretignier et al. , 2022) ) have instead attempted to implement mitigating techniques at the spectral level, before the CCF is applied to extract the RVs, since it is now well established that stellar activity (primarily the inhibition of convective blueshift (CB) by magnetically active regions), affects spectral lines to different extents (e.g., Reiners et al. 2016;Dumusque 2018).
In this study, we investigate how stellar activity affects individual spectral lines parts formed at various temperature regimes in the stellar photosphere.We apply the methodology of Al Moulla et al. (2022), which studied this dependence for the Sun and α Cen B (HD 128621, which is also included in our sample to be analyzed on the long-term), and found that RVs of line parts formed in hotter parts of the photosphere, close to the spectral continuum, are more strongly correlated with the RV contribution of the inhibition of CB compared to line parts formed at cooler temperatures.Here, we extend the work to a sample of 12 late-type (G to early M) dwarfs, in order to investigate whether similar effects can be seen in other spectral types.
In Sect.2, we describe the stellar sample and the available spectral data.In Sect.3, we describe the method used to compute the RV of selected line parts.In Sect.4, we analyze the correlations between stellar activity and formation temperature-binned RV signals.In Sect.5, we discuss our results and conclude.

Spectral pretreatment
The stellar sample consists of 12 G-to early M-type stars observed with the HARPS 1 spectrograph, mounted on the ESO 3.6m Telescope at the La Silla Observatory, Chile, which covers a spectral range from 378 to 691 nm and provides an instrumental resolution of R = 115 000.The targets have been frequently observed at relatively high signal-to-noise ratio (S/N) and exhibit clear S index variations indicative of magnetic activity cycles over the observed time span (except for HD 211038, which appears to have been quiet in recent years).They are listed in Table 1, with their observational data, detailing the number of spectra per target, the median S/N, the span in S index after smoothing with a 200-day rolling average (to be less affected by outliers and rotational modulation of the S index), and the span in log R ′ HK converted from the smoothed S index (Noyes et al. 1984).log R ′ HK is only a re-scaling of the S index; however, it also corrects for the photospheric contribution to the emission in the Ca II H and K lines using the B−V color of each star (see Table 2), permitting a better comparison between spectral types.We list the span in log R ′ HK in Table 1 to compare the relative variations within our stellar sample; however, we use the S index for the remainder of our analysis since we only consider correlations with the activity proxies and not their absolute values.The listed numbers of spectra represent the number of retained and nightly stacked spectra after processing with 1 HARPS stands for High Accuracy Radial velocity Planet Searcher (HARPS).
YARARA (Cretignier et al. 2021(Cretignier et al. , 2023)), during which a qualitybased rejection could be made, and hence the numbers could differ compared to archival data.YARARA performs corrections of stellar activity, telluric absorption, and instrumental systematics using mainly principal component analysis (PCA) applied on the residual of the spectral data matrix (residual flux versus time).
YARARA also provides a master spectrum computed as the inverse-variance weighted average of all individual spectra at each sampled wavelength point.The absolute RV shift of the master spectrum might lead to an arbitrary offset when template matching with the individual spectra (for example, if the RV of the master spectrum is high, low-activity spectra will appear to be slightly more blueshifted relative to the master compared to a master constructed using only the quietest observations).However, since we are only interested in RV variations, this offset should not be of much importance.More importantly, the YARARA corrections work better at higher S/N, which is why we opted to create master spectra using all available observations.These masters serve as the references both to compare the syntheses (Sect.3.1) with true observations, and to perform template-matching when extracting the temperature-binned RV time series (Sect.3.4).After the corrections, we only re-injected the corrected stellar activity of the individual spectra in order to isolate its effect at the spectral level.

Stellar parameters
The adopted stellar parameters are listed in Table 2.The spectral types and B−V colors were retrieved from SIMBAD2 (Wenger et al. 2000).The effective temperatures, T eff , surface gravities, log g, metallicities, [Fe/H], projected rotational velocities, v sin i, masses, M, and radii, R, were retrieved from the references in the table foot.For HD 21693, the mass and the radius could not be found in the literature and were instead computed from the luminosity, L, using the following formulas from Valenti & Fischer (2005): and where quantities with the subscript ⊙ denote solar values, and σ M , σ R , and σ L are the uncertainties of the mass, radius, and luminosity, respectively.It should be noted that the masses and radii were solely used for the computation of the gravitational redshift (see Sect. 3.3) and not needed for the subsequent computations of the temperature-dependent RVs.

Methods
One of the primary ways of investigating the effects of stellar activity at the spectral level is to extract its differential effect on

Target
Sp. type Notes.
The columns indicate the target name, spectral type, B−V color, effective temperature, surface gravity, metallicity, projected rotational velocity, mass, radius, luminosity, and references.
The stars are ordered by effective temperature.For HD 21693, the mass and radius were computed using the effective temperature, surface gravity, and luminosity, as is described in the text.spectral lines as is seen in line-by-line (LBL) RVs (Dumusque 2018;Cretignier et al. 2020;Artigau et al. 2022;Lafarga et al. 2023).In order to further trace the origin of activity, one can, instead of considering entire spectral lines (or spectral regions bounded by two local maxima as in the case of M dwarfs where spectral lines are everywhere), consider various line or CCF parts (Al Moulla et al. 2022;Siegel et al. 2022) formed at certain temperature ranges in the stellar photosphere.Contrary to the LBL method, whose feasibility only relies on observational data, this part-by-part (PBP) method also requires a priori knowledge about where certain spectral features are formed in the photosphere, which can be retrieved from synthetic spectra and thereafter mapped onto cross-matched lines in the observed spectra.
For each target we used a unique line list from the VALD3 3,4 database (Piskunov et al. 1995;Kupka et al. 2000;Ryabchikova et al. 2015), which were queried for a wavelength interval of 3500-7000 Å, slightly larger than the range of the HARPS spectrograph.In the queries we used the stellar parameters in Table 1 and a normalized line depth threshold of 0.01 for all targets except the M dwarf (Gl 229), for which we used a threshold of 0.15 in order to produce line lists with less than a hundred thousand lines (the maximal output of VALD).We used a grid of MARCS (Gustafsson et al. 2008) model atmospheres interpolated for the same stellar parameters as for the queries.The syntheses were broadened with their adopted v sin i (see Table 2), and micro-and macroturbulences adopted from Valenti & Fischer (2005).The microturbulence was fixed to for all stars, and the macroturbulence was set to In Eq. ( 4), we have changed the sign of the T eff -dependent term compared to Valenti & Fischer (2005), since their equation has an apparent typo.This is evident from the plotted relation in their Fig. 3 and was also pointed out by Takeda & UeNo (2017).Although the relations in Valenti & Fischer (2005) were derived for FGK stars, we assume they approximately extend to early M dwarfs as well.Furthermore, minor discrepancies in turbulence broadening are not a limiting factor for our analysis since the subsequent rejection of poorly synthesized lines is relative for each star, and hence if all lines are equally overor under-broadened the final selection remains almost identical.
The spectra were finally convolved with a Gaussian instrumental profile to resemble the HARPS resolution of R = 115 000.For each target, a second synthesis excluding rotational, turbulent, and instrumental broadening was computed for the identification of line blends not otherwise seen.

Average formation temperature
The average formation temperature, denoted as T 1/2 , was then calculated at each sampled wavelength point of these synthetic

Convective blueshift
For each observed line in the master spectra (see Sect. 2.1), defined as a local minimum bounded by two local maxima, we recalibrated the central wavelength by fitting a second-order polynomial to the central-most seven points.Once the formation temperature spectra were computed, we interpolated their value at the line center, denoted as T 1/2,c , of each observed line.
We then identified the species of the observed lines using VALD line lists.The cross-match with VALD was done by querying a VALD line list according to the stellar parameters of each star, and then finding all VALD lines that fall within the left and right borders of the observed lines.For observed lines that only have one possible match in the VALD list, the VALD line was taken as the correct identification.Using these uniquely identified lines, an average wavelength dispersion (caused by uncorrected shifts and/or discrepancy in the wavelength solution) was measured.
For observed lines that have multiple possible matches, the VALD line with the smallest Euclidean norm of the wavelength and depth differences was adopted if the wavelength was within the dispersion measured on uniquely identified lines.Some observed lines were expected to remain unidentified due to having poorly characterized atomic properties, being tellurics or noisy patterns falsely flagged as stellar lines.We thereafter computed the Doppler shift between the observed central wavelengths of identified lines and the restframe wavelengths from VALD, known as the convective blueshift (CB).This shift originates due to the motion of the convective surface layer in late-type stars, where ascending granules emit more light and cover a larger area than the descending intergranular lanes, causing a net blueshift (Dravins et al. 1981).In order to calibrate the CB to an absolute value, we also subtracted a constant gravitational redshift for each star, where G is the gravitational constant and c is the speed of light in vacuum.We note that the uncertainty on the gravitational redshift is not propagated into the individual line CBs, as the gravitational redshift is simply a constant offset for each star that we consider in order to have most of our values be negative, and hence truly blueshifts.For the remainder of our analysis we only focus on the CB slopes and not the zero points.
After rejecting blended lines (as determined both from the master spectra themselves and from the syntheses void of any broadening other than thermal) and poorly synthesized lines (see Al Moulla et al. 2022 for the detailed applied criteria), we were left with our final selection of lines, presented in Fig. 1.We remark that, for most stars, the CB decreases (to larger absolute values) with the average formation temperature in the line core as a result of probing hotter photospheric regions where the convective motions are more turbulent, and hence producing larger blueshifts.For each star, we regressed a linear model and present the best-fitting slopes in Fig. 2, where the distribution of slopes appears to also follow a linear relation with the effective temperature.The negative slope can be explained by the larger extent of the convective zone of later-type stars, where turnover timescales are longer, and thus the turbulent motions expected to be slower.The direction of this trend is also in agreement with the findings of Liebing et al. (2021), who fitted a scaling factor to the solar relation between the CB and the line depth for 810 F-M stars and found the scaling factor to increase with effective temperature, indicating a steeper slope for hotter stars.
A106, page 4 of 15 Fig. 2. Slope of the CB with respect to the average formation temperature of the line cores, as measured from the linear fits seen in Fig. 1.The slope shows a decrease with the stellar effective temperature.
From our limited sample, it remains unclear whether the CB slope plateaus to zero for the coolest stars, or becomes positive.Given the distribution of points in Fig. 1, the linear fits for Gl 229 (the M dwarf) and HD 85512 appear to be poorly constrained.For these stars, the spectral syntheses had larger discrepancies with observations compared to the earlier-type stars, and thus less reliable derived formation temperatures.Another outlier to be remarked upon is HD 211038, which has a steeper slope compared to the stars in the vicinity of its effective temperature.HD 211038 is the only subgiant in our sample of otherwise main-sequence (MS) stars, and it has a significantly lower surface gravity, log g = 3.87±0.06,compared to all the other stars.The steeper slope, indicating a stronger convection strength compared to MS stars of a similar effective temperature, is in agreement with the empirical relations found by Liebing et al. (2023) for evolved stars.

Temperature-binned RVs
For this study we decided to divide the entire range of formation temperatures per star into intervals of equal size, denoted as ∆T 1/2 (N, n), where N is the number of intervals and n = [1, N] the interval index going from cooler to hotter temperatures.We considered the N = 1 (equivalent to no temperature binning, i.e., conventional LBL) and N = 4 cases, in order to have comparable results with Al Moulla et al. (2022).It should be noted that the bounds of the temperature intervals for different stars do not necessarily agree with each other, since the intervals were assessed relative to the line-forming temperature regime of each star, that is, the available range of temperatures as given by the spectral syntheses.
The PBP RV of a formation temperature interval was computed by template-matching segments of the individual observed lines, which happen to fall within the considered temperature interval with the same entire line in the master spectrum (Dumusque 2018).This was then repeated for each temperature interval, spectral line, and spectrum in the time series.After an iterative 4σ outlier rejection on the RMS of each interval and line, we computed a weighted average across all remaining lines, resulting in a single RV time series per temperature interval.We require a minimum of three points for the template-matching; however, three points alone will produce a large uncertainty for a given line and temperature interval, such that it will likely be sigma-clipped when computing the line-averaged RV per temperature interval.Furthermore, the temperature intervals will not necessarily be similar across the wavelength range.The primary reason for this is the influence of continuous absorption on the derived formation temperature (see Fig. 2 in Al Moulla et al. 2022).In the optical range of Sun-like stars, the continuous absorption is dominated by the negative hydrogen ion, which increases its opacity up to about 1 µm (Gray 2008).Therefore, for the HARPS wavelength range, bluer wavelengths are biased to hotter temperatures and redder wavelengths to cooler ones (because the increased opacity allows us to probe less of the photosphere).But it is exactly these kinds of nuanced differences that we aim to achieve with our analysis.
The N = 1 and N = 4 RV time series are presented in Figs.A.1-A.12, together with the nightly averaged S indices and their respective generalized Lomb-Scargle (GLS;Lomb 1976;Scargle 1982;Zechmeister & Kürster 2009) periodograms.The S indices were retrieved from DACE 5,6 , which in turn archives them from the HARPS Data Reduction Software (DRS).The HARPS DRS computes the S index in an equivalent way to the original Mount Wilson survey (see Lovis et al. 2011).

Activity correlations of temperature-binned RVs
We computed the Pearson correlation coefficient between the S index and each temperature-binned RV time series, which are shown in Fig. 3.We find that for most stars, the RVs of line parts formed in hotter temperatures tend to be more positively correlated, and the RVs of line parts formed in cooler temperatures tend to be more negatively correlated.Furthermore, this relation appears to be amplified with the significance of the signal, as is shown both by the median S/N of the stellar spectra for each star and by the power of the largest periodogram peak for the S index relative to its 0.01% false alarm probability (FAP) level.The color-coded variable in Fig. 3 was selected as a quantification of how statistically significant the S index variations are.We believe that this is better than stating a potentially inaccurate estimate of the S index semi-amplitude since the cycles are not fully sampled.We find no clear dependence on the effective temperature.
These correlations are in agreement with the findings of Al Moulla et al. (2022; see their Figs.A.1 and B.3).However, while they attributed the measured signal at hotter temperature ranges to the suppression of CB by faculae, using simulations of the Sun with SOAP-GPU (Zhao & Dumusque 2023), no clear explanation has been established for the seemingly anti-correlated activity pattern observed at cooler temperatures.The S index, derived from the deep Ca II H and K lines, is known to exhibit chromospheric activity; however, the reason why activity signals at hotter to cooler formation temperatures transition from positively correlated, to weakly correlated, to negatively correlated with the S index requires a more fundamental understanding of the types of active regions and the velocity-altering processes to which chromospheric activity is connected.Forthcoming results by Siegel et al. (in prep.)show that RVs measured on the top versus the bottom parts of the CCF (approximately translating to spectral parts formed at hotter versus cooler temperature ranges) correlate with the photometric and convective influences of active regions, respectively.These results strongly indicate that spectral parts formed at different temperatures are primarily influenced by separate physical processes.Their connection to the S index remains the topic of future investigations.In this study, we showcase the consistency of the observed differential correlation for magnetically active stars.
For a few stars in our sample (HD 215152, HD 192310, and HD 128621), we find a striking, near-perfect (≤−0.8)anticorrelation for the ∆T 1/2 (4, 1) time series.Due to the largely deviating covariances of these RVs, we decided to detrend the ∆T 1/2 (1, 1) time series (which contains all the RV information) with proxies composed of the differences of various combinations of the 4-bin RVs.The reason for taking the residual of two RV time series as an activity index and not the RV time series themselves was to cancel out the potential presence of a planetary signal in the proxy, since we expect pure Keplerian signals to be of equal amplitude regardless of the bin configuration.We compared the ratio of the RV root mean square (RMS) before and after detrending to assess the improvement factor of the RV precision.We tried detrending with the S index alone as a reference case, before we detrended with one RV residual at a time, followed by a multi-linear regression (MLR) with all combinations of 4-bin residuals, and then a MLR with all 4-bin residuals and the S index combined.The RMS improvement factors for each star and detrending method are shown in Fig. 4. We find that for a majority of stars in our sample the MLR regression with all bin residuals and the S index is able to decrease the RMS to a greater extent than the S index alone, indicating that the 4-bin RV time series are likely tracing different activity contributions than those of the S index.

Discussion and conclusions
In Al Moulla et al. (2022), the authors demonstrated that for two MS stars for which we can reach extremely high S/N (>500 for the G2V Sun and the K1V α Cen B), it was possible to precisely measure the RV of different layers in the stellar atmosphere and to confidently measure the differential effect of stellar variability as a function of the line formation temperature, which is expected from our understanding of the RV activity signal (e.g., Meunier et al. 2017).In this paper, we further extend the analysis to 12 stars, with spectral types from G2V to M1V and observations that have a more reasonable S/N (∼100−300).We show again the differential effect of stellar variability as a function of depth inside the stellar photosphere, which demonstrates that this technique can be used to probe and correct for stellar activity on a wide variety of MS stars that are observed at a high, but not extreme, S/N.
To confidently measure the differential effect of stellar variability as a function of line formation temperature, there are many aspects to consider and justify.Carrying out the analysis for solar observations has the advantage of theoretical models, such as the spectral synthesis used to derive the formation temperature itself, usually being optimized for the Sun as a reference target.In addition, solar observations enjoy the benefit of rarely being noise limited due to their wealth of photons.Hence, to extend the analysis to solar-like stars requires the assumption of precise enough spectral modeling for other stellar parameters, as well as efficient high-resolution spectrographs that can reach high S/N for stellar observations.
In our work, we justify the use and accuracy of synthetic spectra mapped to high-resolution observations by verifying the consistency of derived relations, such as the net blueshift of line cores (see Sect. 3.3).Selecting unblended lines with relatively A106, page 6 of 15 small observation-synthesis flux residuals, we find a linear relation between the CB and core formation temperature analogous to the third signature of granulation observed when substituting temperature for line depth (Gray 2009;Reiners et al. 2016;Liebing et al. 2023).This also allows us to determine how our synthesis becomes successively less reliable for cooler stars as the relation becomes less constrained (cf.Fig. 1).
By nightly stacking and correcting our spectra for instrumental systematics using YARARA, we were able to extract temperature-binned RVs to an extent where the intrinsic variations still significantly exceed the white noise.Some of the selected stars do have confirmed planetary companions; however, planets with RV semi-amplitudes larger than a few meters per second were removed by YARARA.Any potentially remaining lower-mass planets would not affect our results significantly since they would contribute negligibly to the RV RMS and not be detectable with our limited methodology even after the improvement of the MLR.
We emphasize again the opposing effect of stellar activity on line parts formed at hotter versus cooler photospheric temperatures, as is seen in the correlation coefficients with the S index.The RVs originating from the coolest 25% appear to be the most sensitive to S index variations, and the negativity of the correlation alludes to this temperature regime being more sensitive to chromospheric activity differing from the dominating inhibition of CB at hotter temperatures.After detrending the RV time series computed on the entire range of formation temperatures with proxies composed of residuals of the time series computed on four temperature bins, we found that a MLR of all bin residuals together with the S index yielded an RMS improvement greater than detrending with the S index alone for a majority of the studied stars.On average, the improvement in RV RMS is on the order of 20% when using the MLR and S index.However, when analyzing the results of Fig. 4, it is still challenging to understand why in some cases the RV RMS is strongly reduced while in some other cases the improvement is minimal.For Gl 229, the RMS values after detrending are noticeably unchanged; the ratios are not exactly 1, but close enough to not be discernible.This is not too surprising, however, given the very low activity correlations of all its RV time series (see Fig. A.1).
From the results of this study, the residuals of the temperature-binned RVs are demonstrated as being excellent proxies with which to probe the stellar activity effect, mainly magnetic cycles, and to differentiate planetary long-period signals from stellar ones.However, since the efficiency of those proxies at correcting for the stellar signal changes from star to star, without a clear correlation with spectral type, it remains to be explored how the formalism of using line formation temperatures to extract PBP RVs can be further leveraged to attain RVs that target certain flavors of stellar activity or that are optimally mitigated for planet searches.

3
Fig.1.CB as a function of formation temperature at the line cores, T 1/2,c , for the stellar sample.The stars and their spectral types are labeled in the upper left corner of each subplot, and ordered by effective temperature, from left to right and top to bottom.For most stars, the CB shows a linear trend with T 1/2,c , indicating that lines formed at higher temperatures, i.e., deeper layers of the photosphere where convective motion is stronger, are subject to larger velocity shifts.The individual spectral lines are shown as points and the fitted linear trends are shown as dashed lines.

Fig. 3 .
Fig. 3. Pearson correlation, R, between the temperature-binned RVs and the S index.The stars are ordered by effective temperature.The marker size indicates the median S/N of the studied spectra of each star, and the marker color indicates the power, P, of the largest periodogram peak of the S index relative to its 0.01% FAP level.The color map inverts at 1, such that green points represent periodograms peaks above the 0.01% FAP, and red points represent peaks below it.

Fig. 4 .
Fig. 4. RMS improvement of the ∆T 1/2 (1, 1) RV time series after linear detrending with different proxies.The markers represent detrending with either the S index alone (circles with orange edges), different 4-bin residuals (convex polygons with black edges), a multi-linear combination of all 4-bin residuals (crosses), or a multi-linear combination of all 4-bin residuals and the S index (stars).The marker sizes and colors indicate the same properties as described in Fig. 3.

Table 1 .
Observational specifications relevant to our analysis.
Notes.The columns specify the target name, spectral type, number of spectra, median S/N after nightly stacking, span (peak to peak difference) in the S index after smoothing with a 200-day rolling average (see Figs. A.1-A.12), and span in log R ′ HK converted from the smoothed S index.

Table 2 .
Parameters of the stellar sample.