The CARMENES search for exoplanets around M dwarfs: Stable radial-velocity variations at the rotation period of AD~Leonis -- A test case study of current limitations to treating stellar activity

Context: A challenge with radial-velocity (RV) data is disentangling the origin of signals either due to a planetary companion or to stellar activity. In fact, the existence of a planetary companion has been proposed, as well as contested, around the relatively bright, nearby M3.0V star AD Leo at the same period as the stellar rotation of 2.23d. Aims: We further investigate the nature of this signal. We introduce new CARMENES optical and near-IR RV data and an analysis in combination with archival data taken by HIRES and HARPS, along with more recent data from HARPS-N, GIANO-B, and HPF. Also, we address the confusion concerning the binarity of AD Leo. Methods: We consider possible correlations between the RVs and various stellar activity indicators accessible with CARMENES. We applied models within a Bayesian framework to determine whether a Keplerian model, a red-noise quasi-periodic model using a Gaussian process, or a mixed model would explain the observed data best. We also exclusively focus on spectral lines potentially associated with stellar activity. Results: The CARMENES RV data agree with the previously reported periodicity of 2.23d, correlate with some activity indicators, and exhibit chromaticity. However, when considering the entire RV data set, we find that a mixed model composed of a stable and a variable component performs best. Moreover, when recomputing the RVs using only spectral lines insensitive to activity, there appears to be some residual power at the period of interest. We therefore conclude that it is not possible to determinedly prove that there is no planet orbiting in synchronization with the stellar rotation given our data, current tools, machinery, and knowledge of how stellar activity affects RVs. We do rule out planets more massive than 27M_E (=0.084M_J). We also exclude any binary companion around AD Leo with Msini>3-6M_J on orbital periods<14yr.


Introduction
In the pursuit for Earth-like planets, modern spectrographs are pushing the limits by reaching m s −1 level precision or even tens of cm s −1 (e.g., ESPRESSO; Pepe et al. 2021), which is needed. However, the measurements begin to succumb to unwanted signals for planet searches. Intrinsic stellar variability in the form of dark spots, bright plages, and flares can produce radial-velocity (RV) variations that can conceal true planetary signals, or even masquerade as a fake planet which then can be effectively modeled with a Keplerian orbit. To help mitigate these stellar-activity-induced RV signals, a number of procedures are commonly put in place, such as different modeling approaches, smart data collection strategies, and extraction of particular spectral lines.
Statistical techniques such as Gaussian process (GP) regression have been used by treating stellar activity behavior as a quasi-periodic signal (e.g., Haywood et al. 2014;Rajpaul et al. 2015;Jones et al. 2017). This approach may sometimes lead to better precision and accuracy of planetary parameters, especially those in the lower-mass regime (e.g., Amado et al. 2021). Likewise, these signals can be wavelengthdependent, usually with the amplitude decreasing in the redder regime of the spectrum, but still containing some residual effect depending on the star-spot configuration and temperature difference (Reiners et al. 2010, and references therein). For this reason, there is a push for higher-precision instruments covering the red and near-IR wavelength range such as CARMENES 1 (Quirrenbach et al. 2014(Quirrenbach et al. , 2018, GIARPS 2 (Claudi et al. 2017), HPF 3 (Mahadevan et al. 2012(Mahadevan et al. , 2014, IRD 4 (Tamura et al. 2012;Kotani et al. 2014Kotani et al. , 2018,   (Seifahrt et al. 2018(Seifahrt et al. , 2020, and SPIRou 6 (Donati et al. 2020).
Focusing on certain spectral lines as activity indicators sensitive to chromospheric (e.g., Hα, Ca ii infrared triplet) or photospheric (e.g., TiO) effects on active M dwarfs can oftentimes be successful in determining the star's rotational period (see Fig. 11 in Schöfer et al. 2019). Even then, there is no universal approach because it is still not unique as to which activity indices do peak at the rotational period and under what conditions, as pointed out by Lafarga et al. (2021). Moreover, efforts for identifying which spectral lines, in general (i.e., not focusing on alreadyknown specific lines), seem to be more activity-sensitive than others have been fruitful for a selection of G-K dwarfs (e.g., Wise et al. 2018;Dumusque 2018;Lisogorskyi et al. 2019;Ning et al. 2019;Cretignier et al. 2020; Thompson et al. 2020). Such an approach proves to be rather challenging for M dwarfs, where the spectra contain a forest of blended lines, making it almost impossible to find the continuum (Merrill et al. 1962;Boeshaar 1976;Kirkpatrick et al. 1991;Alonso-Floriano et al. 2015). This approach, however, seems to have been successful for stars that show a clear stellar activity impact on the RVs (e.g., EV Lac, Lafarga priv. comm.). 1 Calar Alto high-Resolution search for M dwarfs with Exo-earths with Near-infrared and optical Échelle Spectrographs, http://carmenes. caha.es 2 GIAno & haRPS-n, https://www.tng.inaf.it/news/2017/ 04/04/giarps/ 3 The Habitable Zone Planet Finder, https://hpf.psu.edu/ 4 InfraRed Doppler instrument 5 Red-optical, high-resolution spectrograph with focus on mid-to late-type M dwarfs, https://www.gemini.edu/instrumentation/ maroon-x/ 6 A near-infrared spectropolarimeter/velocimeter, https://spirou. omp.eu/ We turn our attention specifically to the active mid-type M dwarf AD Leo, a star whose stellar rotation period of 2.23 d presents itself both in photometry and RVs Tuomi et al. 2018;Carleo et al. 2020;Robertson et al. 2020). Despite its strong flaring activity manifesting itself at many wavelengths (Buccino et al. 2007;Rauer et al. 2011;Tofflemire et al. 2012;Vidotto et al. 2013), AD Leo has been included in a number of studies addressing the existence of planets orbiting this star. Tuomi et al. (2018), referred to as T18 hereinafter, first suggested that a planet may be orbiting AD Leo in a 1:1 spinorbit resonance since it proved to be difficult to simultaneously explain both the photometry and RV measurements using a variety of star-spot scenarios. Furthermore, they claimed that the RV measurements were time-and wavelength-independent, and the putative planet exhibited a semi-amplitude of ∼19 m s −1 . Despite some evidence for solely stellar activity behavior, T18 concluded that AD Leo is an active M dwarf hosting a hot Jupiter (≥ 0.2M Jup ) in a 1:1 spin-orbit resonance.
However, the existence of this hot Jupiter around AD Leo was challenged. Carleo et al. (2020, referred to as C20 hereinafter) investigated the 2.23 d signal and obtained observations with GIARPS at the 3.56 m Telescopio Nazionale Galileo (Claudi et al. 2017). Even though the 2.23 d signal was persistent, the amplitude heavily diminished as a function of wavelength and of time. Simultaneous photometric data from STELLA showed a shift of ∼0.25 in phase (∼0.6 d) in comparison to the HARPS-N RV curves. Therefore, C20 disputed the argument posed by T18, concluding that the RV modulation is not compatible with a planetary companion. Shortly after, the conclusions by Robertson et al. (2020), referred to as R20 hereinafter, were also in line with C20, as they likewise observed a decrease in amplitude between the two observing seasons using HPF spectroscopic data. At the time of submission, additional data from SPIRou and SOPHIE 7 also indicated no evidence for a corotating planet (Carmona et al. in prep.).
This concept of such a hot Jupiter around an M dwarf in a synchronized rotation is unique. Many studies addressing the correlation between the orbital period (P orb ) and the stellar rotation period (P rot ) have thus far been focused on more solarlike transiting host stars rather than M dwarfs, using Kepler data (Borucki et al. 2010). Findings from McQuillan et al. (2013) and Walkowicz & Basri (2013) suggest a clear absence of close-in planets (P orb ≲ 2-3 days) around rapidly rotating stars (P rot ≲ 5-10 days), where planets with shorter periods were nearly synchronous (P orb ∼ P rot ). Teitler & Königl (2014) proposed that the cause was orbital decay of close-in planets by their host stars. Besides, with regards to M dwarfs, Newton et al. (2016) showed that typical stellar rotation periods would often match the periods of planets in their habitable zone. Regardless, there is a dearth of close-in massive planets around M dwarfs where only a handful of such are known, such as NGST-1 b (Bayliss et al. 2018) or TOI-519 b (M ≲ 14 M Jup ; Parviainen et al. 2021). This is not a result of observational bias as these types of planets would be simple to detect through transits and precise RVs. AD Leo was even used as a test case to check the effects of XUV irradiation on planet atmospheres (Chadney et al. 2016), considering that the ionsphere can play an influential role in regulating the stability of upper atmospheres on giant planets. Therefore, understanding the case of AD Leo could aid in finding the key methods for other potential planet detections that suffer from ambiguities due to 1:1 spin-orbit resonances.
In our focused study using optical and near-IR RV measurements, we claim that one cannot prove (or disprove) the suggested planetary signal in the AD Leo system, as the problem is completely degenerate and we are limited by our current state-of-the-art machinery and measurements. We organize the paper as follows. In Sect. 2, we first introduce the active M dwarf AD Leo and investigate its presumed binarity for the first time. All available RV data used for the analysis are presented in Sect. 3. We turn our focus on the chromaticity of the RV signal for the CARMENES data in Sect. 4. Still concentrating on the CARMENES data, in Sect. 5 we perform a spectral line analysis for the target. Then combining all available spectroscopic data for the first time, we explore various modeling techniques within the Bayesian framework in Sect. 6. Finally, in Sect. 7, we present a discussion on our results as well as suggestions to break the degeneracy of this situation. We summarize our conclusions in Sect. 8.

Stellar parameters
AD Leo (GJ 388), an M3.0 V star at a distance of slightly less than 5 pc and with V ∼ 9.5 mag, is one of the closest and brightest M dwarfs. Already tabulated in the Bonner Sternverzeichnis by Argelander (1861), AD Leo has been the subject of numerous investigations in the last century (e.g., Abell 1959;Engelkemeir 1959;Lang et al. 1983;Saar & Linsky 1985;Hawley & Pettersen 1991;Hawley et al. 2003;Osten & Bastian 2008;Hunt-Walker et al. 2012).
In Table 1, we list the stellar properties of AD Leo. In particular, we tabulate equatorial coordinates, proper motions, and parallax from the Gaia Early Data Release 3 (EDR3; Gaia Collaboration et al. 2021) and absolute RV from Lafarga et al. (2020, uncorrected of gravitational redshift for consistency with the previous literature), from which we recompute Galactocentric space velocities as in Cortés-Contreras (2016). The spectral type of Alonso-Floriano et al. (2015) superseded previous determinations (e.g., Johnson & Morgan 1953;Bidelman 1985;Stephenson 1986;Keenan & McNeil 1989), while the photosphere parameters of Passegger et al. (2019) match previous CARMENES publications and are similar, but not identical, to those of Rojas-Ayala et al. (2012), Lépine et al. (2013), Gaidos et al. (2014), or Mann et al. (2015. With the effective temperature of Passegger et al. (2019), the bolometric luminosity of Cifuentes et al. (2020), and the Stefan-Boltzmann law we derived the stellar radius and, with the radius-mass relation of Schweitzer et al. (2019), the stellar mass. For compiling the most precise parameters of the activity indicators, we used the SVO Discovery Tool 8 . The projected rotational velocity v sin i was computed by us exactly as in Reiners et al. (2018), but on the newest CARMENES template spectra (Sect. 3).
As first reported by Montes et al. (2001), the Galactocentric space velocity of AD Leo is consistent with it belonging to the Galactic young disk (Leggett 1992). Later, López-Santiago et al. (2010) and Klutsch et al. (2014) proposed AD Leo as a candidate member of the Castor moving group, in agreement with our latest kinematic data. The age of the Castor moving group, of about 300-500 Myr (Barrado y Navascués 1998; Mamajek et al. 2013), is consistent with age determinations for AD Leo by Shkolnik et al. (2009), Brandt et al. (2014, and Meshkat et al. (2017).
Such a young age partly explains the flares frequently observed in AD Leo. The star has been known to exhibit activity 8 http://sdc.cab.inta-csic.es/SVODiscoveryTool/ 2006) 9 . Frequent flaring activity was further observed during an extreme ultraviolet (EUV) 1-month monitoring of AD Leo (Sanz-Forcada & Micela 2002). Later on, emission in the Xray and radio regimes has also been observed (Gurzadyan 1971;Robinson et al. 1976), and which has shown bursting radio emission at GHz frequencies and below (Osten & Bastian 2008;Villadsen & Hallinan 2019). Muheki et al. (2020) and Namekata et al. (2020) have presented the most recent analyses on highresolution optical spectroscopy and X-ray observations of flares on AD Leo. The young age of AD Leo also explains the moderately large rotational velocity and short rotational period, of 2.23 d, as well as X-ray, Ca ii H&K, and Hα emission (see references in Table 1). In addition, this star presents large RV variations: it shows a standard deviation larger than 20 m s −1 , (∼ 1.5× median absolute deviation RV as in Tal-Or et al. 2019 andGrandjean et al. 2020), which can be connected to stellar activity, the presence of a planetary companion, or both. For the derivation of log R ′ HK , we used the method described by Perdelwitz et al. (2021), albeit with a slight modification of the k1 passband used for normalization. A total of 316 archival spectra from ESPaDOnS (232 spectra), HARPS (58 spectra), HIRES (24 spectra) and NARVAL (2 spectra) were analyzed, yielding a mean value of log R ′ HK = −3.97 ± 0.05, which is in agreement with the result published by Astudillo-Defru et al. (2017) based on the HARPS spectra alone. We compute the stellar inclination, i, to be ∼ 13 • from 10 6 MCMC realizations using v sin i, R, and P rot as provided in the stellar parameters table. Previous papers quoted values of ∼ 20 • ) and ∼ 15 • (Tuomi et al. 2018).

Photometric rotational period
T18 took a closer look at the available photometry from ASAS-North and ASAS-South, as well as from short-cadence MOST observations (see T18 for photometry references). To summarize, the ASAS photometry shows a long-term periodicity of around 4070 d, most likely due to a stellar activity cycle. During the brightness minimum of this cycle, the signal at the rotation period of 2.23 d seems to disappear turning up only during the brightness maximum, in which the MOST observations were taken. These short-cadence MOST data, taken over the course of 9 d , showed fluctuations with periodicity 2.23 d, but also demonstrated behavior of slight phase shifts and amplitude and period variations, indicating that the stellar surface should be experiencing rapid evolution. Other observations taken with different instruments are in agreement with the 2.23 d period (see Table A.1).
AD Leo was additionally observed in three of the TESS sectors 10 , particularly in sector 45 (camera #3 CCD #4; 6 November 2021 to 02 December 2021), sector 46 (camera #2 CCD #2; 02 December 2021 to 30 December 2021), and in sector 48 (camera #1 CCD #4; 28 January 2022 to 26 February 2022). All sectors constitute short cadence 20-second and 2-minute integrations. Unfortunately, the pixel column was saturated for sectors 45 and 46 leading to unusable data. The sector 48 data were still valuable for determining a periodic variation. In doing so, we removed strong flare events (up to 275 ppt) and any outliers on the 20-second integration data which were binned every 10 points for computational reasons. After applying a sinusoid model, the resulting amplitude and P rot are ∼800 ppm and P rot of 2.2304 ± 0.0014 d, respectively.

Hypothetical binary
There is some confusion in the literature regarding the hypothetical multiplicity of AD Leo. First of all, in spite of being listed as "γ Leo C" in some astronomical databases (e.g., SIMBAD), AD Leo is not a wide (ρ ∼ 5.7 arcmin) physical companion to Algieba, which is a binary system of two G-and K-type giant stars visible to the naked-eye located eight times further (γ 01 Leo + γ 02 Leo; van Leeuwen 2007;Gaia Collaboration et al. 2021).
Next, the Washington Double Star (WDS; Mason et al. 2001) catalog tabulates AD Leo as a close binary candidate. The presence of a hypothetical companion around AD Leo was first inferred by Reuyl (1943) from measurements of photographic plate images. He suggested an eccentric (e = 0.6) and close-in (a ∼ 0.54 au, projected angular separation of ρ ∼ 0.11 arcsec) orbit with an orbital period of 26.5 years, making the companion a brown dwarf (M ≈ 0.032M ⊙ ). Later, van de Kamp & Lippincott (1949) found indications of an astrometric trend from photographic plates that did not fit those orbital elements. Following up two decades later, Lippincott (1969) found an ambiguous deviation from linear proper motion, thus deciding that it was inconclusive to determine whether there could be a variable proper motion due to a companion.
With speckle imaging at 7800 Å at the 3.6 m Canada-France-Hawai'i Telescope, Balega et al. (1984) resolved a companion candidate to AD Leo at ρ = 0.078 ± 0.010 arcsec (r ≈ 0.39 au) and position angle of θ = 39 ± 4 deg in 1981, which received the WDS name 10200+1950 and discoverer code BAG 32. Two additional measurements in 1983 were provided by Balega & Balega (1985) at ρ ∼ 0.11 arcsec, but with θ ∼ 330 deg. Afterwards, with lucky imaging in the I band at the 1.5 m Telescopio Carlos Sánchez, Cortés-Contreras et al. (2017) resolved a candidate source, about 2.0 mag fainter than AD Leo, at ρ = 0.195 ± 0.061 arcsec (r ≈ 0.97 au) and θ = 23.8 ± 3.7 deg in 2012. However, it fell on the first Airy disk, and they were not able to detect it again in observations in 2015 with the same instrument setup. From ∆I ≈ 2.0 mag and the projected physical separation of about 0.97 au, Cortés-Contreras et al. (2017) estimated an M5 V spectral type for the companion candidate and an orbital period of P ∼ 1.5 yr, which might be the source discovered by Balega et al. (1984). Because of this hypothetical companion, AD Leo was initially discarded from the CARMENES guaranteed time observation target list (Caballero et al. 2016;Reiners et al. 2018).
Further attempts at resolving the companion proved to be unsuccessful, such as additional speckle imaging by Docobo et al. (2006). With adaptive optics at 8 m-class telescopes, Daemgen et al. (2007) and Brandt et al. (2014) imposed very strict upper limits to the presence of companions, of ∆K s ∼ 7.8 mag at ρ ≥ 0.5 arcsec and ∆H ∼ 10.9 mag at ρ ≥ 1.5 arcsec, respectively. From Fig. 7 in Daemgen et al. (2007), their Altair/Gemini North observations discarded any companion of ∆K s ∼ 2.0 mag at ρ ≥ 0.1 arcsec. No trace of companions is seen either in archival images obtained in 2005 with the Hubble Space Telescope (ACS/HRC in F330W) and in 2018 with the Very Large Telescope (SPHERE in H) 11 . Therefore, the findings of Balega et al. (1984), Balega & Balega (1985), and D. Kossakowski et al.: AD Leo: stellar activity in radial velocity Cortés-Contreras et al. (2017) are most probably not related to a real physical companion, unless the orbital motion can explain the later nondetections. Furthermore, the Gaia EDR3 renomalized unit weight error (RUWE) for AD Leo is 1.15, below the critical threshold of 1.40, the along-scan observations (astrometric_gof_al) is 3.22, and the excess noise of the source (astrometric_excess_noise) is 0.16, all of which indicate that it is most likely a single source (Lindegren et al. 2018). To rule out any wide companions, we searched for objects with common parallax and proper motions up to a projected physical separation of 100 000 au, as in Caballero et al. (2022), and found no hints of any wide potential companion. Figure 1 illustrates all the mentioned findings from the literature. The binary issue was not addressed by T18, C20, or R20. The absence of any close companion to AD Leo is further investigated with our spectroscopic data in Sect. 2.4.

Evidence for single star in spectroscopic data
As a first approach to search to address the potential stellar binarity and look for traces of it, we combined all the RV data (presented in Sect. 3) to look for long-term trends. However, this method presents its challenges as there is no temporal overlap between the older and newer data sets and each instrument has its own, unknown zero-point offset. We nonetheless attempted a grid search for a Keplerian signal for which we stepped through the period, amplitude, eccentricity and argument of periastron parameter space, where each instrument also had its own offset and jitter term 12 (e.g., Baluev 2009). One then adapts the values with the best log likelihood. Unfortunately, the lack of temporal overlap led to strong ambiguities and degeneracies. We therefore did not find any conclusive periodicity or indicative linear trend in this RV analysis.
We continued to investigate whether there is evidence for the presence of a companion within our CARMENES spectra that could affect the analysis presented later in this work. We computed a 1D cross-correlation function (CCF) with a binary mask over a large RV range and did not find any hint of a clear secondary peak, neither in the VIS nor the NIR data. To be certain, we also ran todcor (Zucker & Mazeh 1994), which computes a 2D CCF to get the RVs of the two components simultaneously, and the results showed no evidence for a companion as well. The 2D CCF method with todcor is well suited in the case of double-lined binaries, but since the secondary signal seems to be either too weak if expected to be ∼10 km s −1 away from the primary, or too hidden if expected to be very close to the primary, we concluded that no secondary heavily distorts the CCF profile and, therefore, would not cause any noticeable effects on the CCF parameter values. Therefore, we first discard the presence of a nearly equal brightness double-lined spectroscopic binary from the todcor analysis. Secondly, based on our RV time series analysis, we rule out companions with masses greater than 3-6 M Jup using orbital periods of less than ∼2-14 yr and amplitudes 3-5 greater than the rms of the earlier data set instruments (i.e., HARPS, HIRES), while assuming sin i 0. We therefore conclude that there is no evidence for a stellar binary or brown dwarf companion of AD Leo within this parameter regime. 12 Jitter terms are added in quadrature to the given error bars for each respective instrument.

Bursting radio emission from the AD Leo system
Radio bursts from AD Leo have been known for quite a long time, and have been ascribed to (coherent) plasma emission (Stepanov et al. 2001;Osten & Bastian 2008). More recently, Villadsen & Hallinan (2019) detected both short-duration (seconds to minutes) and long-duration (30 minutes or more) bursts of AD Leo, using ultrawide-band VLA observations (in the ranges 0.2-0.5 GHz and continuous frequency coverage between 1 and 6 GHz) in several epochs between 2013 and 2015. The timescale of the short-duration events was consistent with the duration of "space-weather" events, such as the impulsive phase of a flare or a coronal mass ejection (CME) crossing the corona. The long-duration bursts, which lasted hours, and possibly extended for up to days or even longer times, between observing epochs (similarly to the case of Proxima Centauri; see, e.g., Pérez-Torres et al. 2021 and references therein) require an ongoing electron acceleration mechanism during the burst. Candidate acceleration processes are found within the paradigms of solar radio bursts and periodic radio aurorae produced by ultracool dwarfs and planets.
The emission mechanism responsible for the aurorae from stars and planets alike is the electron cyclotron maser (ECM) instability (Melrose & Dulk 1982), whereby plasma processes within the star (or planet) magnetosphere generate a population of unstable electrons that amplifies the emission. The characteristic frequency of the ECM emission is given by the electron gyrofrequency, ν G = 2.8 B MHz, where B is the local magnetic field in the source region, in Gauss. ECM emission is a coherent mechanism that yields broadband (∆ ν ∼ ν G /2), highly polarized (sometimes reaching 100%), amplified nonthermal radiation. The long-duration bursts seen for AD Leo by (Villadsen & Hallinan 2019) between 1 and 6 GHz are characteristic of ECM emission, implying electron densities ≲ 10 9 cm −3 , and originate in regions with field strengths of 0.36-2.1 kG (in fundamental emission) or 0.18-1.1 kG (if second harmonic). Since the surface magnetic field of AD Leo is estimated to be of ≃ 3.36 kG (Reiners et al. 2022), the above field strengths are likely to occur relatively close to the star, at a height less than about one stellar radius.
Auroral cyclotron maser emission is powered by the acceleration of confined electrons with energies between ∼10 keV and up to 1 MeV energies. For substellar objects, the currents could be driven by the breakdown of rigid corotation of magnetospheric plasma with the object's magnetic field, for example, due to the interaction between a rotating magnetosphere and the interstellar medium. Corotation breakdown observed in a number of ultracool (UCD) dwarfs with periods less than about 3 hr seem to be rotation powered, with radio powers of up to a few times 10 13 erg s −1 Hz −1 (e.g., Turnpenney et al. 2017). AD Leo has a rotation period of ∼54 hr. Assuming coronal parameters similar to radio-loud UCDs, any corotation breakdown in AD Leo will generate a radio power that is ∼ 6 × 10 11 erg s −1 Hz −1 , or less than about 100 µ Jy. Therefore, corotation breakdown, while is likely to contribute to the overall radio power observed from AD Leo, cannot account for the observed radio bursts, which reach peaks of tens of mJy at GHz frequencies (Villadsen & Hallinan 2019). A discussion of further potential mechanisms, such as star-planet interaction, to explain the measurements is continued in Sect. 7.2.3.

Spectroscopic data
The Doppler data that are used in our analysis are described below. All data and spectrograph descriptions are summarized in Table 2. Data used for the analysis are displayed in Fig. 2, along with the generalized Lomb-Scargle (GLS) periodograms (Zechmeister & Kürster 2009) of the data in Fig. 3.

CARMENES
The CARMENES instrument is located at the 3.5 m telescope at the Calar Alto Observatory in Spain. It is a dual arm instrument to produce RV measurements at optical (VIS) and near-infrared (NIR) wavelengths (Quirrenbach et al. 2014(Quirrenbach et al. , 2018 Table 2 for the instrument specifications). We obtained 26 spectra for AD Leo (Karmn J10196+198) over the time span of 25 days in 2018, notably 12 years after the majority of the HIRES and HARPS data were taken. On some nights, we observed AD Leo multiple times to better sample the short period of 2.23 d. We call this subset of data "VIS1" and "NIR1" given that AD Leo was observed again later in 2020 as part of a DDT 13 program, for which 20 data points were obtained over ∼3 d, which we call subset "VIS2" and "NIR2". The RVs and several activity indicators (e.g., CRX -chromatic index, dLW -differential line width) from both subsets were first extracted with serval (Zechmeister et al. 2018), so that they are corrected for barycentric motion, secular acceleration, and instrumental drift. The final RVs we use had nightly zero-points applied (Tal-Or et al. 2019;Trifonov et al. 2020). Figure 4 displays the CARMENES data and the best respective models (detailed in Sect. 6). For calculating the CCF parameters within CARMENES, weighted binary masks, which depend on spectral type and v sin i, are produced by coadding spectra corrected for tellurics & RV shifts and then selecting pronounced minima (Lafarga et al. 2020). From there, the CCF 13 Director's Discretionary Time parameters, namely, the bisector span (BIS), contrast, and fullwidth half-maximum (FWHM), are obtained for both the VIS and NIR channel.

Archival RV data
HARPS-South. Spectroscopic data for AD Leo are in the public archive from the High-Accuracy Radial velocity Planet Searcher (HARPS, Mayor et al. 2003), a high-resolution échelle spectrograph located at the ESO 3.6 m telescope at La Silla Observatory in Chile. We retrieved a total of 52 spectra that were taken over a span of ∼4500 d. Due to a fiber-upgrade intervention (Lo Curto et al. 2015), we considered the HARPS data as coming from two separate instruments, before and after this upgrade. After the intervention, only five data points were taken close to one another (within 2 d) and much later than the earlier data. These data also have a very low root-mean-square (RMS) deviation from their mean. We did not consider them for the further analysis since they could fit anywhere in a model with a large offset and therefore do not provide new insight. Thus, for the rest of the analysis, we used the 47 HARPS spectra from before the intervention. Additionally, the majority of data (i.e., 33 points) were taken over the course of ∼115 d (January-May 2006). The spectra were first processed using serval ) to obtain the RVs and helpful stellar activity indicators (e.g., CRX, dLW, Hα, see Zechmeister et al. 2018, for more explanation). Nightly zero point corrections provided by Trifonov et al. (2020) were applied to the RVs 14 and are used in our analysis (Table A.

4).
HIRES. We collected 43 spectra from the high-resolution spectrograph HIRES (Vogt et al. 1994) mounted on the 10 m Keck-I telescope located at the Mauna Kea Observatory in Hawai'i, which has been in service since 1994. The data were taken over the course of 4500 days. A majority of them (i.e., 22 points) were taken in 2005/2006 over 120 d and, moreover, overlap with the higher cadence HARPS data (see Fig. 2). Butler et al. (2017) released a large RV database of 64 480 observations for a sample of 1699 stars, which was later reanalyzed by Tal-Or et al. (2019) for minor, though significant systematic effects, such as an RV offset due to the CCD upgrade in 2004, long-term drifts, and slight intra-night drifts. Therefore, we continued with the corrected HIRES RVs provided by Tal-Or et al. (2019), presented in Table A.5. Data from HIRES are also specifically addressed further under Sect. 5.2.

Additional RV data from the literature
Radial velocities are taken directly from C20 and R20, as described in the following paragraphs, and are used for the analysis as well.
GIARPS. We include spectroscopic data taken in GIARPS mode (Claudi et al. 2017), where high resolution spectroscopic measurements are obtained simultaneously with HARPS-N (Cosentino et al. 2012, extracted with the TERRA pipeline) and with GIANO-B (Oliva et al. 2006) reduced with the on-line DRS pipeline and the off-line GOFIO pipeline. Both instruments are located at the 3.58 m Telescopio Nazionale Galileo (TNG) at the Roque de los Muchachos Observatory in La Palma, Spain. The GIARPS mode is similar to CARMENES in the sense it also addresses potential variations of a signal amplitude over a wide wavelength range. There were two runs of HARPS-N, with only the first run having simultaneous GIANO-B data; four of the HARPS-N data overlap with the CARMENES data. To stay consistent with C20, we also consider two separate runs for HARPS-N and GIANO-B, and designate them as "HT1", "HT2", "G1", and "G2" for our analysis. Due to the high uncertainty on the data points from the GIANO-B instrument, we do not include them for our analysis.
HPF. Our last data set comes from the Habitable-zone Planet Finder (HPF), a stable NIR Doppler spectrograph that is designed to reach 1-3 m s −1 RV precision for M dwarfs with the help of wavelength calibration via a custom NIR laser frequency comb (Mahadevan et al. 2012(Mahadevan et al. , 2014. The spectrograph is in- Notes. (a) A subset is treated effectively as an independent data set.  stalled at the 10 m Hobby-Eberly Telescope at McDonald Observatory in Texas. The HPF data are of high quality providing an RV precision of 1.5 m s −1 on AD Leo. A total of 35 HPF RVs were obtained, five of which during HPF commissioning and 30 afterwards, namely "HPF1" and "HPF2". The HPF data in this paper are tabulated in R20. These data overlap with optical RV data from HARPS-N and show also an amplitude decrease between observing seasons.

Wavelength dependence of RV signal in CARMENES data
Signals in RV measurements induced by dark spots corotating on the stellar surface are expected to be more pronounced in the bluer wavelength regime for M-dwarf stars (e.g., Desort et al. 2007;Reiners et al. 2010;Mahmud et al. 2011;Sarkis et al. 2018). In contrast, a true planetary signal should have an amplitude and phase that are consistent both in time and as a function of wavelength.

Chromatic index
The chromatic index is a photospheric activity indicator that measures the RV-log λ correlation where a straight line is fit to the RV values computed from individual échelle orders as a function of log λ Tal-Or et al. 2018 Different slopes at different wavelengths are not necessarily unexpected. In principle, this depends on which mechanism is dominating, namely either the star-spot temperature contrast or Zeeman broadening, where a more negative RV-CRX correlation suggests the former to be predominant. This is also equivalent to the amplitude of the RVs due to stellar features decreasing with longer wavelengths (Sect. 4.2). The fact that the slopes of CARMENES and HARPS are in contradiction could simply indicate that different mechanisms are prevalent. Likewise, taking into consideration that we are probing very different spectral lines at different wavelengths may introduce a trend with wavelength that is not yet well understood. When considering only the orders of CARMENES VIS that overlap with the HARPS wavelength range, the correlation would still stay the same. Here, the wavelength range considered is now equivalent, nonetheless the signs of the slopes are inconsistent. No firm conclusion can be drawn, though this could be a result of different mechanisms dominating during different time periods of long-term stellar activity.
Furthermore, the correlation plot for the CARMENES VIS demonstrates a "closed-loop" (circular) behavior (see Sect. 7.2.5 for a discussion). Removing the CRX trend on just the CARMENES data results in two-fold decrease in the rms of the corrected RVs (from 18.5 m s −1 to 9.2 m s −1 for first season and from 15.6 m s −1 to 9.3 m s −1 for the second). This subtrac- . Radial-velocity times series and phase-folded plots for the CARMENES VIS1, NIR1, VIS2, and NIR2 data (top to bottom) with residuals from the best-fit model. In some cases, especially in the panels with the VIS1 and VIS2 data, the uncertainties are too small to be seen in the plots. The best-fit models are overplotted (gray line) along with the 68% posterior bands gray shaded regions).
Article number, page 9 of 27 A&A proofs: manuscript no. output tion was not performed on the CARMENES NIR data since the Pearson-r coefficient was not large enough (Jeffers et al. 2022).

Wavelength segments
The CARMENES VIS channel records 55 échelle orders, 42 of which are considered when computing the RV measurement via a weighted mean . We chose to combine these orders into four wavelength segments where each segment consists of ten (or 11) échelle orders in order to preserve some precision, and the RVs are then recomputed for each respective wavelength coverage. Similarly, the CARMENES NIR channel has 28 RV orders over the Y, J, and H photometric bands. Due to telluric contamination, especially in the J band , only a selected few orders are considered (Bauer et al. 2020). We then use two wavelength segments, one for the Y and another for the H band, consisting of 12 and seven individual orders, respectively. For this wavelength segment analysis, we consider just the season one CARMENES data because the data from season two exhibit an amplitude decrease (see Sect. 6.3), and likewise, breaking up the orders would introduce noise given the few data points over a short time span. Each wavelength segment is treated as an individual data set. We fit a simple sinusoid to obtain K, the semi-amplitude of the signal and 1σ errors. When doing this, the semi-amplitude clearly decreases with increasing wavelength in the optical regime, but then reaches a plateau when continuing in the near-IR (Fig. 5). This behavior of decreasing but then constant RV semi-amplitude is in agreement with Reiners et al. (2010) for a dark spot on the surface of an M dwarf; specifically, a spot covering 1.5% of the projected surface, with a temperature 200 K cooler than the star (assumed T eff = 3700 K), and a stellar v sin i of 5 km s −1 , in line with AD Leo's stellar parameters 15 (Table 1). Simulations show a linear relation between spot coverage and RV semi-amplitude for low spot coverage values. A large spotstar temperature difference (T spot = 2 3 T eff ), in comparison, would not cause a notable semi-amplitude dependency as a function of wavelength; however, it is not likely for cooler stars, such as AD Leo, to have large spot-star temperature differences (Bauer et al. 2018).
Likewise, the HARPS instrument covers 72 spectral orders, of which 61 produce reliable RVs after being processed by serval (signal-to-noise is too low for the others). Similarly, we computed six wavelength chunks with ten spectral orders each (11 for the reddest chunk) and followed the same methodology as for the CARMENES wavelength chunks. Reiners et al. (2013) and T18 performed similar analyses but differed in interpretation. T18 suggested there is no dependence on wavelength, whereas Reiners et al. (2013) claimed otherwise and mentioned that this is the first case of a known active star with increasing amplitude with wavelength. Here, we find a slight positive incline, which is plotted for comparison to the CARMENES data in Fig. 5. The positive slope can be interpreted as being due to the Zeeman effect, which has the opposite effect compared to a spot-temperature difference where the RV amplitude is predicted to increase for redder wavelengths . For the overlapping wavelengths, the amplitudes of the wavelength chunks do not particularly agree which can simply be an artifact that HARPS data were taken at a time where AD Leo 15 We adopt the higher value of v sin i (i.e., 5 km s −1 ) rather than 2 km s −1 (as also computed in Reiners et al. 2010) because AD Leo is relatively active for its v sin i.

Identifying individual spectral lines affected by stellar activity
We further investigate the possible effect of stellar activity on the spectral lines themselves. We follow an approach similar to Dumusque (2018), which has proven to be fruitful for a selection of G-K dwarfs whose RVs are dominated by activity. Essentially, we use the CCF technique, where we obtain a binary mask that contains all available, reliable lines for a spectrum (following Lafarga et al. 2020). We then compute individual RVs of the few thousands of spectral lines we have identified, and obtain an RV time series for each line. To classify the lines according to their sensitivity to activity, we correlate their RVs to an activity indicator obtained from the same spectra (such as the CRX, BIS, or the total RV). We then select a subsample of spectral lines that are least affected by stellar activity (those that do not show a strong correlation) and recompute the RVs using this subsample to mitigate stellar activity. The recomputed RVs then have a smaller scatter and the modulation due to stellar rotation decreases. This technique has been tested for other M dwarfs similar to AD Leo (spectral types 3.0 to 4.5 V, relatively low rotational velocities and high activity levels) that are well-known to exhibit strong stellar activity signals (e.g., YZ CMi, EV Lac) and appears to perform well and as expected (Lafarga priv. comm.). We found similar results regardless of the activity indicator (total RV, CRX or BIS) used to compute the correlations with the individual line RVs, with the total RV yielding slightly smaller RV scatters.

Individual spectral lines using the CARMENES data
Specifically for AD Leo, we compare the computed individual line RVs to the CRX or the BIS and estimate the strength of stellar activity based on these correlations. We focus on just the first season in 2018, namely, VIS1, since the second season only covers a bit more than one rotation cycle and this then introduces too much scatter due to the photon noise being too high. To quantify the correlation strength, we used the Pearson correlation coefficient r. We considered lines as 'inactive' with r close to 0. We also discarded lines that had time series scatter larger than 400 m s −1 , as measured from their weighted standard deviation (WSTD) RV; these are weak lines that mostly add noise to the recomputed RVs. Our results for the correlation with the CRX are shown in Fig. 6, where we generated three line subsamples: (1) all spectral lines -black, (2) where |r| ≤ 0.30orange , and (3) where |r| ≤ 0.20 -blue. We obtained similar values using the correlation with the BIS. We recomputed the RVs from both the orange and blue subsamples, thus using the least activity-affected lines. Assuming the scatter of the RV data to be dominated by the stellar activity, we then would expect a decrease in the time series RV scatter, which we do not observe. Regarding the modulation, the periodograms of these two data sets show a decrease in the power of the 2.23 d peak compared to the all-lines data set, but there is still some power left (false alarm probability (FAP) > 10 % and almost 10 % for the orange and blue subsets respectively).
This would be clear if the RV scatter decreased for the orange or blue data sets. However, since we do not observe such a decrease, it is difficult to discern what is causing the periodogram behavior. The increase in the RV scatter could be caused by increasing photon noise in the RVs of the orange and blue data sets (because we are using a smaller number of lines to compute the RVs, hence less signal). Then, the decrease in the periodogram power could be due to this increasing noise, and not to a decrease in the activity signal.
Compared to the results from other stars (e.g., EV Lac, YZ CMi), the correlations found for AD Leo are much weaker (the mean correlation coefficient r is ∼ -0.2 and some of them show |r| ≥ 0.8 for the correlation with the CRX, while for AD Leo, the mean is at 0 and very few lines have |r| ≥ 0.6). This lack of clear correlations indicates that the correlations that we find for AD Leo do not have much information related to the activity of the star, and this could be why we are not able to effectively mitigate the stellar activity signal in the RVs.
This difference in the correlation strength could be due to the different RV amplitudes of the stars. AD Leo shows a small RV amplitude compared to the other considered stars: K ∼ 25 m s −1 in comparison to EV Lac with ∼100 m s −1 . Both stars have similar spectral types and activity levels (for EV Lac, pEW(Hα) = -4.983 ± 0.021, as computed from the CARMENES observations, Schöfer et al. 2019), so the difference in RV amplitudes seems to be caused by different spot configurations. AD Leo has a relatively low inclination (i ∼ 13 • ) in comparison to EV Lac or the other considered stars (≥ 60 • , see e.g., Morin et al. 2008). This close to pole-on inclination could cause any visible corotating spots to induce a smaller modulation in the RVs simply because v sin i is smaller. Also, the photosphere of AD Leo could be more homogeneously spotted, also inducing smaller RV modulations.
To summarize, we recomputed the RVs using only the lines least affected by activity in the AD Leo spectra, and observed a decrease in the periodogram peak at 2.23 d (0.1 % to 10 %); however, there was still some significant residual power and we did not observe a significant decrease in the RV scatter. This could be due to AD Leo having a different activity signal in the RVs than other stars with similar characteristics, for which we observe a clear mitigation of activity.

Other spectroscopic activity indicators
To aid in disentangling the origin of a signal, a variety of stellar activity indicators, such as photometric variability, bisector inverse slope (BIS, the measurement of asymmetry of the CCF), Ca ii infrared triplet (IRT), Hα, CRX, the dLW, or the TiO bands, has been shown to be successful in identifying and indicating a nonplanetary signal. In the literature, periodic RV signals are often considered planetary in nature if the same signal is not found in activity indicators. Therefore, it is important to examine as many activity indicators as possible.
Correlations between the mentioned stellar activity indicators (when available) and the RVs along with the GLS periodograms for these indicators can be found in Appendix C for both the CARMENES VIS and NIR channel as well as for the HARPS instrument. The strong correlations found with the CRX was already discussed in Sect. 4.1. The BIS-RV anticorrelation for the VIS channel (r = -0.79) is inline with the same anticorrelation found by C20 in the HARPS-N data (r = -0.74). Additionally, we look at the dLW However, in our case, there seems to be no correlation. Along with the dLW, other expected stellar activity indicators such as the emission of Hα (Kürster et al. 2003;Hatzes et al. 2015 2021) for a high-mass and Hα-active star as AD Leo. In addition, we also consider the pseudo-equivalent widths (pEW') of various potential stellar activity lines as described in Schöfer et al. (2019), though find no correlations or significant peaks other than suggestive bumps at 1.8 d (i.e., the daily alias of 2.23 d).

Modeling the time dependence of the RV signal
Discontinuous sampling over the course of a long time baseline impedes obtaining key information as to what occurs in large data gaps -as is the case between the majority of the HIRES, HARPS data and CARMENES, GIARPS, HPF data for AD Leo being 12 years apart. Remarkably, the 2.23 d signal persisted over all this time, however, it is evident that amplitude changes and phaseshifts occurred. Therefore, given such a rich RV data set, we test out the current approaches for modeling stellar activity on a real-life example and present our results below.

Model setup
We model the RV data using two components: a stable component (in amplitude, period, and phase) and a variable component. The former is described by a deterministic model, typically a circular or eccentric Keplerian orbit. We would like to emphasize that when we use a Keplerian model, we consider it more as a "stable model" which may include both a planetary companion and a stable aspect of stellar activity, or only one of the two. Therefore, parameters that are typically called P planet and K planet to represent the period and semi-amplitude of the Keplerian signal are rather referred as P stable and K stable to indicate that we are considering a "stable mode" that appears persistent. From here on, we use the terms Keplerian model and stable model interchangeably, though we are not claiming that this component is solely due to a planetary companion. As for the variable component, typically comprising the quasi-periodic behavior of stellar activity, we do not currently possess a deterministic model, therefore, we use a nonparametric GP model to describe these modulations. We employ two different kernels.
The first is an exponential-squared-sine-squared kernel, that is, quasi-periodic (QP-GP) kernel, provided by george (Ambikasaran et al. 2015), where τ = |t i − t j | is temporal distance between two points, σ GP is the amplitude of the GP modulation, α GP is the inverse lengthscale 16 of the GP exponential component, P GP, rot corresponds to the recurrence timescale, and Γ GP is the smoothing parameter. The former term is an exponential that can model the decorrelation due to the changes in phase and amplitude as active regions grow and decay over time, whereas the latter term accounts for the reoccuring periodicity. The second kernel is a sum of two stochastically driven, damped harmonic oscillator (SHO) terms, or a double SHO GP (dSHO-GP), where the power spectrum of one SHO term is given by Anderson et al. (1990), SHO 1 (ω GP ) = 2 π S 0 ω 4 1 (ω 2 GP − ω 2 1 ) 2 + ω 2 1 ω 2 GP /Q 2 1 (2a) 16 The parameter α GP is defined as 1/l 2 , where l is the timescale of variations. The original relation in Espinoza et al. (2019) incorrectly defined α = 1/2l 2 . This has since been corrected. and SHO 2 (ω GP ) = 2 π S 0 ω 4 2 (ω 2 GP − ω 2 2 ) 2 + ω 2 2 ω 2 GP /Q 2 2 , for which we applied a reparametrization using the hyperparameters and where σ GP is again the amplitude of the GP kernel, P GP, rot is the primary period of the variability, Q 0 is the quality factor for the secondary oscillation, δQ is the difference between the quality factors of the first and second oscillations, and f represents the fractional amplitude of the secondary oscillation with respect to the primary one. We investigate multiple models to see which one is preferred for each individual data set as well as for the combined data set comprising all available RVs. The models being tested are: 1. Keplerian-only models, i.e, both circular and eccentric 2. GP-only model, either with a QP-GP or a dSHO-GP, as a proxy to describe the quasi-periodic nature of the stellar activity 3. Mixed (Keplerian + GP) models, to describe stable components with variable ones We use the model-fitting python package, juliet (Espinoza et al. 2019) in order to compare these models using a Bayesian framework. For our purposes, the RVs are modeled by radvel 17 (Fulton et al. 2018), and the GP models with the help of george 18 (Ambikasaran et al. 2015) (2008), we follow the rule of thumb that a ∆ ln Z greater than 5 between two models indicates strong evidence in favor for the model with the larger Bayesian log evidence (odds are ∼150 to 1), whereas a ∆ ln Z more than 2.5 indicates moderate evidence for the winning model, and anything less we consider to be inconclusive.

Prior setup
The RV periodograms (Fig. 3) indicate the interesting signal to be around 2.23 d. Thus, for all models with a Keplerian, the prior of the period of the stable component, P stable , was kept uniform, but relatively narrow, U(2.0 d, 2.5 d), in order to avoid picking up aliases (i.e., 1.8 d due to the daily sampling), and similarly for the time of transit center 21 , t 0 , U(2458200 d, 2458202 d), in order to cover solely one cycle and to avoid picking up other potential modes. The semi-amplitude prior was simply uniform as well, U(0 m s −1 , 50 m s −1 ).
As for the GP models, the rotational period, P GP, rot follows the same prior as P stable for consistency. Then, we consider wide log-uniform priors for the other QP-GP hyperparameters: Γ GP (between 10 −2 and 10 2 ) and α GP (between 10 −10 d −2 and 10 −1 d −2 , corresponding to timescales of ∼3 d -∼70 000 d). For the dSHO-GP hyperparameters: Q 0 was log-uniform (from 10 2 to 10 5 ), as well as δQ (from 10 −1 to 10 5 ), and f was kept uniform (between 0 and 1). For both GP kernels, the sigma of the GP was kept log-uniform (between 0.1 and 70 m s −1 ) where each instrument had its own, noted as σ GP, inst , because each instrument has its own characteristics, such as noise level, zero point offset, wavelength range, or intrinsic stellar jitter. The other GP hyperparameters are shared. Additional instrumental jitter terms (log-uniform from 10 −2 to 30 m s −1 ) and offsets for each individual instrument were considered as well.

Results
We applied this recipe to the three following data sets: CARMENES VIS only, CARMENES NIR only, and the entire RV data set. The resulting posteriors on the P stable (P GP, rot if using a GP), semi-amplitude K stable (σ GP, inst if GP), added instrumental jitter terms (σ inst ) and eccentricity along with the differences in Bayesian log evidence can be found in Table 3.

CARMENES only
Focusing on just the CARMENES season one data, an eccentric (e ∼ 0.19) Keplerian model is preferred over a circular model for the VIS channel (∆ ln Z ∼ 5). This is in agreement with the optical HARPS-N data where C20 also found a similar eccentricity when applying a Keplerian-only model. As for the NIR channel, a circular model performs best. Moreover, the NIR data does not actually have a clear model preference, simply attributed to the lower precision. The amplitude decreased from 28.43 +0.94 −0.98 m s −1 in the VIS1 data to 18.0 +2.6 −2.7 m s −1 in the NIR1 data. When introducing combined stable plus GP models, the stable component dominates and the GP component is not needed, in other words, σ GP becomes consistent with zero. This finding is anticipated given that the time span of the data is roughly ∼10 cycles of the periodicity, which would be too short for any noticeable changes of the stellar spot pattern assuming a long spot lifetime.
For the season two data, we slightly readjusted the priors as introduced in Sect. 6.2. The time of transit center was moved to U(2458893 d, 2458895 d) to comply with the pertinent time stamps. Likewise, the period was narrowly constrained, U(2.22 d, 2.24 d), to ensure the correct periodicity is recognized given the short time baseline of less than two periodic cycles. Likewise for this reason, just the circular Keplerian-only models were performed, and as a result, the posterior values are not included in Table 3. Similarly to the season one data, there was also an amplitude decrease between the instruments, namely from 17.02 +0.63 −0.62 m s −1 to 9.2 +3.0 −2.6 m s −1 between VIS2 and NIR2, respectively, both of which are lower with respect to their season one counterparts. Figure 4 shows how the time series and phasefolded plots using the best model fits for the four CARMENES data sets appear reasonable, showing a uniformly distributed scatter in the model residuals. We additionally tested out combining the two seasons with each respective instrument (i.e., VIS1+VIS2, NIR1+NIR2) and applied all the models (including those with GPs), finding that the circular Keplerian was favored for both.

Whole data set
Next we considered the whole AD Leo RV data set. We expect the Keplerian-only models, where the amplitude for all given instruments is shared, to perform poorly given that the amplitude is clearly decreasing as a function of wavelength (shown in Sect. 4). The idea of combining all available data sets covering a wide wavelength range is that the data set with the smallest amplitude, specifically the NIR data, will act as an upper limit for any stable signal, whereas the GP component will adapt the difference. Within the framework of juliet, the data is first attempted to fit the deterministic model as best as possible, where the GP is then applied to the residuals. Likewise, we expect that the added jitter terms and GP amplitude parameters for the redder instruments will be consistent with zero in the mixed model because the stable model will most likely saturate for the red-der instruments. The role of the GP is to account for the excess amplitude observed with the bluer instruments.
The preferred model among all is the mixed circular Keplerian + dSHO-GP model. When compared to a Keplerianonly model, it exceeds in Bayesian log evidence tremendously (∆ ln Z ∼ 180), but it does not prevail against a dSHO-GPonly model that enormously (∆ ln Z ∼ 3.5). As anticipated, the Keplerian-only models fail to explain the data in the sense that finding a common amplitude is not feasible such that the jitter values become too high, and likewise, the phase shifts are quite strong that a simple sinusoid can not describe this behavior. Surprisingly, even though the dSHO-GP-only model is flexible and could model the data well, the evidence suggests that a model containing an additional stable periodic component is favored. Likewise, the same concept applies for the models including the QP-GP kernel, though including the additional stable component produced a ∆ ln Z ∼ 10 rather than ∼ 3.5 for the dSHO-GP models, that is, very strong evidence for one model versus bordering just moderate evidence. This delicate boundary could change the interpretation. Tests on simulated RV data based on StarSim (Herrero et al. 2016) indicate that a Keplerian signal is in rare cases more efficient than the QP-GP in modeling a coherent activity signal (Stock et al. 2022).
We find that the stable signal amplitude K stable is 16.6 ± 2.2 m s −1 , where both the P stable and P GP is rounded to 2.23 d. All posteriors for the different models can be found in Table 3 and the phase-folded plots for the best model, namely the mixed model, are shown in Fig. 7.
To summarize our findings, we have the following takeaway points: -The CARMENES-VIS1 data alone prefer an eccentric Keplerian (e ∼ 0.19) model. This is in agreement with the HARPS-N data from C20. Whereas, the CARMENES-NIR1 data alone prefer a circular Keplerian, that is, a sinusoid.

Discussion and future outlook
In this analysis, we presented new CARMENES optical and near-IR data that cover a wide wavelength range well into the red part of the spectrum to test the origin of the 2.23 d signal found in the RVs for AD Leo. The presence of stellar activity is supported with CARMENES through the proof of wavelength-and time-dependence. With the CARMENES data alone, the shape of the wavelength dependence of the RV amplitude can be attributed to a star-spot configuration following Reiners et al. (2010). Between the two CARMENES seasons, the amplitude decreased in both instruments. Specifically, from 28.43 +0.94 −0.98 m s −1 to 17.02 +0.63 −0.62 m s −1 in the VIS channel, and from 18.0 +2.6 −2.7 m s −1 to 9.2 +3.0 −2.6 m s −1 in the near-IR channel when modeling the signals as pure sinsuoidal variations. C20 also found an amplitude decrease between the HARPS-N and GIANO-B data, that is, from 33 m s −1 to less than 23 m s −1 between the instruments and from 33 m s −1 to 13 m s −1 between seasons. However, given the large errorbars of the GIANO-B data (∼20 m s −1 ), it was only possible to identify an amplitude decrease, but not to constrain whether there even is a signal present in these data at 2.23 d. Likewise, the high-precision near-IR data from HPF show a much smaller RV amplitude (R20), inconsistent with the higher RV amplitude in the optical, and there is also a discrepancy even between HPF seasons, namely, the RMS scatter dropped from 23 m s −1 to 6.4 m s −1 . This decrease in amplitude is also seen in the overlapping, optical HARPS-N data and second season of CARMENES data (see again Fig. 2). Nonetheless, both our results and those in C20 and R20 agree that the signal's amplitude does indeed decrease rather than increase, indicating that AD Leo is entering a lower-activity phase. We conclude that the effect of the spots on the RVs is dominated by temperature differences rather than the Zeeman effect (see Reiners et al. 2013), in accordance with the other evidence at hand such as photometric variability (presented in T18).

Possibility of stable stellar-spot signal
It is common practice to model activity-induced RVs with GPs (e.g., Rajpaul et al. 2015, and further citations), as they are nondeterministic models that can sufficiently fit stochastic modulations. Provided that the evidence against a planet outweighs that in favor of a planet, our understanding was to expect that the GPonly model should have been adequate enough in accounting for the stellar activity, with the assumption that the signal in AD Leo is purely stellar activity induced. The modeling instead showed that neither a stable-only (i.e., sinusoid) model nor red-noiseonly (i.e., GP-only) model can correctly describe the whole data set (Table 3). A mix of a stable and variable component is rather the preferred model, where the stable model has an amplitude of K stable = 16.6 ± 2.2 m s −1 , lower than the value of 19 m s −1 proposed by T18. We do not claim that the stable component is solely due to a planetary companion as it could perhaps be also due to the persistent presence of a circumpolar stellar spot imposing a constant behavior, or even a mixture of both. Thus, we set an 3σ upper limit (i.e., from K stable = 16.6 ± 3 × 2.2 m s −1 ) of 27 M ⊕ (or 0.084 M Jup ) for a putative planetary companion, in comparison to the proposed mass of ∼0.2 M Jup by T18. Focusing on the amplitudes of the GP kernel of the mixed model, we found that those of the earlier optical RV data, specifically HARPS and HIRES (σ GP, inst ∼ 20 m s −1 ), were higher than most of the more recent optical data (σ GP, inst ≲ 13 m s −1 ). These values support the assumption that the star was more active earlier and is entering a less active phase (Sect. 2.2). The only exception was with HT2, with σ GP, HT2 ∼ 22 m s −1 , an even higher value. The GP amplitudes of near-IR spectrographs were consistent with σ GP, inst ∼ 0 m s −1 , except for HPF2 (σ GP, HPF2 ∼ 8 m s −1 ).  Phase-folded plots for the optical (left) and near-infrared (right) instruments using the circular Keplerian + dSHO-GP model after subtracting the GP component out. The first subset of any given instrument is represented with a circle, whereas the second subset by a triangle, when applicable. The y-axis ranges are consistent between the two plots, except for the residual plots to better visualize the scatter around the fit.
The data from the two instruments that stand out (i.e., HT2 and HPF2) were contemporaneously taken, indicating that there was some phase shift between the stable and variable component during this time period. A better understanding of how model comparison can behave given a purely stellar-activity induced signal would be necessary to determine if such a result is even expected (see also Sect. 7.2 for future suggestions). Additionally, our analysis of recomputing the RVs, considering only those spectral lines unaffected by stellar activity, still showed some strong significant periodicity at 2.23 d (Fig. 6), unlike for other known active stars where the signals had disappeared (Lafarga priv. comm.).
Nonetheless, this should raise the question of whether a signal with an amplitude of ∼17 m s −1 that is fully stable over all these observations (∼19 yr) in such an M dwarf is possible. In fact, it is not so surprising that spot-induced RV fluctuations for M dwarfs are long-lived (e.g., Günther et al. 2020;Quirrenbach et al. 2022). Though evidence of stable stellar activity behavior had previously been shown in photometry (e.g., GJ 1243, Davenport et al. 2020) and RVs (e.g., α Tau, Hatzes et al. 2015) over time, our paper demonstrates the first case for modeling RVs over time and wavelength. Future studies performing simulations with software packages such as StarSim 2.0 (Herrero et al. 2016) or SOAP 2.0 (Dumusque et al. 2014) using various star-spot configurations may also shed light on why a signal can stay persistent over many years (Herrero et al. 2016;Rosich et al. 2020). Utilizing StarSim 2.0 to compare RV-CRX correlations between simulated and real data as performed by Baroch et al. (2020), with YZ CMi as a case study, could be beneficial in determining the star-spot temperature difference, the star-spot filling factor, and the location of the spot. A first test with a simple assumption (i.e., one big polar spot) could reproduce well the RVs and CRX values of AD Leo for the CARMENES VIS data (Baroch priv. comm.). But such an approach can become degenerate when considering so many instruments and various star-spot configurations, and a planetary signal would act as an achromatic offset to the RVs and CRX.

Current limitations
AD Leo is a prime case study for M dwarfs to explore the possibility of a planetary signal with an orbital period indistinguishably close to its stellar host's rotation period. Given its rich multi-wavelength spectroscopic monitoring over a large time baseline, there are no other targets with such an extensive data set coverage. Even then, we showed that there are apparent holes in our current state of research in this field in which studies to address them would require proposals, dedicated surveys and telescope time. In essence, we are either limited by our modeling approaches, by our astrophysical knowledge of star-spot configurations, or a combination of both. In order to progress and disentangle this problem, we propose the following, in no particular order.

Simultaneous multiband photometry
The shape of the curve in Fig. 5 is determined by stellar activity behavior, whereas the presence of a planet would simply act as an offset across the wavelength space. Even when doubling (or halving) the number of available spectroscopic data points, disentangling the contribution of a planetary signal within this curve still persists as a degenerate problem. Obtaining continuous and simultaneous multiband photometry would help in painting a better picture of the stellar-spot map distribution. In doing so, this can then be translated by forward modeling to determine the stellar-activity induced RVs (e.g., via Starsim 2.0 or SOAP 2.0; Herrero et al. 2016;Dumusque et al. 2014). The residual between these computed RVs and those obtained would indicate the contribution of the Keplerian component. However, the caveat is, that spot configurations on the stellar surface can vary with time. Thus, continuous photometry would only apply to simultaneously taken RVs. Table 3. Model comparison of RV fits done with juliet comparing a Keplerian model, a red-noise model, and a mixture of both for the various instrument data sets: CARMENES VIS1, CARMENES NIR1, and the combined data set. The selected best model is boldfaced (refer to Sect. 6 for details on the priors used, model choice, and discussion of results). The columns refer to which model was being used, whereas the rows correspond to the model parameter. A larger, positive ∆ ln Z value indicates a better model.

A better-suited RV data set for spectral line analyses
In the absence of available photometric facilities, the approach of selecting spectral lines not affected by activity seems promising. This line-by-line approach has been successfully applied to other stars (e.g., Dumusque 2018; Wise et al. 2018;Cretignier et al. 2020). In the case for AD Leo, it is difficult to decipher whether it was the target itself or the relatively new method that resulted in an inconclusive outcome. Perhaps the RVs were not as affected by the stellar spots (i.e., due to the pole-on inclination), the RV amplitude was not as large as the others, or there were not enough data points to strengthen the correlation plots by establishing a statistically stronger correlation. Definitely, if over-all activity decreases, then it would be more difficult to pick up the activity within certain spectral lines. Hence, observing a star when it is in its most active phase would constitute the ideal data set. Additionally, our spectral line analysis should be applied to a wider wavelength range (e.g., HARPS) as well to determine if the signal disappears. Another possible explanation could be that all lines are inherently affected by activity to some degree. Nonetheless, determining the ideal targets for this technique as well as developing the method further will be paramount in suppressing activity for future M-dwarf studies.

Magnetic star-planet interaction
The radio bursts thus far observed for the AD Leo system (outlined in Sect. 2.5) suggest an excess that is unaccounted for. A possible scenario to justify the discrepancy could be due to sub-Alfvénic interaction of the stellar magnetosphere with an orbiting body, in this case with a putative planet around AD Leo, which has P orb = P rot . This magnetic interaction between the star and a planet (Magnetic Star-Planet Interaction; MSPI) would be responsible for producing: (1) an enhancement of the magnetic activity of the star, and (2) the long-duration radio emission as that detailed in Sect. 2.5, both effects of which are modulated in phase with the orbital period of the planet. In a star like AD Leo, the difficulties for detecting this modulation are that this target is seen nearly pole-on (i ∼ 13 deg) and that a putative planet orbit is synchronized with the stellar rotation. This produces a weaker MSPI, as it depends on the relative velocity of the planet with respect to the magnetosphere of the star (Lanza 2009(Lanza , 2012 and makes it difficult to disentangle it from stellar intrinsic activity. However, the presence of such a planet could explain the need for the stable component and part of its semi-amplitude (K stable < 16.6 m s −1 ) seen in the RVs in this work. MSPI would show as a hot spot in the chromosphere that would mainly affect the chromospheric activity indicators such as log R ′ HK from the Ca ii H&K lines, not available in the CARMENES data set but covered by HARPS.
Within this scenario, we followed the prescriptions in Appendix B of Pérez-Torres et al. (2021) to estimate the flux density expected to arise from the interaction between a putative planet and its host star. We assume for simplicity an isothermal Parker wind with a temperature of (2-3)×10 6 K, and a density of 10 7−−8 cm −3 at the base of the corona. Since the hypothetical planet would be very close to its host star, we find that the planet is in the sub-Alfvénic regime, which makes it possible that energy and momentum can be carried upstream of the flow along Alfvén wings. We estimate that the radio flux in the few GHz band, arising from star-planet interaction, is in the range from less than about 0.1 mJy and up to a few mJy, or even higher, depending on the stellar wind parameters at both the corona and at the orbital position of the planet. We therefore encourage a campaign for simultaneous monitoring of AD Leo in both optical (i.e., spectroscopy, including Ca ii H&K lines, and photometry) and in radio wavelengths to test the star-planet (sub-Alfvénic) interaction scenario. The detection of a correlation in the radio signal with the orbital period would be the smoking gun of sub-Alfvénic interaction, serving as an argument in support of the presence of a planet in the system.

Application of GP models for red-noise correlated signals
Additional data may not be what is necessary as much as is a better understanding of the astrophysical effects stellar spots impose on the star and how we choose to model them through nonphysically motivated means. The usage of GPs is currently the most popular approach in the RV field for modeling stellar activity. That being said, the most commonly used kernels (i.e., QP-GP, dSHO-GP) are not completely physically motivated, but nonetheless serve as a very good approximation to our aim of describing the quasi-periodic behavior. They can play a vital role, especially when trying to search for low-mass planets with low amplitudes among noisy data (e.g., Demangeon et al. 2021;Faria et al. 2022).
The problem is, however, that we do not know just how effective the GPs are in serving their purpose and what limitations and consequences on (planetary) parameters they particularly have. Perhaps, the choice of the GP for different time baselines could be also of importance, as is the case for AD Leo (see Sect. 6), considering there are, for the time being, too many unknowns that can lead to misinterpretations. Newer implementations of GP kernels with a potentially better physical interpretation of stellar processes could become essential (e.g., Perger et al. 2021;Luger et al. 2021). Applying these models to real data sets turns out to not be the solution, as it becomes complicated to correctly interpret the results. For this reason, these models should be compared against each other on controlled, simulated RV time series. A more in-depth study could include simulating RV time series by injecting solely stellar-activity-induced signals and a mixture of stellar activity and planetary components for a general grasp (Stock et al. 2022). One step further would be to specifically simulate the stellar characteristics of AD Leo, together with the time stamps and various wavelengths. Even here, the limitation of our knowledge of spot distribution maps and correctly choosing the correct one persists.
In future potential cases where there is a presumed planetary and stellar activity signal with two instruments covering a wide wavelength range, for instance, with CARMENES, it is crucial to adequately model the wavelength dependence of the stellar activity to ensure the most precise planetary parameters. More so for M-dwarf stars considering that their stellar rotation periods would often correlate to the orbital periods of planets residing in the habitable zone (Newton et al. 2016). AD Leo in this regard is an interesting case study to perform such fits when considering real data, however, simulated data sets may offer a more controlled environment to test the effectiveness of our current modeling techniques. Either way, AD Leo certainly serves as a particularly intriguing system for studying the impact of a stable spot periodicity on the search for planetary signals.

Using auxiliary data to correct for stellar activity
Another promising procedure is to include the stellar activity indicators as auxiliary data to try to detrend the RVs, commonly done so as linear components (e.g., Anglada-Escudé et al. 2016) or within the GP kernel (e.g., Suárez Mascareño et al. 2020;Cale et al. 2021;Barragán et al. 2022) in regards to planet detection. Likewise, this technique is being carried out in a sole stellar activity context (Jeffers et al. 2022;Cardona Guillén et al. 2022), with particular attention to not only detrending with various indicators, but also to characterizing the behavior of the correlations. Specifically, this "closed-loop" relation is becoming more prevalent, as already exhibited with known active M-dwarf stars such as YZ CMi Baroch et al. 2020), EV Lac (Jeffers et al. 2022), and GJ 674 (Bonfils et al. 2007), and now for AD Leo as displayed in the RV-CRX and RV-BIS correlations (Fig. C.1). Subtracting out the CRX dependency can improve the rms of the RVs, though it is nearly impossible to obtain a straight line from the correction, that is, completely removing all of the stellar activity, even in cases with very high signal-to-noise ratios and very dense sampling (e.g., factor of three improvement for EV Lac; Jeffers et al. 2022). For our case of AD Leo (Sect. 4.1), the CRX detrending improved the rms roughly by a factor of two for VIS1 (from 18.5 m s −1 to 9.2 m s −1 ).
In this line, a relatively recent technique is that of using spectropolarimetric measurements to correct the RVs of the variability introduced by the activity. The data are used to derived precise RVs and to obtain maps of the magnetic field and the dis-tribution of surface inhomogeneities. These maps are then used to obtain a model RV curve which are the subtracted from the observed RVs. The removal of activity variations with this technique decreased the amplitude of the RV curve in V830 Tau by a factor 10, which allowed the detection of a hot Jupiter orbiting this very young T Tauri star (Donati et al. 2016).
While these approaches are encouraging, they still poses limitations when considering a putative planetary signal whose orbital period is equal to the rotational period because the same issue, as mentioned already above, of degeneracy still holds true. Again, a systematic simulated analysis of forward modeling RVs and stellar activity indicators could be key in understanding what correlations to expect and how to properly detrend them.

Conclusions
In this paper we presented new CARMENES optical and nearinfrared spectroscopic data for a known nearby, active M dwarf, AD Leo. The stellar rotation period of 2.23 d is clearly present in the RVs, and we address the question of whether there could be a planet orbiting with the same periodicity, P orb = P rot . Taking advantage of the wide wavelength range of the CARMENES instrument, we demonstrated the shape of the wavelength dependency of the RV semi-amplitudes to be in agreement with what is expected for a star-spot temperature difference configuration. The strong anticorrelation found between CRX-RV and BIS-RV for the CARMENES VIS data were additional signs of stellar activity. In addition, we recomputed RVs using spectral lines unaffected by stellar activity, which has not been done before for AD Leo, and found that there was still some significant residual power at the rotational period.
When incorporating all available RVs (HARPS, HIRES, CARMENES, and the more recent HARPS-N, GIANO-B, and HPF data), it became evident that the signal has undergone various amplitude fluctuations as well as phase shifts, and this is behavior that cannot be attributed to the presence of a planetary companion. However, a closer look into the model comparison showed that a mixed model of a stable plus a quasi-periodic red-noise model best explained the data, where the stable component had a semi-amplitude of K stable = 16.6 ± 2.2 m s −1 , setting a 3σ upper limit on the mass of a putative planet, 27 M ⊕ (= 0.084 M Jup ). Based on all this evidence, our current machinery, and the given data at hand, we conclude we cannot unequivocally prove, or even disprove, the 2.23 d periodic signal found in AD Leo to be solely due to stellar activity. We therefore suggest to obtain simultaneous photometry and spectroscopy in various wavelengths accompanied by radio observations, and complemented with a clearer insight of stellar activity behavior on measurements via simulated time series to break the degeneracy.
A&A proofs: manuscript no. output  BJD (TDB * ) RV (m s −1 ) σ RV (m s −1 ) Table A.4. HARPS RV and accompanying data for AD Leo used in this paper, first processed using serval ) and then corrected for nightly zero points (Trifonov et al. 2020