© The Authors 2026

1 Introduction

Since the discovery of 51 Pegasi b in 1995 (Mayor & Queloz 1995), the radial velocity (RV) method has allowed for the detection of over a thousand exoplanets, moving over time from large Jupiter-like planets to smaller rocky worlds orbiting in the temperate zone of their host stars (e.g. Hatzes et al. 2000; McArthur et al. 2004; Anglada-Escudé et al. 2016). This progress has been largely due to technical improvements on spectrographs, which have gradually pushed detection limits from signals in the tens of metres per second down to the m s⁻¹ domain, with instruments such as HARPS (Mayor et al. 2003), HARPS-N (Cosentino et al. 2012), CARMENES (Quirrenbach et al. 2016), and SPIRou (Donati et al. 2018). The newest spectrographs, for example, ESPRESSO (Pepe et al. 2010), with an instrumental precision reaching down to around 10 cm s⁻¹, represent a big leap towards the detection of even smaller Keplerian signals induced by Earth-like rocky planets (Figueira et al. 2025).

Despite these instrumental advancements, stellar magnetic activity remains a fundamental barrier to detecting low-amplitude RV signals. Surface inhomogeneities such as spots, faculae, and granulation can induce RV variations greater than 1 m s⁻¹, even in relatively quiet stars (e.g. Dumusque et al. 2012; Perger et al. 2017). These stellar activity effects distort spectral line shapes, potentially mimicking or obscuring planetary signals. The two main RV extraction methods are affected by these wavelength-dependent distortions, and the community has developed different strategies for each technique to disentangle planetary signals from instrumental or stellar activity variability.

In the cross-correlation function (CCF) method, the observed spectrum is cross-correlated with a weighted mask tailored to the stellar spectral type (Baranne et al. 1996). The resulting CCF represents an average stellar absorption line. For most stars, the CCF is well described by a Gaussian profile. The centroid yields the RV, while the contrast (CON) and the full width at half maximum ($\text{FWHM}=2\sqrt{\ln 2}\sigma$, where σ² is the variance of the Gaussian) are line-shape activity indicators. The bisector inverse slope (BIS), defined as the velocity difference between the upper (60–90%) and lower (10–40%) parts of the CCF, is used as a measurement of line asymmetry (Queloz et al. 2001). Several alternative approaches have been developed to extract activity information from the CCF, including bi-Gaussian fitting with asymmetric widths (Figueira et al. 2013), Fourier decomposition (Zhao & Tinney 2020), principal component analysis (PCA) of the auto-correlation function (Collier Cameron et al. 2021), and PCA on shape-driven CCFs orthogonalised with respect to the first derivative of a template (Klein et al. 2024).

In spectral-level methods, RVs are computed directly from the observed spectra using high signal-to-noise templates. The optimal wavelength shift that aligns an observation with the template is typically obtained through one of two main formalisms: (i) a numerical approach that performs a least-squares minimisation by exploring the velocity parameter space or (ii) an analytical approach based on the Bouchy et al. (2001) formalism, where the shift is inferred from a first-order Taylor expansion of the spectrum. The latter requires calculating the derivative of the flux with respect to wavelength from a nearly noise-free template. Either of these formalisms can be applied over different wavelength ranges, thus broadly dividing their practical application into two techniques. (1) Global template matching applies the shift measurement to the entire spectrum at once (Anglada-Escudé & Butler 2012; Zechmeister et al. 2018; Silva et al. 2022). Line-width variations are captured in the differential line width activity indicator (Zechmeister & Kürster 2009), and the difference in RV variability as a function of wavelength is measured with the chromatic index (Zechmeister et al. 2018; Baroch et al. 2020). (2) The line-by-line RV method allows us to measure the RV signal from each individual line with respect to the template, rather than only producing a single global RV value (Dumusque 2018; Artigau et al. 2022; Lafarga et al. 2023). This allows for detailed studies of differential line responses to activity, making it possible to select subsets of lines based on their line-by-line RV sensitivity to parameters such as line depth (Cretignier et al. 2022), formation temperature (Al Moulla et al. 2022), or spot-to-photosphere temperature contrast (Larue et al. 2025). More data-driven approaches, such as PCA-based line selection, have also been explored (Cretignier et al. 2023).

Several techniques have been developed to correct for stellar activity using either CCF-based or spectral-level diagnostics. Gaussian process (GP) regression has become a widely adopted method for modelling the quasi-periodic variability induced by stellar magnetic activity in RV time series (Haywood et al. 2014; Ambikasaran et al. 2015; Foreman-Mackey et al. 2017; Perger et al. 2021). More recently, multi-dimensional GP frameworks have been introduced, combining RVs with ancillary indicators such as photometry, CCF diagnostics, or chromospheric line indicators (e.g. Barragán et al. 2022; Delisle et al. 2022). Neural network (NN)-based methods have also gathered significant interest. Convolutional neural networks (CNNs) have been applied to solar CCFs (de Beurs et al. 2022) and to other input data, such as the spectral-shell representations (Cretignier et al. 2022; Zhao et al. 2024), though without explicit temporal modelling. Perger et al. (2023) extended CNNs to time series modelling of CCF indicators using physically motivated training data generated with the StarSim code (Herrero et al. 2016). Architectures such as autoencoders have also been proposed, for example, to distinguish real from apparent spectral-line shifts in simulated data (Liang et al. 2023).

In this work, we propose an alternative approach based on separating line shifts from line shape distortions with an orthogonal basis function whose derived activity indicators are used to train a time-aware NN to mitigate stellar activity in the RV data. We apply the methodology to two test stars: ϵ Eridani (ϵ Eri) and TZ Arietis (TZ Ari). These represent two distinct cases with different CCF shapes, so we aim to prove the capability of our method to derive activity indicators on CCFs with Gaussian and non-Gaussian shapes.

The paper is organised as follows. The test stars are presented in Sect. 2. In Sect. 3, we introduce the method to extract line-shape-indicator coefficients from the CCF, and in Sect. 4 we apply this theoretical framework to the test cases. In Sect. 5, we explain how we fed these indicators into a NN architecture for time series modelling (convolutional attention network) and compare its performance to a network agnostic to time information (fully connected network). We trained both networks on synthetic data produced with the StarSim code. In Sect. 6, we evaluate the best-performing network on the two test cases and study the improved sensitivity to planetary signals provided by the network’s stellar activity mitigation. In Sect. 7 we discuss the results, and in Sect. 8 we summarise the main findings and implications of the study.

2 Test cases ϵ Eri and TZ Ari

The decomposition framework of this study was validated on two test stars. These specific targets were chosen to showcase the framework’s ability to capture both Gaussian-like and non-Gaussian-like CCF profiles.

2.1 Target characteristics

2.1.1 ϵ Eri

ϵ Eri is a young (400–800 Myr; Mamajek & Hillenbrand 2008) K2 dwarf located at a distance of 3.22 pc, making it one of the closest known exoplanet hosts. With a rotation period of 11.2 d (Fröhlich 2007) and v sin i = 2.4 ± 0.5 km s⁻¹ (Valenti & Fischer 2005), it serves as a representative example of a moderately active star of this age. The system hosts a prominent debris disk (Greaves et al. 1998) and a Jupiter-mass planet, ϵ Eri b, with an orbital period of approximately seven years (Hatzes et al. 2000; Llop-Sayson et al. 2021). Giguere et al. (2016) estimated a stellar inclination of the rotation axis to be ${69.5}_{-7.6}^{+5.6}$ deg by modelling spot modulation in photometric and RV data. We summarise the relevant stellar parameters in Table 1.

For this study, we used 205 publicly available spectra from the HARPS spectrograph, on the ESO 3.6 m telescope at La Silla Observatory, Chile. HARPS covers the optical wavelength range from 0.38 to 0.69 µm at a resolving power of R ≈ 120,000. Our dataset covers a baseline of 88 days between October 5, 2019 (JD 2458762) and January 1, 2020 (JD 2458850). To mitigate short-term variations, we binned the observations nightly, resulting in 66 final data points.

2.1.2 TZ Ari

TZ Ari (Gl 83.1) is a nearby M5.0 dwarf located at 4.47 pc from Earth. Quirrenbach et al. (2022) reported a rotation period of 1.96 d based on spectroscopic indicators and ground-based photometric data from multiple facilities. This is confirmed by our analysis of TESS photometry (Ricker et al. 2015) from September to November 2023 (Appendix A). The rotation period combined with the upper limit on the projected rotational velocity v sin i < 2 km s⁻¹ (Reiners et al. 2018) implies an upper limit on the stellar inclination of 28°. Quirrenbach et al. (2022) also identified a Saturn-mass planet, TZ Ari b, with an orbital period of ${772.05}_{-1.84}^{+2.41}$ d, a minimum mass of 0.21 ± 0.02 M_Jup, and an orbital eccentricity of $e={0.46}_{-0.07}^{+0.06}$.

For this work, we used 92 spectra (after discarding three clear outliers as identified on the CCF activity indicators time series) obtained with the high-resolution CARMENES spectrograph on the 3.5 m telescope at the Calar Alto Observatory. We employed CARMENES visual channel, which covers the optical wavelength range from 0.52 to 0.96 µm at a resolving power of R ≈ 94,600. The observations span 92 nights between January 31, 2016 (JD 2457419), and November 29, 2019 (JD 2458817). The last observing season, starting on July 17, 2019 (JD 2458682), was conducted at a higher cadence to properly sample the stellar rotation period.

2.2 Spectroscopic data

The data for each target star were analysed with the CCF method using the raccoon pipeline (Lafarga et al. 2020). We employed masks tailored to the spectral type of each star (Table 1), adopting the detector pixel size as the velocity sampling step for the CCF computation. We computed the templates for each target star by time-averaging the flux-normalised CCFs. We normalised the flux by dividing each CCF by its baseline continuum and subtracting the result from unity, ensuring that the final profiles act as strictly positive, probability-like functions. The final templates are shown in Fig. 1. For the ϵ Eri template (red), a Gaussian function provides an accurate description of the shape. The TZ Ari template (blue) shows a Gaussian shape in the core of the CCF, but displays two humps at either side, which is typical of CCFs of cool M dwarfs measured at visible wavelengths (Lafarga et al. 2020). General characteristics of the RV, CON, FWHM and BIS time series for both targets are presented in the table of Appendix B.

Fig. 1 Time-averaged CCFs (templates) of ϵ Eri (red) and TZ Ari (blue) with the corresponding Gaussian model fit for ϵ Eri (orange) and the two-component Gaussian model fit for TZ Ari (cyan).

2.3 Basic data analysis

To identify periodicities in the time series, we employ the Generalised Lomb–Scargle (GLS, Zechmeister & Kürster 2009) periodogram. The significance of the detected peaks is assessed using false-alarm probability (FAP) levels derived via bootstrap resampling (Efron & Tibshirani 1985). We classify signals as tentative or significant if their power exceeds the 1% or 0.1% FAP levels, respectively.

We present the RV and activity indicator time series together with their GLS periodograms for ϵ Eri in Fig. 2. All time series exhibit significant power at the stellar rotation period and its first harmonic, both highlighted in purple. Similarly, Fig. 3 shows the corresponding data for TZ Ari. In this case, the RV, CON, and BIS time series display significant variability at the stellar rotation period as well as at its daily alias (1 − 1/1.96 = 1/2.04 d⁻¹). In addition, we observe long-term variability in the RVs associated with the orbital motion of TZ Ari b. Since GLS models are sinusoidal, the planet’s eccentricity shifts the highest peak slightly with respect to the literature value (Quirrenbach et al. 2022). The CON time series exhibits a downwards trend, introducing excess power at low frequencies in the periodogram, while the FWHM shows a significant signal at 40.17 d.

In Appendix B, we list the significant periods and amplitudes retrieved via a pre-whitening procedure. This process involves iteratively identifying the most significant periodogram peak, fitting the RVs with a sinusoidal function of the same frequency as that peak, and subtracting it from the time series until no power crossing the 10% FAP threshold remains.

Fig. 2 Time series (left) and GLS periodograms (right) of the RV and activity indicators derived from the CCFs for ϵ Eri. We mark in purple the stellar rotation period (11.2 d) and its first harmonic (5.6 d). The 0.1, 1, and 10% FAP levels are shown with dashed orange lines.

2.4 Synthetic data

2.4.1 StarSim modelling

We generated synthetic time series data using the StarSim¹ code (Herrero et al. 2016; Rosich et al. 2020; Gomes et al., in prep.) in order to model the test stars. StarSim is a physically motivated stellar activity simulator that produces synthetic photometric and spectroscopic observables by modelling the inhomogeneous surface of a rotating star. The stellar disk is discretised into surface elements corresponding to three distinct components: quiet photosphere, cool spots, and bright faculae. Each component is represented using a high-resolution synthetic spectrum from the PHOENIX library (Husser et al. 2013), selected to match the appropriate temperature, surface gravity, and metallicity. The library also includes spectra computed at different limb angles to account for centre-to-limb variations across the stellar disk. Rather than synthesising full high-resolution spectra from the surface maps, the simulator computes the CCF corresponding to each surface element directly, resulting in a considerable reduction of the computing time. This CCF computation process explicitly emulates the raccoon pipeline, employing the identical masks and velocity sampling steps applied to the observational data.

We constructed tailored stellar models for each target star using constraints on stellar parameters from the literature (see Table 1). To estimate the contrast temperature difference between photosphere and spot and the convective blueshift, we adopted the empirical relations from Herbst et al. (2021) and Liebing et al. (2021). However, we allowed for wide uniform priors on these parameters rather than fixing them to a single best-fit value. This ensured that the relationships do not impose overly strong constraints on the simulated values, providing the NN with a diverse training set required to generalise well.

We also simulated a large ensemble of stellar surface maps with active regions randomly distributed across the visible disk. Each active region is modelled as a circular spot, characterised by seven parameters: (1) appearance time, (2) lifetime, (3) colatitude, (4) longitude, (5) spot radius, (6) linear growth rate, and (7) linear decay rate. Surrounding each spot is a facular region, modelled as a concentric ring. The facular area is parametrized using a global facula-to-spot area ratio, Q, which is fixed for each simulation. As shown in previous studies (Lanza et al. 2003; Silva-Valio et al. 2010; Dumusque et al. 2014; Herrero et al. 2016), modelling circular spots surrounded by a facular corona successfully reproduces spot maps for high precision photometry of the Sun and other stars.

However, several studies combining Solar Dynamics Observatory (SDO) images and magnetograms of the Sun with RVs have demonstrated that the Sun is facula-dominated, possessing an extensive network of standalone faculae that persists even during solar minimum when spots are absent (Collier Cameron et al. 2021; Haywood et al. 2022). Consequently, an independent treatment of spots and faculae is typically preferred to accurately simulate such networks (Zhao et al. 2025). Given the high activity levels of our specific test targets (Sect. 2) and our objective to restrict the dimensionality of the parameter space, we keep the Q-factor approximation for this study. We note that this simplification is no longer present in the upcoming solar-benchmarked branch of the StarSim code (SunSim, Stucki et al., in prep.).

For simplicity, we assumed that each active region maintains a constant maximum size over most of its lifetime, with linear growth and decay timescales lasting a few days. Although more complex evolution laws for active regions could be implemented, the large number of simultaneously evolving regions, each with different lifetimes, locations, and appearance times, effectively captures the stochastic nature of surface inhomogeneities and approximates the effects of more realistic evolution scenarios.

We calibrate the simulations by adjusting the parameters of the active regions – specifically their number and size distributions – in order to reproduce the typical observed amplitudes and scatter in the spectroscopic time series. It is important to note that our objective was not to recreate the exact epoch-by-epoch RV pattern of the observations, nor did we perform a simultaneous fit (e.g. via Markov chain Monte Carlo, MCMC) of the spot contrast, convective blueshift, and active region geometries to the exact time series. Instead, our goal is to generate a robust training dataset that statistically matches the macroscopic properties of the observed stellar jitter, such as the overall root-mean-square (RMS) scatter, typical periodogram amplitudes, and the amplitude ratios between different CCF indicators. By maintaining a wide range for these parameters, we ensure the simulations retain sufficient generality to account for unexplained variability in the observed RV data, such as instrumental noise or unknown planetary signals. In Sect. 7.5, we discuss the current shortcomings of our simulations in matching some of these statistical metrics.

We simulated between 15 and 25 active regions for ϵ Eri, while we used between 80 and 110 for TZ Ari. This number is significantly larger in the latter case due to its longer observational time baseline of over 1400 d. Spot radii are drawn from a uniform distribution in the range of 3–6° for ϵ Eri and 6–8° for TZ Ari. We also vary the Q parameter, using values between 0 and 9 for ϵ Eri. For TZ Ari, we fix Q = 0 as it is an M5 dwarf, where spot-induced variability dominates (Herrero et al. 2016; Baroch et al. 2020) and faculae have been found to be dark, rather than bright, in 3D magnetohydrodynamic simulations (Johnson et al. 2021; Bhatia et al. 2026).

Fig. 3 Time series (left) and GLS periodograms (right) of the RV and activity indicators derived from the CCF for TZ Ari. The stellar rotation period (1.96 d) is marked in purple; the planetary orbital period from the literature (722.05 d; Quirrenbach et al. 2022) is marked in light blue; and the 0.1%, 1%, and 10% FAP levels are shown with dashed orange lines. The frequency axis of the periodogram show two scales split at 0.035 d−1 to better visualise long-term variability. The RVs display low-frequency signals associated with the Keplerian motion of TZ Ari b, while the CON and FWHM show lower-frequency stellar activity signals.

2.4.2 Noise injection

In addition to modelling stellar activity using StarSim, we simulate the impact of photon noise on the CCF. To ensure our synthetic observations are realistic, we inject noise patterns that match the statistical characteristics of the actual observations. We quantify the noise budget for each target star by computing the ratio of the mean observational error of each spectroscopic data to the total root-mean-square (RMS) of its observed time series (Sect. 3.4). This allows us to estimate the fractional contribution of noise to each time series and injects Gaussian noise according to this ratio. A summary of the estimated noise levels for each time series is provided in Appendix B.

3 Theoretical framework

When not considering the influence from instrumental effects, we assume that a spectroscopic observation O(v, t) in the exo-planet domain (e.g. spectrum, single line, CCF), depending on the vector of RV v, and the time of observation t, can be influenced by the following two effects: (1) an orbiting planetary companion that induces a non-chromatic Doppler shift in v, and which we denote in the following with ϵ, and (2) by stellar activity, which induces distortions on the spectral line shapes, including apparent line shifts, and which we denote in the following with κ.

We can decompose any spectroscopic observation, O(v, t), as $$O(v,t)={\displaystyle \sum }_{n=0}^{N-1}{a}_n(t)B_n(v)+E(v,t),$$(1)

where B_n(v) is a set of functions, hereafter basis functions, parametrized by time-dependent coefficients a_n(t). The term E(v, t) represents the residual uncertainty of the decomposition, and N is the number of decomposition terms, where N ≤ N_v. N_v refers to the number of elements at which u is evaluated. Strictly speaking, no information is lost if N = N_v. However, our aim is to condense the original observation into a smaller number of components with minimal loss of physical information. In the following, we explain how to isolate an orthonormal component from B_n(v), which contains only information on the velocity shifts, from other components only containing information on the shape distortions.

3.1 Isolating the line shift component

If we separate the contribution from Doppler shifts (O_ϵ) from the effects of stellar activity (O_∗)², we can write $$O(v,t)=c_0(t)I(v)+O_{ϵ}(v,t)+O_{\ast }(v,t)+E(v,t),$$(2)

where I(v) is a function of ones and c₀(t) is a constant. In the case where the observation is only affected by ϵ, this induces a shift: $$v_{ϵ}(t)=v+ϵ(t)I.$$(3)

Therefore, we can apply the Taylor expansion of O_ϵ with respect to ϵ around v: $$O_{ϵ}(v_{ϵ}(t))=O(v+ϵ(t)I)=O_r(v)+ϵ(t)O_r^{\prime }(v)+O({ϵ}^2),$$(4)

where O_r(v) is the undistorted rest frame, and ${O^{\prime }}_r(v)=\frac{d{O}_r(v)}{dv}$. In practical applications, this theoretical undistorted rest frame is approximated by our empirically derived time-averaged template CCF, T(v), as discussed in Sect. 4.1. Second-order expansion terms 𝒪(ϵ²) are negligible for the case of small ϵ, typical of planetary companions.

On the other hand, the stellar activity changes in v can be written in the most general way as $$v_{\ast }(t)=v+\kappa (t)I+s(t).$$(5)

Here, s(t) is a vector with an unknown number of components describing pure line shape distortions. Since we do not know the functional dependence of these components on physical parameters, a multivariate Taylor expansion is not feasible. The only constraints we imposed on s(t) is that they are orthogonal to the line shift component, ensuring that κ(t) captures the net line shift. In the most general way, we can then decompose the stellar activity part O_∗(v, t) of the observation with a Taylor expansion with respect to κ and capture the pure line shape distortions with G_n(v), hereafter the distortion basis function: $$O_{\ast }(v_{\ast },t)=O_r(v)+\kappa (t)O_r^{\prime }(v)+{\displaystyle \sum }_{n=2}^{N-1}{g}_n(t)G_n(v),$$(6)

where g_n(t) are the time-dependent coefficients describing G_n(v). If we combine Eqs. (2)–(6), we can write $$\begin{array}{c}O(v_{ϵ}+v_{\ast },t)=c_0(t)I(v)+c_1(t)O_r(v)+l(t)O_r^{\prime }(v)\\ +{\displaystyle \sum }_{n=2}^{N-1}{g}_n(t)G_n(v)+E(v,t),\end{array}$$(7)

where l(t) includes both the Doppler shift and the shift induced by the distorted line, i.e. l(t) = ϵ(t) + κ(t), and c₁(t) is a scaling factor that multiplies the undistorted rest frame O_r(v), which in this case is equal to 1³.

We thereby decomposed our observation O(v, t) into a basis function B that contains a line-shifted component captured by ${O^{\prime }}_r(v)$ and a distortion basis function, G(u). Comparing the above with Eq. (1), we found $$B(v)=[I(v),O_r(v),O_r^{\prime }(v),G_2(v),\dots ,G_{N-1}(v)],$$(8) $$a(t)=[c_0(t),c_1(t),l(t),g_2(t),\dots ,g_{N-1}(t)].$$(9)

We note, that ϵ and κ cannot be determined independently and are mathematically degenerate. This degeneracy can only be broken by relating the physical or statistical properties of the c₀, c₁ and g_n coefficients to κ (see Sect. 5).

3.2 Distortion basis functions

For the general framework many different choices for the distortion basis functions G_n(v) can be used. PCA isolates the principal modes of variability of O(v, t). But this method does not distinguish between the different sources of variability (Doppler, activity and instrumental) nor does it ensure the orthogonality of its components to ${O^{\prime }}_r(v)$. PCA was applied in Doppler-free observations on an orthogonalised O(v, t) with respect to ${O^{\prime }}_r(v)$ in Klein et al. (2024), and on the autocorrelation function of O(v, t), in Collier Cameron et al. (2021). An alternative option would be to use the Fourier basis, as proposed in Zhao & Tinney (2020). The basis consists of sine and cosine functions of different frequencies that are orthogonal to each other. However, since the basis is periodic overall v, it puts equal weight on the central part of the CCF and in the wings, which are heavily affected by noise.

For a Gaussian-like observation O(v, t), such as a CCF or a specific spectral line, we propose using the Hermite basis as the distortion basis function, which is defined as $$G_n(x)=\frac{e^{-x^2/2}}{\sqrt{2^{n}n!}\sqrt{\pi }}{H}_n(x),H_n(x)={(}^{-}{e}^{x^2}\frac{d^n}{d{x}^n}{e}^{-x^2},$$(10)

where x = (v − v₀)/σ is the normalised velocity coordinate centred at v₀, scaled by the width σ, and H_n(x) are the Hermite polynomials. These are defined by the recurrence relation $$H_{n+1}(x)=2x{H}_n(x)-2n{H}_{n-1}(x),$$(11)

where H₀(x) = 1 and H₁(x) = 2x.

A Hermite basis can be understood as a weighted generalisation of the monomial basis xⁿ, which is used to compute the central moments of a function. However, using pure monomials as basis functions leads to their divergence at large |x|. By introducing a Gaussian weight function to the monomial basis function, the Hermite basis effectively suppresses these divergences and delivers an estimate of the orthonormal central moments of a Gaussian-like function O(v, t). Crucially, a Gaussian multiplied by Hermite polynomials forms a basis that is mathematically orthogonal and complete. This completeness is a primary motivation for our choice, as it guarantees that any smooth deviation or asymmetry in the line profile can be accurately represented.

In principle, the distortion basis functions can be constructed using the n-th order derivatives of the Gaussian fit to the observation O(v), with n > 2, given that the G₀(x) and G₁(x) components are already accounted for in our general basis function B(v) by the O_r(v) and ${O^{\prime }}_r(v)$ components (Eq. (8)). For non-Gaussian profiles, such as M dwarfs similar to TZ Ari (Fig. 1), we apply a multi-Hermite basis. We model the reference O_r(v) by fitting a sum of two Gaussians, prioritising the simplest solution with the fewest free parameters. In the case of TZ Ari, the optimal fit combines a narrow positive Gaussian for the core with a broad negative Gaussian for the wings (Fig. 1).

3.3 Calculation of coefficients

We determined the optimal coefficients a_n through a χ² minimisation process: $${\chi }^2=\frac{1}{2}{\langle O(t)-a(t)B\rangle }^2.$$(12)

This process can be rewritten as $${\chi }^2=\frac{1}{2}{\langle O(t)\rangle }^2+\frac{1}{2}{\displaystyle \sum }_{n=0}^N{\displaystyle \sum }_{n^{\prime }=0}^N\left(a_n(t)a_{n^{\prime }}(t)\langle B_n,B_{n^{\prime }}\rangle -2a_n(t)\langle O(t),B_n\rangle \right).$$(13)

If we orthonormalise B_n, so $\langle B_n,B_{n^{\prime }}\rangle ={\delta }_{n{n}^{\prime }}$, this simplifies to $${\chi }^2=\frac{1}{2}{\langle O(t)\rangle }^2+\frac{1}{2}{\displaystyle \sum }_{n=0}^N\left(\langle a_n^2\rangle -2a_n\langle O(t),B_n\rangle \right).$$(14)

The coefficients that minimise χ² are then obtained by solving $$\frac{d{\chi }^2}{d{a}_n}=a_n(t)-\langle O(t),B_n\rangle =0$$(15)

and hence can be directly determined from the basis function B_n and the observation O(t). The different a_n(t) coefficients are independent from each other if B(v) is orthogonal. As explained in Sect. 3.2, the Hermite basis is an orthogonal basis, but it may not be orthogonal to the other components of the basis B(v), namely I(v), O_r(v), or ${O^{\prime }}_r(v)$.

The orthogonality of B(v) can be achieved using the algebraic Gram-Schmidt process, which transforms any set of linearly independent functions into an orthogonal set. Given a non-orthogonal function F_N, we orthogonalised it with respect to B_n using the iterative procedure $$f_N=F_N-{\displaystyle \sum }_{n=0}^N\langle B_n,F_N\rangle B_N,$$(16) $$B_N=\frac{f_N}{\Vert f_N\Vert },$$(17)

where f_N is the intermediate unnormalised function and ∥ · ∥ denotes its norm.

3.4 Coefficient errors

The uncertainties in the coefficients are obtained by propagating the errors from two independent sources: (1) the flux uncertainties of the observations, and (2) the uncertainties in the single-or multi-Gaussian fit to the template, from which the distortion basis functions are constructed. The total uncertainty in each coefficient can therefore be expressed as $$\delta a_n(t)=\sqrt{{\displaystyle \sum }_{i=1}^{N_v}[B_n{(v_{i)}\delta O(v_{i,}t)]}^2+[O(v_{i,}t)\delta B_n{(v_{i)}]}^2}.$$(18)

The uncertainty in the observed profile O(v, t) includes both photon noise and read-out noise (Appendix C), as described in Lafarga et al. (2020). To compute realistic noise estimates per pixel, O(v, t) must be sampled at the physical pixel size of the detector. However, the Gram-Schmidt orthonormalisation used to construct the basis functions B_n makes an analytical expression for δB_n intractable. Instead, we propagated the uncertainties from the Gaussian-fit parameters into the orthonormalised basis by computing the formal covariance from the fit and propagating it through the coefficient definition $$\delta a_{n,\text{basis}}^2(t)=J_n^T{\Sigma }_p{J}_n,$$(19)

where J_n is the Jacobian vector of partial derivatives of the coefficient a_n(t) with respect to the set of fitted parameters p, and Σ_p is the full covariance matrix of the fit. The Jacobian is defined as $$J_n=\frac{\partial a_n(t)}{\partial p}=\left[\frac{\partial a_n(t)}{\partial p_1},\frac{\partial a_n(t)}{\partial p_2},\dots ,\frac{\partial a_n(t)}{\partial p_k}\right],$$(20)

where k is the number of fit parameters. Since these derivatives are not analytically tractable, we estimated them numerically using finite differences. For each parameter p_k, we computed $$\frac{\partial a_n(t)}{\partial p_k}\approx \frac{a_n{(t)}_{|p_k+\Delta }-a_n{(t)}_{|p_k}}{\Delta },$$(21)

where ∆ is a small perturbation to the parameter p_k.

3.5 Explained variance ratio

To evaluate the contribution of each coefficient a_n to the reconstruction of the original time series, we compute the explained variance ratio (EVR), a metric analogous to that used in PCA. The EVR quantifies the fraction of the total variability in the full time series of O(v, t) captured by a given time series coefficient: $$Ev{R}_n=\frac{{{\displaystyle \sum }}_{j=0}^{N_T}{a}_n{(t)}^2}{1/N_v{{\displaystyle \sum }}_{i=1}^{N_v}{{\displaystyle \sum }}_{j=0}^{N_T}{\left[w(v_{i)}\left(O(v_i,t_j)-{\langle O(v_i)\rangle }_t\right)\right]}^2},$$(22)

where, in the most general case, a weight function w can be used for the observations. Further, ⟨⟩_t is the expected value over t, and N_T is the number of observations.

4 Application of the method

Several considerations must be taken into account when applying the theoretical framework to real observational data. The observations shall be flux normalised before applying this framework using, for example, AFS (Xu et al. 2019), rassine (Cretignier et al. 2020), or raccoon (Lafarga et al. 2020). If the continuum is not properly corrected to match the baseline flux level of O_r(v), this propagates to the computation of all a_n(t), given that this continuum level is a multiplicative factor in Eq. (1). Moreover, real observational data are affected by noise, which can bias the estimate of the EVR. We adopt the Gaussian fit to O_r(v) as a weight function, so we can properly weigh down the wings of the CCF, which are predominantly dominated by noise.

4.1 Template construction

As seen in Eq. (7), the decomposition is performed with respect to O_r(v), which represents the observation in an undistorted reference frame. This reference is also used to compute its first derivative, ${O^{\prime }}_r(v)$, in order to measure the line shift l(t). In general, however, we do not have access to an undistorted reference observation O_r(v). We therefore decomposed our observations with respect to a template T(v) that is constructed by averaging all observations of the time series, t, as explained in Sect. 2.2. We aim to reduce the influence of time-variable distortions and obtain a more stable representation of the stellar signal. The template provides a valid approximation to the undistorted reference profile required for the decomposition. Given that we use a template, generally, c₁ ≠ 0, and it will estimate the changes in amplitude.

Consequently, we approximated ${O^{\prime }}_r(v)$ by the first derivative of the template T(v). We note that T′(v) was calculated from the first derivative of a cubic spline on T(v). The first element of the Gram–Schmidt process is T′(v), so we could ensure that all other elements of our basis function are orthogonal to line shifts.

4.2 ϵ Eri

As can be seen in the left panel of Fig. 1, the template of ϵ Eri is in good agreement with a single Gaussian profile. We therefore adopted a single Hermite basis for its decomposition. In Fig. 4 the first six components of the decomposition of the ϵ Eri CCFs are shown in red, including the Hermite basis from G₂ to G₄. In Fig. 5, we show in red the individual and cumulative EVR per component. The individual EVR peaks at the c₁ component, connected to the template T, and decreases monotonously as we go to higher order distortion coefficients. The first five coefficients explain 91.7% of the variance, while 11 components account for 97.8%. Notably, our method enables the quantification of the contribution to the variance by line shifts via the l coefficient, which alone accounts for 19.3% of the total variance.

In Fig. 6, we show (i) the reconstruction of the differential (template-subtracted) CCFs at six random epochs by adding the contribution from the different basis function components (left column), (ii) the time series of the first eight coefficients (middle column), and (iii) their corresponding GLS periodograms, (right column). Coefficients up to g₆ trace the stellar rotation period, marked with purple lines, with l, g₂ and g₃ also showing sensitivity to the first harmonic at ∼5.6 d. The slight shifts in the exact periodogram peaks across different coefficient orders can be physically attributed to stellar differential rotation (as reported in Croll et al. 2006) and active region evolution. We also note that the baseline offset term, c₀, also exhibits sensitivity to stellar rotation. While mathematically defined as a constant, c₀ naturally absorbs changes in the overall integrated area of the CCF, whether caused by imperfect flux normalisation or, as in this case, stellar activity driving variations in both the continuum level and line amplitude.

As expected from their lower EVRs, higher-order coefficients contribute less to the signal reconstruction and are increasingly affected by noise, therefore reducing their sensitivity to rotationally modulated stellar activity. For a summary of the general statistics of the distortion coefficients and the identified signals from a pre-whitening process we refer to Appendix B.

Figure D.1 shows the correlation between decomposition coefficients and classical CCF indicators. The strongest correlations are CON with c₁, RV with l, FWHM with g₂, and BIS with g₃, confirming that the decomposition recovers the classical CCF activity indicators. The strongest correlation between decomposition coefficients is between l and g₃, and between c₁ and g₂. Beyond this, higher-order coefficients extend the sensitivity to distortions not captured by classical CCF activity indicators.

Fig. 4 First six basis components for ϵ Eri (Hermite, red) and TZ Ari (multi-Hermite, blue).

Fig. 5 Explained variance ratio (solid) and cumulative variance (dashed) per component for ϵ Eri (Hermite, red) and TZ Ari (multi-Hermite, blue).

4.3 TZ Ari

Figure 1 shows the template used to compute the basis functions. It displays the characteristic shape of cool M dwarfs (Lafarga et al. 2020). As a result, a double Hermite basis is used to fit this template, consisting of a narrow positive-amplitude Gaussian for the core and a broader negative-amplitude Gaussian for the wings. In Fig. 4, the first six components of the decomposition of the TZ Ari CCFs are depicted in blue, showing a distinct shape for the multi-Hermite case. Figure 5 illustrates the EVRs for the multi-Hermite decomposition coefficients applied to TZ Ari. The EVR peaks at the c₁ coefficient, explaining 47.1% of the EVR. The l coefficient alone accounts for 19.6% of the total variance. The g₂ coefficient contributes only marginally, while g₅ provides a comparatively large contribution to the variance. The cumulative EVR reaches 86.8% at this coefficient, reaching the 97.1% mark if we extend it up to the g₉ coefficient.

Figure 7 shows the differential CCF reconstructions (left), coefficient time series (middle), and their corresponding GLS periodograms (right). The stellar rotation period (1.96 d) and its daily alias are clearly recovered in the low-order coefficients; c₁, l, g₂ and g₃. Long-period signals are also observed in c₁, l, g₂, and g₄ (Appendix B). The exoplanetary signal at 772.05 d is significantly detected in the GLS periodogram of l only, as expected.

Figure D.2 displays correlations between decomposition coefficients and classical CCF activity indicators. As for the case of ϵ Eri, the strongest correlations are found between CON and c₁, RV and l, FWHM and g₂, and BIS and g₃. However, the correlations between FWHM and g₂, and between BIS and g₃ are not as strong as in the ϵ Eri case. The strongest correlation between decomposition coefficients is between c₀ and c₁, and between l and g₃.

Fig. 6 Cross-correlation function decomposition of ϵ Eri with a Hermite basis. Left: contribution of consecutive basis components to the differential CCF reconstruction. We displayed six illustrative epochs randomly selected to visualise the temporal variability of the line-shape distortions. The observed differential CCFs are shown in grey, while the reconstructed components are colour-coded according to the observation time. Middle: time series of the corresponding coefficients. Right: GLS periodograms of the coefficients. The 0.1%, 1%, and 10% FAP levels are indicated in orange, and the stellar rotation period and its first harmonic are marked in purple.

5 Neural networks for stellar activity correction

In this section we introduce the NN framework. The goal is to predict and remove the activity-induced RV shifts using the information from the line-shape distortion coefficients.

5.1 Architectures

We have developed a deep-learning model that we refer to as Convolutional-Attention Network for STellar Activity Removal (CANSTAR). An overview of the CANSTAR architecture is provided in Fig. 8. Several (N_c) convolutional layers apply N_f localised filters to the input time series, which consist of the amplitude coefficient c₁(t) and the distortion coefficients g_i(t). We explicitly exclude the baseline offset term, c₀(t), from the inputs because it is highly susceptible to flux normalisation issues and sensitive to noise sources not accounted for in our simulations. These filters extract short-term, high-frequency features in the input according to the filter size N_s. The convolutions are applied with appropriate padding and stride, ensuring that the output feature maps keep the same temporal dimension as the input. Each convolutional layer uses multiple filters, producing a set of feature maps that capture different local attributes of the data.

Following the convolutional layers, the extracted feature maps are passed into N_te transformer encoder layers. These layers consist of a multi-headed self-attention mechanism followed by a position-wise feed-forward network, as defined in Vaswani et al. (2017). The self-attention mechanism captures dependencies across the full temporal domain, allowing the network to model both short- and long-term correlations in the data. The embedding dimension of the attention layer is matched to the number of convolutional filters, and N_h attention heads are used to extract different temporal features. Finally, the output of the transformer block is flattened and passed through a linear projection layer, producing the predicted κ(t) time series. This hybrid architecture follows principles similar to those proposed by Liu et al. (2020) and Guo et al. (2022), who demonstrated the effectiveness of combining convolution and attention in the analysis of time series and image data.

We also tested a simpler fully connected network (FCN). This architecture processes only instantaneous snapshots of the line-shape distortion coefficients, g_i, without explicitly modelling their temporal evolution. The FCN architecture consists of a sequence N_ff of feed forward layers. All input coefficients from a single observation are concatenated into a one-dimensional feature vector, which is fed into the network. The outputs of each layer are fully connected to the inputs of the next by N_n number of neurons. The final layer of the FCN is a linear layer which predicts the final κ.

Fig. 7 Same as Fig. 6 but for the multi-Hermite basis applied to TZ Ari. We also mark in light blue the planetary orbital period from the literature (722.05 d; Quirrenbach et al. 2022).

Fig. 8 Schematic representation of CANSTAR (bottom), a convolutional attention network designed to predict the stellar activity contribution to line shifts, κ(t), from the time series of distortion coefficients, gi(t). The input coefficients are treated as independent channels and processed through a series of one-dimensional convolutional layers that extract local temporal features. The convolutional outputs are then passed to a transformer encoder module, which consists of a multi-head self-attention part (top) with different ‘heads’, where each one of them focuses on different temporal dependencies within the data, and their outputs are combined, normalised, and propagated through additional attention layers. Finally, the resulting features are flattened and projected through a linear layer to produce the predicted κ(t) time series.

5.2 Neural network training

Training is performed using synthetic datasets generated with the StarSim code (Sect. 2.4.1). For each target star, we generate a dataset of 500 000 time series simulations. The dataset is divided into training (80%), validation (10%), and test (10%) sets. The training set is used to find the optimal parameters of the NN based on the minimisation of the Mean Squared Error (MSE) between the true labels, i.e. the κ(t) component of l(t), and the prediction by the network. We apply z-score normalisation to the input and output dataset, meaning we subtract the mean value of each simulated dataset and rescale them according to the standard deviation of the full dataset.

We optimised all network architectures using the Adam optimiser (Kingma 2014) implemented in PyTorch (Paszke et al. 2019). We selected as the final model the model just before the improvement in validation MSE plateaus while the training MSE continues to decrease significantly, ensuring we avoid the over-fitting regime. We reserved this test set for the evaluation of the performance of this final network on unseen data. We saved the residual relative error (RRE) from the final model, defined as the ratio of the residual RMS after correction to the original RMS (Perger et al. 2023).

5.3 Hyper-parameter choices

Each architecture is defined by several hyper-parameters. We performed a systematic hyper-parameter search to determine these values. The tested and final selected configurations are summarised in the upper and lower parts of Table 2 for CANSTAR and the FCN, respectively. The final selection was based on minimising the validation set MSE, prioritising the simplest effective choices.

6 Results

We trained separate CANSTAR models for each target star (ϵ Eri and TZ Ari) for different noise-level regimes (noiseless and similar to the observations; see Sect. 2.4.2) and including different numbers of distortion coefficients. The aim with this training strategy was to evaluate the ability of CANSTAR to mitigate stellar activity as a function of the target star and S/N of the observations. In order to test the capability of modelling stellar activity by the different distortion coefficients, we iteratively trained CANSTAR networks using more coefficients as inputs, following the natural order of the coefficients from the Hermite and multi-Hermite formalism. We repeated this same training strategy for the FCN.

After the training phase, we selected the best-performing network trained on noise similar to the observations, defined as the model that minimises the RRE while requiring a smaller number of distortion coefficients. To quantify the stability of the solution, we trained five independent instances of this architecture. Although the variance in the predictions on the synthetic test set is negligible, it can have an effect when predicting on the real observational data. To account for this, we applied the full ensemble of five networks to the data. Furthermore, to propagate the observational uncertainties, we generated 1000 Monte Carlo realisations of the input distortion coefficients, perturbed according to their measurement errors. Each of the five networks generates predictions for all 1000 realisations. We adopted the median of this combined distribution as the final activity model and used the standard deviation as the total uncertainty.

6.1 Simulated data

We show the performance on the simulated test set of CANSTAR and FCN in Fig. 9. In the noiseless case, CANSTAR achieves an RRE of 3.7% for ϵ Eri and 3.3% for TZ Ari, when trained with 15 coefficients. Notably, the performance already saturates at 7 coefficients, reaching an RRE 5.9% for ϵ Eri and 3.9% for TZ Ari. By contrast, the FCN achieves an RRE of 12.7% for ϵ Eri and 4.9% for TZ Ari after training with 15 coefficients. CANSTAR outperforms the FCN in RRE reduction, requiring fewer coefficients to reach the same level of activity correction.

With noise matched to the observed RV uncertainties, CANSTAR reaches an RRE of 12.9% for ϵ Eri and 35.7% for TZ Ari, when trained with 15 coefficients. On the other hand, the FCN network reaches an RRE of 38.8% for ϵ Eri and 49.1% for TZ Ari. In both regimes, the RRE reduction plateaus at a smaller number of coefficients than in the noiseless case, highlighting the fact that higher-order coefficients are increasingly dominated by noise and contribute little to the correction.

There is also a difference in the performance of the network between the two target stars. For ϵ Eri, the most significant coefficients when correcting for activity are the c₁, g₂ and g₃ coefficient. In the case of TZ Ari, c₁ and g₂ are not as correlated to the line shift component as g₃, which is the most informative for training both CANSTAR and the FCN network. We also note that we reached a lower overall RRE reduction for the cases of ϵ Eri with observational noise when compared to TZ Ari. This is because the former is a much brighter star with significantly higher S/N observations.

6.2 ϵ Eri

For ϵ Eri, we selected the network trained on the first five distortion coefficients as the best-performing configuration. We applied z-score normalisation to the real observed datasets in order to minimise the discrepancy in the scale between observed and synthetic datasets (Appendix E). We explored the inclusion of a linear scaling parameter to fine-tune the amplitude of the predictions given the loss of the absolute scale due to this. The scaling parameter converged to a value of α = 0.9. This value is close to unity, as expected, given that the RV variability of ϵ Eri is dominated by stellar activity on these timescales. Therefore, the z-score normalisation used during training accurately reflects the scale of the observed activity signal (Appendix E).

CANSTAR’s correction for ϵ Eri is shown in Fig. 10. The predicted RV time series closely matches the observed data, and the GLS periodogram reveals that the model effectively captures the dominant stellar rotation signal and its first harmonic. After subtracting the predicted κ component from the observations, we obtain an RRE of 52.5%, resulting in a final RMS of 2.46 m s⁻¹ (Table 3). The power of signals associated with stellar activity, which peak at 12.53 d and 6.26 d, is significantly diminished, decreasing below the 10% FAP level. No signal remains in the periodogram above the 10% FAP level. The residuals are in agreement with the ϵ Eri b solution presented in Thompson et al. (2025), which just produces a trend for the narrow observing window.

We injected sinusoidal signals directly into the RV time series to quantify the detection limits at different frequency ranges. We injected sinusoidal signals with frequencies ranging from 0.5 to 1/T_base (where T_base is the time baseline of the observations), in steps of 1/5T_base. The semi-amplitude grid spans from 0.5 m s⁻¹ to the RMS of the observations, 4.5 m s⁻¹, in steps of 0.1 m s⁻¹. For each combination of frequency and semi-amplitude, we injected ten evenly spaced phases from 0 to 2π. The phase-averaged results are illustrated in Fig. 11 as heat maps showing the relative difference between the injected and recovered semi-amplitudes (K) and frequencies for the CANSTAR activity-corrected residuals. We defined the detection limit as the minimum injected amplitude required to recover the semi-amplitude with a specific precision (e.g. 10% or 30%), provided the frequency is also retrieved within a 10% relative error. Using this criterion, the detection limits are 2.51 ± 0.43 m s⁻¹ for a 10% precision threshold and 1.99 ± 0.41 m s⁻¹ for a 30% threshold. As expected, the recovery of planetary signals becomes more challenging around the stellar rotation period and its harmonics due to the inherent degeneracy between Keplerian and activity-induced variability. We also show in Appendix F the detection limits for the null hypothesis case. The detection limits are 3.99 ± 0.18 m s⁻¹ (10% precision) and 3.87 ± 0.18 m s⁻¹ (30% precision).

Fig. 9 Residual relative error of ϵ Eri (top) and TZ Ari (bottom) StarSim data after correction with CANSTAR (solid line with circle marker) and an FCN (dashed line with square marker) for the noiseless case (blue) and for the case with noise equivalent to the observations (green).

Fig. 10 Top: Radial velocity time series (left) and GLS periodogram (right) of the HARPS ϵ Eri data (blue) compared to the stellar activity signal predicted by CANSTAR (green). Bottom: Residuals after subtracting the prediction (grey), with their corresponding periodogram (right). The 10%, 1%, and 0.1% FAP levels are indicated by dashed orange lines and the stellar rotation period (11.2 d), and its first harmonic (5.6 d) is marked in purple. The residuals are compatible with the ϵ Eri b solution (red; Thompson et al. 2025).

Fig. 11 Detection limits in the ϵ Eri residuals after applying CANSTAR’s correction, shown as a function of the difference between injected and retrieved frequencies (left) and semi-amplitudes (right) for the different injected frequencies and semi-amplitudes sinusoidal signals.

6.3 TZ Ari

For TZ Ari, we similarly select the network trained on the first three order coefficients as the best-performing configuration. Unlike ϵ Eri, the RV variability of TZ Ari is dominated by the high-amplitude Keplerian signal of the planet (K = 21.11 m s⁻¹; Quirrenbach et al. 2022). Consequently, we need to perform a joint fit for the Keplerian solution of the planet together with the CANSTAR activity correction. We use the MCMC sampler emcee (Foreman-Mackey et al. 2013), simultaneously solving for the orbital parameters of TZ Ari b and the optimal value of the scaling parameter of the CANSTAR activity correction, α. For the MCMC parameter optimisation, we employed 700 walkers, each run for 32 000 steps to ensure convergence of the chains. The first 2000 steps were discarded as burn-in. The adopted priors and the resulting optimised parameter values are listed in Table 4, together with the Kep-lerian solution we retrieve if we do not apply the CANSTAR activity correction. We also compare these results to a GP framework using a simple harmonic oscillator (SHO) kernel, reproducing the strategy of Quirrenbach et al. (2022) but restricting their multi-instrument dataset to CARMENES only (Appendix G).

We compare in Fig. 12 the posterior distribution of the common parameters between the CANSTAR and GP model. CANSTAR activity correction allows for a better determination of the period and semi-amplitude of the planet, compared to the GP modelling. The GP solution outperforms CANSTAR in reducing the RV scatter, as illustrated by its lower fitted jitter term, which captures additional noise sources not captured by the model. This is expected given the strong flexibility of GPs at filtering white noise and being prone to overfitting (Blunt et al. 2023). The posterior distributions of the CANSTAR + Keplerian model parameters are shown as corner plots in Appendix H.

Figure 13 presents the time series and GLS periodograms of the CARMENES RV observations (blue), together with the CANSTAR activity correction (green), the Keplerian solution (red), and the residuals (grey). In the time series plot (left), the green curve represents the sum of the CANSTAR correction and the Keplerian orbit. This allowed us to visualise the total fit to the observations. The resulting model successfully captures both the Doppler and stellar activity signals present in the data. Notably, the planetary origin of the long-period signal is confirmed by the periodogram (right panel), where the CANSTAR model (green) is shown in isolation and is clearly insensitive to the planetary signal. After applying the correction, the RV residuals show an RRE of 62.4%, measuring the ratio with respect to the variability after subtracting the fitted Keplerian model (Table 3). The residual GLS periodogram of the CANSTAR correction is largely flat, showing only a minor insignificant peak at 64.7 d (FAP = 11.7%). In contrast, the residuals of the GP + Keplerian model (Fig. 14) display two distinct signals above the 1% FAP level at 37.2 d and 41.4 d. These periodicities correspond to yearly aliases of one another (1/37.2 ≈ 1/41.4 + 1/365) and are likely caused by uncorrected stellar activity, given the significant variability observed in the FWHM and g₂ indicators at 40.18 d (Figs. 3 and 7).

During the final CARMENES observing season, when the star was monitored more intensively, the temporal sampling becomes sufficient to mitigate rotationally modulated stellar activity more effectively with CANSTAR. A zoom-in of this last season is shown in Fig. 13, where the RRE is further reduced to 33.4%.

Fig. 12 Histograms showing the posterior distribution of the shared parameters between the CANSTAR and GP correction methods used to model the TZ Ari RV data.

7 Discussion

7.1 The value of orthogonal decomposition

Standard activity indicators derived from the CCF rely on a simple Gaussian fit. This approach is suboptimal for capturing complex line deformations, requiring ad hoc diagnostics such as the BIS to quantify asymmetries. Furthermore, as shown in Fig. 1, the Gaussian model fails to accurately reproduce the non-Gaussian, double-humped CCFs characteristic of many M dwarfs.

A key advantage of the theoretical framework presented here is its generality, allowing for the construction of basis functions tailored to the specific spectral type of the star. We demonstrate that the standard Hermite basis is optimal for describing Gaussian-like CCFs, while M dwarf CCFs exhibiting pronounced double humps are better modelled by the multi-Hermite basis. For the specific case of TZ Ari, the humps are largely symmetric, making a two-Gaussian fit sufficient to construct G₀. However, more complex or asymmetric CCF morphologies might necessitate three or more Gaussians. A robust and automated approach to determine the optimal number of Gaussians for the mother function is to perform a systematic model selection – for instance, by minimising the Bayesian Information Criterion (BIC) – when fitting the time-averaged template CCF.

This theoretical framework enables the recovery of the classical CCF activity indicators, as seen in the high level of correlation between c₁-CON, g₂-FWHM and g₃-BIS (Sects. 4.2 and 4.3), in addition to extending the analysis to higher order variability. Moreover, the EVR allowed us to quantify the actual number of distortion basis components needed to describe the observed CCFs.

While a mathematical degeneracy remains between true Doppler shifts induced by planets (ϵ) and apparent shifts caused by stellar activity (κ), the orthogonal decomposition ensures that our inputs to the NN, the coefficients g_n(t), capture shape information that is mathematically independent of the target variable l(t). Using the first three distortions coefficients, (c₁, g₂ and g₃), we are effectively using similar information to the classical CCF activity indicators, but extending the analysis to higher order coefficients allowed us to improve the activity mitigation performance.

More generally, the distortion coefficients are not limited to the specific architecture proposed. They could serve as inputs for other stellar mitigation techniques, such as the multi-dimensional GP framework. Finally, this parametrisation is computationally efficient and agnostic to the specific pipeline used. It can be easily implemented in standard reduction codes such as raccoon or the ESO DRS.

7.2 Temporal awareness in activity correction

The comparison between CANSTAR and the FCN highlights the critical role of temporal information. As shown in Fig. 9, the FCN, which treats observations as independent snapshots, fails to reach the same correction performance as CANSTAR, particularly in noisy regimes.

Stellar activity is inherently a time-correlated process. Active regions evolve, migrate, and reappear over the stellar rotation period. The inclusion of self-attention mechanisms allows CANSTAR to ‘learn’ these temporal correlations, effectively using the past and future context of the distortion coefficients to mitigate stellar activity occurring at different timescales. This capability is illustrated by the fact that CANSTAR effectively saturates its performance with fewer input coefficients than the FCN. By including the temporal context, the network can extract more information from the lower-order coefficients (c₁, g₂, g₃), reducing the reliance on higher-order terms that are more affected by photon noise.

For exoplanet host stars with sparse RV monitoring, temporal correlations may be lost. In such cases, activity mitigation strategies that do not rely on temporal continuity remain necessary. Our results with the FCN demonstrate the advantage of using the full set of distortion coefficients as inputs compared to classical CCF activity indicators, which are equivalent to using only the lower-order terms (c₁, g₂, g₃).

7.3 Noise estimation in neural network predictions

We emphasise the critical importance of correctly propagating photon noise from the spectral level to the CCF and subsequently to the derived activity indicators (see Sect. 3.4 and Appendix C). These precise noise estimates allowed us to inject realistic noise into the synthetic dataset (Sect. 2.4.2). Training on data with realistic noise properties is essential to prevent the network from over-fitting to specific noise patterns. Furthermore, these estimates enable us to quantify the uncertainty of the NN predictions. As detailed in Sect. 6, we achieve this by evaluating the network on noisy realisations of the input data (propagating measurement error) and by training an ensemble of independent networks to estimate the model variance.

Fig. 13 Radial velocity time series and GLS periodograms for TZ Ari showing (a) the complete multi-season dataset and (b) a detailed view of the final observing season. For each panel, the top row displays the CARMENES data (blue), the CANSTAR prediction (green), and the Keplerian solution (red). In the time series, the green curve includes the Keplerian signal to illustrate the full fit, whereas in the periodogram it represents the activity model alone. The bottom row shows the residuals after subtracting both the CANSTAR activity correction and the Keplerian model (grey). The 10%, 1%, and 0.1% FAP levels are indicated by dashed orange lines. The stellar rotation period (1.96 d) is marked with a purple dashed line, and the newly derived planetary orbital period (789.92 d) is marked in light blue (not shown in panel (b) due to its shorter time baseline).

7.4 Comparison with previous studies

This study builds upon the previous work developed in Perger et al. (2023). We have explored different input data, transitioning from classical CCF activity indicators (CON, FWHM and BIS) to time series of distortion coefficients, and different NN architectures, from a CNN to a convolutional attention network. Both studies relied on StarSim data for training. However, in this work, we utilised a newer version of the code that has been improved, both in its physical model as well as on the practical implementation (Gomes et al., in prep.).

We can directly compare the performance of Perger et al. (2023) against CANSTAR regarding stellar activity correction on the ϵ Eri dataset (Fig. 15). The CANSTAR correction matches the observed data better, resulting in an RRE reduction to 52.5% of the original variability, compared to the 67% reduction level of the Perger et al. (2023) results⁴. The residuals from both approaches show significant differences. The GLS periodogram of the CANSTAR residuals shows a clear mitigation of the rotational period (FAP > 10%). On the other hand, the peak at the rotation period remains highly significant (FAP < 0.1%) in the GLS periodogram of the Perger et al. (2023) residuals. The residuals from both methods display a low degree of correlation, with a Pearson correlation coefficient of 0.38.

In Sect. 6.3, we compared the CANSTAR activity correction on TZ Ari against the GP framework with the SHO kernel from Quirrenbach et al. (2022), but we limited the analysis to CARMENES. CANSTAR achieves compatible but more precise determinations of the semi-amplitude and period of the Keplerian parameters (Fig. 12). The GP achieves a larger reduction of RRE, as illustrated in the lower fitted jitter term. This is expected due to the flexibility of GP modelling solely on RV data, which can lead to over-fitting leaving small residual amplitudes (Blunt et al. 2023). However, two residual signals in the GLS periodogram remain above the 1% FAP level after the GP solution, with periods 37.2 d and 41.4 d. These periods are likely yearly aliases of a stellar activity signal, as evidenced by significant variability in the g₂ and FWHM indicators at 40.18 d (Table B.2). In contrast, CANSTAR effectively removes these signals because its correction is explicitly conditioned on the line-shape distortion coefficients. This demonstrates that CANSTAR provides a physically motivated correction that avoids the over-fitting pitfalls of pure RV modelling.

Fig. 14 Residual time series (left) and GLS periodogram (right) after subtracting the GP + Keplerian model to the TZ Ari RV data (grey). The 10%, 1%, and 0.1% FAP levels are indicated by dashed orange lines. The stellar rotation period (1.96 d) is shown in purple, the newly derived planetary orbital period (789.92 d) is in light blue, and the two signals exceeding the 1% FAP level (37.2 d and 41.4 d) are marked in black. These likely correspond to stellar activity given the significant (40.18 d) periodicity observed in the FWHM and g2 indicators.

Fig. 15 Top: radial velocity time series (left) and GLS periodogram (right) of the HARPS ϵ Eri data (blue) compared to the stellar activity signal predicted by CANSTAR (green) and Perger et al. (2023, magenta). We note that the Perger et al. (2023) did not provide uncertainties to their NN predictions. Bottom: Residuals after subtracting the prediction with their corresponding periodogram (right). The 10%, 1%, and 0.1% FAP levels are indicated by dashed orange lines.

7.5 Bridging the synthetic gap

We find a difference in performance when evaluating the CANSTAR network on synthetic versus real data. In our noiseless simulations, CANSTAR achieves nearly perfect correction, and even in datasets with injected noise matching the observations, the network achieves an RRE of 12.9% for ϵ Eri and 35.7% for TZ Ari. However, when applied to the real datasets, the RRE becomes 52.5% and 61.6%, respectively. This performance gap indicates that while StarSim provides a physically realistic output, it still lacks certain stellar activity and instrumental effects present in reality.

First, our noise injection procedure currently accounts only for photon noise measured at the spectral level (white noise). However, real observations are affected by instrumental systematics and atmospheric effects that introduce correlated noise into the data. Since our network is trained solely on white noise, it may be less effective at filtering out these complex, non-Gaussian systematic trends. Future synthetic datasets should therefore integrate instrumental systematics to better simulate the red noise floor of the specific spectrograph.

Second, current StarSim simulations rely on 1D PHOENIX stellar spectra (Husser et al. 2013) combined with modified bisector line shapes to simulate the photosphere, spots, and faculae. However, spots and faculae are direct manifestations of the interplay between magnetic fields and the complex 3D structure of the photosphere (Witzke et al. 2022). This dependence on 1D atmospheric models is a shared limitation among current simulators, such as SOAP-GPU (Zhao & Dumusque 2023) and SOAP 4.0 (Cristo et al. 2025).

A significant step towards bridging this physical gap is the use of synthetic spectra derived from 3D magnetohydrodynamic (MHD) simulations. The upcoming version of StarSim (Gomes et al., in prep., Stucki et al., in prep.) is designed to incorporate spectra from the MURaM MPS-ATLAS 3D MHD models (Witzke et al. 2024). By explicitly modelling the 3D structure of the stellar atmosphere driven by convection, this approach will improve the physical prescription of the photosphere, spots and faculae, while enabling the inclusion of other convective phenomena such as granulation and super-granulation. Simulating these phenomena will be essential for extending this framework to less active stars (Reinhold et al. 2019; Meunier & Lagrange 2020).

On the other hand, accurately modelling very active, flaring M dwarfs would require StarSim to simulate coronal heating driven by high-intensity magnetic fields. The self-consistent simulation of stellar chromospheres and flare events remains a significant limitation shared across the entire field (Allen et al. 2026).

8 Conclusions

We have introduced CANSTAR, a novel framework that exploits information on CCF line-shape distortions and their temporal evolution. The framework consists of two main steps: (1) separating pure Doppler shifts from line-shape distortions using an orthonormal basis expansion, which enables the subtraction of distortion coefficients time series, and (2) modelling the temporal evolution of these stellar activity distortions with a convolutional attention network trained on synthetic datasets generated with StarSim.

We have demonstrated that CANSTAR achieves near-perfect stellar activity correction on simulated data, outperforming an FCN that does not explicitly model temporal correlations. The performance worsens when noise resembling that of real observations is injected, primarily because higher-order coefficients become increasingly noise dominated.

When trained on synthetic noisy data and applied to real observations of ϵ Eri and TZ Ari, CANSTAR successfully mitigated a significant fraction of the stellar activity signal. For ϵ Eri, we improved the activity correction provided in Perger et al. (2023) by reducing the RRE from 67% to 52.5%. Applying the CANSTAR correction improves the detection limit for a 10% relative precision in K to 2.51 ± 0.43 m s⁻¹, compared to the baseline of 3.99 ± 0.18 m s⁻¹ obtained for the uncorrected (null hypothesis) case.

For TZ Ari, the network successfully disentangles planetary Doppler shifts from activity-induced signals, yielding improved orbital parameter estimates for TZ Ari b compared to a GP fit. This demonstrates that CANSTAR outperforms the current state-of-the-art solution, enabling both robust activity correction and precise planetary characterisation. After subtracting the Keplerian solution, CANSTAR achieves an RRE of 62.4% on the full dataset. Notably, this residual scatter is further reduced to 33.4% for the final observing season, where the higher cadence allowed for proper sampling of the 1.96 d stellar rotation period.

Beyond the integrated framework, the individual components of CANSTAR offer independent utility. The distortion coefficients obtained from the orthonormal basis expansion can reproduce classical activity indicators while providing sensitivity to additional stellar phenomena via higher-order terms. Furthermore, the convolutional attention network has a general time–aware architecture capable of capturing variability on multiple timescales. Its inputs are not limited to CCF coefficients, but it can be readily adapted to use photometry, chromospheric indices, or line-by-line activity indicators.

Looking forward, there is a clear path to bridge the performance gap between simulations and real data. Future StarSim developments will incorporate spectra from 3D MHD simulations, allowing us to model less active stars where granulation and faculae dominate. By integrating these physical effects along with instrumental systematics into the training sets, we expect to further enhance the network’s predictive power.

Ultimately, this work demonstrates that NNs are not just a future prospect but a present capability, and they already outperform established state-of-the-art solutions, such as GPs, in regimes dominated by complex stellar activity. By effectively using temporal context and high-order line shape distortions, CANSTAR offers a robust pathway to disentangling planetary signals from stellar activity effects. This approach, with further maturation and work closing the gap between synthetic and real data, can be instrumental in pushing detection limits towards the photon-noise floor, which would enable the discovery of lower-mass planets and eventually unlock the domain of exo-Earths.

Acknowledgements

The authors thank the anonymous referee and the editor for their constructive feedback and careful reading of the manuscript. We also warmly thank Pedro Figueira for providing valuable comments that helped improve this work. J.B.-P., M.P., G.A.-E., I.R., J.C.M., O.P. and S.S. acknowledge financial support from Spanish grants PID2021-125627OB-C31 funded by MCIU/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”, PID2024-158486OB-C31 funded by MCIU/AEI, by the programme Unidad de Excelencia María de Maeztu CEX2020-001058-M and by the MaX-CSIC Excellence Award MaX4-SOMMA-ICE, by the Generalitat de Catalunya/CERCA programme, and by the European Research Council (ERC) under the European Union’s Horizon Europe programme (ERC Advanced Grant SPOTLESS; no. 101140786). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. J.B.-P., M.P. and G.A.-E. also acknowledge financial support from Spanish grant PID2020-120375GB-I00, funded by MCIU/AEI, and Consolidación 2022 CNS2022-136050. J.B.-P. also acknowledges financial support from Spanish grant PRE2022-101942 funded by MICIU/AEI/10.13039/501100011033 and ESF+. M.L acknowledges support by the UKRI (Grant EP/X027562/1). The data production, processing and analysis tools for this paper have been developed, implemented and operated in collaboration with the Port d’Informació Científica (PIC) data center. PIC is maintained through a collaboration agreement between the Institut de Física d’Altes Energies (IFAE) and the Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT).

References

Appendix A TESS photometry

The Transiting Exoplanet Survey Satellite (TESS) is performing an all-sky survey in search for transiting exoplanets around the closest and brightest stars (Ricker et al. 2015). TESS first observed TZ Ari in Sectors 70 and 71 at the start of the mission’s year 6, spanning from September to November 2023. We analyse the Pre-search Data Conditioning Single Aperture Photometry (PDCSAP) flux provided by the mission, which has been corrected for trends specific to each CCD and observing sector. The light curves and GLS periodograms are dominated by a stable signal induced by the rotation period. We also see a large flare event at around BJD ≈ 2460210.

Based on the identification of the highest power peak in the periodograms of the two sectors and estimating the uncertainty from the FWHM of the dominant peak, we report a rotation period of 1.96±0.07 d, which is in agreement with the reported period from high-resolution spectroscopy and ground-based photometry analysis in Quirrenbach et al. (2022).

Appendix B Statistics of CCF products

We show the main statistics on the derived distortion coefficients and classical CCF activity indicators of the test targets of the study ϵ Eri and TZ Ari. These statistics consist of the root mean square (rms) of the time series, the mean error, the EVR, the cumulative EVR, and the statistics from the pre-whitening process, consisting of the RRE, the identified periods and the fitted amplitudes.

Appendix C CCF error propagation

The noise budget of the CCF is determined by the uncertainty of the measurement of electrons on the spectrograph detector. The CCF value at a given RV point is computed as the weighted sum of the flux from N_p pixels that overlap with the N_m spectral lines defined in the weighted mask. This can be expressed as $$\begin{array}{l}\text{CCF}(v)={\displaystyle \sum }_{l=1}^{N_m}{\displaystyle \sum }_{x=1}^{N_p}{m}_l{f}_x{\text{Δ}}_{xl},\\ \end{array}$$(C.1)

where m_l is the weight of the mask for spectral line l, f_x represents the flux of pixel x in the analog-to-digital Unit (ADU), and ∆_xl is the fraction of pixel x covered by the mask line l after shifting it by a v shift.

To estimate the uncertainty on the CCF profile, we consider the noise variance of each individual pixel. The flux in ADU, f_x, is related to the number of photo-electrons, $n_{{\text{e}}^{-}}$, via the detector gain g (in ADU/e⁻): $$f_x=g\cdot n_{e^{-}}.$$(C.2)

Given that the photo-electrons follow Poisson statistics, the uncertainty in the electron count is $\delta n_{e^{-}}=\sqrt{n_{e^{-}}}$. Converting this to ADU and including the detector readout noise, RON, the total variance for a single pixel x is given by: $$\delta f_x^2=g^2\delta n_{e^{-}}^2+g^2{\text{RON}}^2=g^2{n}_{e^{-}}+g^2{\text{RON}}^2=g{f}_x+g^2{\text{RON}}^2.$$(C.3)

Since the noise contributions from individual pixels are statistically independent, the total error on the CCF is obtained by adding the individual variances in quadrature. Applying standard error propagation to Eq. (C.1) results in $$\delta \text{CCF}=\sqrt{{\displaystyle \sum }_{l=1}^{N_m}{\displaystyle \sum }_{x=1}^{N_p}{\left(m_l{\Delta }_{xl}\right)}^2\delta f_x^2}.$$(C.4)

Substituting the expression for pixel variance: $$\delta \text{CCF}=\sqrt{{\displaystyle \sum }_{l=1}^{N_m}{\displaystyle \sum }_{x=1}^{N_p}{\left(m_l{\Delta }_{xl}\right)}^2\left(g{f}_x+g^2{\text{RON}}^2\right)}$$(C.5)

Appendix D Correlation between distortion coefficients and classical CCF activity indicators

We show the correlation between decomposition coefficients and classical CCF indicators.

For the ϵ Eri case, the strongest correlations are CON with c₁, RV with l, FWHM with g₂, and BIS with g₃. This confirms that the decomposition recovers the classical CCF activity indicators. The strongest correlation between decomposition coefficients is between l and g₃, and between c₁ and g₂.

Similarly, for the TZ Ari case, the strongest correlations are also between CON with c₁, RV with l, FWHM with g₂, and BIS with g₃. However, we see a smaller degree of correlation between FWHM and g₂, and between BIS and g₃, compared to the ϵ Eri case. The strongest correlation between decomposition coefficients is between c₀ and c₁, and between l and g₃.

Appendix E Scale mismatch between synthetic and observed RMS time series coefficients

We generate synthetic datasets using StarSim aimed at reproducing the stellar activity signals present in the observed data. However, a discrepancy remains between the simulations and observations, as illustrated by the difference in RMS between the synthetic and observed coefficient time series. In order to minimise this discrepancy, we apply z-score normalisation to both datasets so that their scales matches.

Appendix F Injection and retrievals for the null hypothesis

We perform injection and retrievals for the null hypothesis case. The detection limits are 3.99 ± 0.18 m s⁻¹ (10% precision in K) and 3.87 ± 0.18 m s⁻¹ (30% precision).

Appendix G Gaussian process kernel

In order to benchmark the stellar activity correction capabilities of CANSTAR, we also test the GP regression strategy used in Quirrenbach et al. (2022). In agreement with such work, we employ the Simple Harmonic Oscillator (SHO) kernel. This kernel is physically motivated, describing a stochastically driven, damped harmonic oscillator, which acts as an effective approximation for the quasi-periodic variability characteristic of stellar rotation and active region evolution. This kernel takes the form $$k(\tau )=S\omega Q\exp \left(-\frac{\omega \tau }{2Q}\right)\left[\cos (\eta \omega \tau )+\frac{1}{2Q\eta }\sin (\eta \omega \tau )\right],$$(G.1)

where S is the amplitude, ω is the angular frequency, and Q is the quality factor of the oscillator. ω is related to the rotational period ω = 2π/P and Q is also related to the typical coherent lifetime, l of the signal, Q = lω/2. In our analysis, we fit for the hyper-parameters utilising the celerite2 implementation (Foreman-Mackey et al. 2017).

Fig. A.1 Time series (left) and GLS periodograms (right) of the TESS photometry of TZ Ari for sectors 70 (top) and 71 (bottom). We mark in purple the stellar rotation period (1.96±0.07 d).

Appendix H CANSTAR + Keplerian solution posteriors

We show the posterior distributions from the CANSTAR + Keplerian model.

Fig. D.1 Correlations between the first 8 decomposition coefficients, colour-coded according to the absolute value of their Pearson correlation coefficient, where the exact value is printed over the image. Strongest correlations with classical CCF indicators appear at CON–c1, RV–l, FWHM–g2, and BIS–g3 (top right).

Fig. D.2 Same as Fig. D.1 but for the multi-Hermite basis applied to TZ Ari.

Fig. E.1 Histograms showing the distribution of RMS values for the coefficients derived from 10 000 StarSim simulations for ϵ Eri. The RMS values of the actual observed HARPS coefficients are shown by vertical red dashed lines. We apply z-score normalisation to bridge the gap between the time series RMS from the simulated distributions and the observed ones.

Fig. F.1 Same as Fig. 11 but for the null hypothesis case.

Fig. H.1 Posterior distribution from emcee for the CANSTAR + Keplerian solution

All Tables

All Figures

	Fig. 1 Time-averaged CCFs (templates) of ϵ Eri (red) and TZ Ari (blue) with the corresponding Gaussian model fit for ϵ Eri (orange) and the two-component Gaussian model fit for TZ Ari (cyan).
In the text

	Fig. 2 Time series (left) and GLS periodograms (right) of the RV and activity indicators derived from the CCFs for ϵ Eri. We mark in purple the stellar rotation period (11.2 d) and its first harmonic (5.6 d). The 0.1, 1, and 10% FAP levels are shown with dashed orange lines.
In the text

Fig. 3 Time series (left) and GLS periodograms (right) of the RV and activity indicators derived from the CCF for TZ Ari. The stellar rotation period (1.96 d) is marked in purple; the planetary orbital period from the literature (722.05 d; Quirrenbach et al. 2022) is marked in light blue; and the 0.1%, 1%, and 10% FAP levels are shown with dashed orange lines. The frequency axis of the periodogram show two scales split at 0.035 d−1 to better visualise long-term variability. The RVs display low-frequency signals associated with the Keplerian motion of TZ Ari b, while the CON and FWHM show lower-frequency stellar activity signals.

In the text

	Fig. 4 First six basis components for ϵ Eri (Hermite, red) and TZ Ari (multi-Hermite, blue).
In the text

	Fig. 5 Explained variance ratio (solid) and cumulative variance (dashed) per component for ϵ Eri (Hermite, red) and TZ Ari (multi-Hermite, blue).
In the text

Fig. 6 Cross-correlation function decomposition of ϵ Eri with a Hermite basis. Left: contribution of consecutive basis components to the differential CCF reconstruction. We displayed six illustrative epochs randomly selected to visualise the temporal variability of the line-shape distortions. The observed differential CCFs are shown in grey, while the reconstructed components are colour-coded according to the observation time. Middle: time series of the corresponding coefficients. Right: GLS periodograms of the coefficients. The 0.1%, 1%, and 10% FAP levels are indicated in orange, and the stellar rotation period and its first harmonic are marked in purple.

In the text

	Fig. 7 Same as Fig. 6 but for the multi-Hermite basis applied to TZ Ari. We also mark in light blue the planetary orbital period from the literature (722.05 d; Quirrenbach et al. 2022).
In the text

Fig. 8 Schematic representation of CANSTAR (bottom), a convolutional attention network designed to predict the stellar activity contribution to line shifts, κ(t), from the time series of distortion coefficients, gi(t). The input coefficients are treated as independent channels and processed through a series of one-dimensional convolutional layers that extract local temporal features. The convolutional outputs are then passed to a transformer encoder module, which consists of a multi-head self-attention part (top) with different ‘heads’, where each one of them focuses on different temporal dependencies within the data, and their outputs are combined, normalised, and propagated through additional attention layers. Finally, the resulting features are flattened and projected through a linear layer to produce the predicted κ(t) time series.

In the text

	Fig. 9 Residual relative error of ϵ Eri (top) and TZ Ari (bottom) StarSim data after correction with CANSTAR (solid line with circle marker) and an FCN (dashed line with square marker) for the noiseless case (blue) and for the case with noise equivalent to the observations (green).
In the text

Fig. 10 Top: Radial velocity time series (left) and GLS periodogram (right) of the HARPS ϵ Eri data (blue) compared to the stellar activity signal predicted by CANSTAR (green). Bottom: Residuals after subtracting the prediction (grey), with their corresponding periodogram (right). The 10%, 1%, and 0.1% FAP levels are indicated by dashed orange lines and the stellar rotation period (11.2 d), and its first harmonic (5.6 d) is marked in purple. The residuals are compatible with the ϵ Eri b solution (red; Thompson et al. 2025).

In the text

	Fig. 11 Detection limits in the ϵ Eri residuals after applying CANSTAR’s correction, shown as a function of the difference between injected and retrieved frequencies (left) and semi-amplitudes (right) for the different injected frequencies and semi-amplitudes sinusoidal signals.
In the text

	Fig. 12 Histograms showing the posterior distribution of the shared parameters between the CANSTAR and GP correction methods used to model the TZ Ari RV data.
In the text

Fig. 13 Radial velocity time series and GLS periodograms for TZ Ari showing (a) the complete multi-season dataset and (b) a detailed view of the final observing season. For each panel, the top row displays the CARMENES data (blue), the CANSTAR prediction (green), and the Keplerian solution (red). In the time series, the green curve includes the Keplerian signal to illustrate the full fit, whereas in the periodogram it represents the activity model alone. The bottom row shows the residuals after subtracting both the CANSTAR activity correction and the Keplerian model (grey). The 10%, 1%, and 0.1% FAP levels are indicated by dashed orange lines. The stellar rotation period (1.96 d) is marked with a purple dashed line, and the newly derived planetary orbital period (789.92 d) is marked in light blue (not shown in panel (b) due to its shorter time baseline).

In the text

Fig. 14 Residual time series (left) and GLS periodogram (right) after subtracting the GP + Keplerian model to the TZ Ari RV data (grey). The 10%, 1%, and 0.1% FAP levels are indicated by dashed orange lines. The stellar rotation period (1.96 d) is shown in purple, the newly derived planetary orbital period (789.92 d) is in light blue, and the two signals exceeding the 1% FAP level (37.2 d and 41.4 d) are marked in black. These likely correspond to stellar activity given the significant (40.18 d) periodicity observed in the FWHM and g2 indicators.

In the text

Fig. 15 Top: radial velocity time series (left) and GLS periodogram (right) of the HARPS ϵ Eri data (blue) compared to the stellar activity signal predicted by CANSTAR (green) and Perger et al. (2023, magenta). We note that the Perger et al. (2023) did not provide uncertainties to their NN predictions. Bottom: Residuals after subtracting the prediction with their corresponding periodogram (right). The 10%, 1%, and 0.1% FAP levels are indicated by dashed orange lines.

In the text

	Fig. A.1 Time series (left) and GLS periodograms (right) of the TESS photometry of TZ Ari for sectors 70 (top) and 71 (bottom). We mark in purple the stellar rotation period (1.96±0.07 d).
In the text

	Fig. D.1 Correlations between the first 8 decomposition coefficients, colour-coded according to the absolute value of their Pearson correlation coefficient, where the exact value is printed over the image. Strongest correlations with classical CCF indicators appear at CON–c1, RV–l, FWHM–g2, and BIS–g3 (top right).
In the text

	Fig. D.2 Same as Fig. D.1 but for the multi-Hermite basis applied to TZ Ari.
In the text

	Fig. E.1 Histograms showing the distribution of RMS values for the coefficients derived from 10 000 StarSim simulations for ϵ Eri. The RMS values of the actual observed HARPS coefficients are shown by vertical red dashed lines. We apply z-score normalisation to bridge the gap between the time series RMS from the simulated distributions and the observed ones.
In the text

	Fig. F.1 Same as Fig. 11 but for the null hypothesis case.
In the text

	Fig. H.1 Posterior distribution from emcee for the CANSTAR + Keplerian solution
In the text