© The Authors 2024

1 Introduction

Understanding the processes that govern planet formation and evolution is one of the key questions in planetary science. To address this, several surveys aiming at detecting planets through photometric transits (e.g. Kepler and TESS collaborations) have been operational. Most planet discoveries thus far have predominantly centred around main sequence stars (Pepe et al. 2004), with a smaller but still significant proportion found around post-main sequence giants (Veras 2016), brown dwarfs (Chauvin et al. 2004), intermediate age stars (Quinn et al. 2014), and pulsars (Wolszczan & Frail 1992). Interestingly, only a few dozen planets have been discovered around younger stellar systems (≲ 10 Myr), most of them being through direct imaging (e.g. Lafrenière et al. 2008; Hammond et al. 2023; Wagner et al. 2023), a few being through the transit method (e.g. Rowe et al. 2014; Mann et al. 2016), and only one by the radial velocity method (V830 Tau, Donati et al. 2016). Another method leading to planet detection around young T Tauri stars (TTSs) was made possible by imaging gaps in their disks (e.g. Keppler et al. 2018; Pinte et al. 2019; Wang et al. 2020). This method is limited to detecting them at relatively large distances from their host stars. This limitation underscores the current constraints of the latest generation telescopes, which have yet to achieve the necessary high resolutions to detect planets within the innermost active regions (≲0.1 au). Consequently, it remains uncertain when and where close-in planets form; also uncertain is the pace at which they migrate either inwards or outwards through the disk.

The relative paucity of confirmed planetary companions to very young stars highlights the difficulty in detecting them around such stars, particularly in the T Tauri class, which are among the most variable of the young stars. TTSs are divided into two classes: weak line TTSs (WTTSs) and classical TTSs (CTTSs). WTTSs are systems that lack spectroscopic accretion signatures and generally display rather stable sinusoidal light curves induced by surface features (e.g. Stassun & Wood 1999; Grankin et al. 2008; Frasca et al. 2009; Donati et al. 2014). Detecting planets using the radial velocity (RV) method in such stable systems is less challenging than in CTTSs, which display complex characteristics in their light curves (and RVs) as a result of intense stellar activity and/or active accretion phenomena (e.g. Vrba et al. 1993; Herbst et al. 1994; Rucinski et al. 2008; Siwak et al. 2010; Cody & Hillenbrand 2010; Alencar et al. 2010; Fischer et al. 2023).

CI Tau remains the only CTTS in which a disk-embedded close-in planet with ~ 11 M_Jup and a period of ~ 9 d was claimed to be detected by Johns-Krull et al. (2016) and later by Biddle et al. (2018). However, such claims were challenged following Zeeman-Doppler Imaging (ZDI) analysis of ESPaDOnS spectra, which demonstrated that the 9 d modulation corresponds to the stellar rotation period and is induced by stellar activity rather than a purported planet (Donati et al. 2020a).

CI Tau is estimated to be approximately 2 Myr old (Guilloteau et al. 2014). Simon et al. (2019) estimated a mass of 0.9 ± 0.02 M_⊙ based on pre-main sequence (PMS) evolutionary tracks, while an earlier study by Guilloteau et al. (2014) reported a slightly lower mass of 0.8 ± 0.02 M_⊙, obtained using rotation of its Keplerian disk. The central star has an effective temperature of 4250 ± 50 K (Donati et al. 2020a) and is positioned at 160.3 ± 0.4 pc from the Sun in the Taurus molecular cloud (Gaia Collaboration 2023). CI Tau is recognised for hosting a robust predominantly radial magnetic field of up to 3.7 kG and displays a varying mass accretion rate of around 2 × 10⁻⁸ M_⊙ yr⁻¹ (Donati et al. 2020a). On larger scales, the central star is surrounded by a circumstellar disk observable in millimetre continuum images, extending up to 200 au and inclined at ~55 ± 5°, as reported by Guilloteau et al. (2011) and Clarke et al. (2018). This disk exhibits a sequence of dusty rings, with discernible gaps at radii of approximately 13, 39, and 100 au, suggesting ongoing planet formation processes (Clarke et al. 2018). On smaller scales, CI Tau is surrounded by an inner disk with an inclination of ~70^° (GRAVITY Collaboration 2023), mis-aligned with the outer disk.

In this paper, we report the detection of a significant periodic signal at ~25.2 d in both photometric (K2 and LCOGT) and spectroscopic time series (ESPaDOnS and SPIRou RVs). The periodic signal in the RVs becomes clearer when stellar activity due to spot modulation is mitigated. In Sect. 2 we describe the datasets (photometric and spectroscopic) used in this study. Sections 3 and 4 respectively detail our photometric and spectroscopic time series data analysis. We discuss our results in Sect. 5, and conclude in Sect. 6.

2 Observations

2.1 K2 photometry

We used publicly available data from the K2 mission (campaign 13, Howell et al. 2014). Long-cadence simple aperture photometry (SAP) flux light curves were downloaded using the LIGHTKURVE package (Lightkurve Collaboration 2018) to convert the raw data into a FITS file, from which the usable data were extracted.

2.2 LCOGT multiband photometry

Additional photometric time series data for CI Tau were obtained at Las Cumbres Observatory Global Network (LCOGT, Brown et al. 2013) within two epochs separated by two years apart. The first epoch (20B) covered from October 30, 2020 to January 31, 2021 (JD 2459 153–JD 2459245) where 539 images were acquired in the Sloan 𝑔′r′i′ filters over three months with a sub-day sampling rate (runs LCO2020B-019, PI: L. Rebull; CLN2020B-005, PI: A. Bayo, for a total of 11 h). The second epoch (22B/23A) extended from September 27, 2022, to March 31, 2023 (JD 2459 850–JD 2460033), during which 564 images were acquired over six months in the Sloan 𝑔′r′i′ filters, with a sampling of typically one measurement per day (runs DDT2022B-003, DDT2023A-001; PI: L. Rebull, for a total of 11 hours). At both epochs, the 𝑔′r′i′ images were obtained with the 0.4 m SBIG telescopes offering a field of view of 29.2 × 19.5 arcmin, with exposure times of 60, 30, and 30 s, respectively. We retrieved the BANZAI reduced images from the LCOGT archive service and the non-calibrated photometric catalogs provided in the image headers for all detected stars in the field.

In order to compute differential photometry, we considered three stars, 2MASS J04334286+2255065 (JH104), 2MASS J04334414+2256179 (JH 105), and 2MASS J04340717+2251227, located within 6 arcmin of CI Tau and of similar brightness. We used JH 105 as a comparison star to compute the differential light curve (CI Tau-JH 105), while 2MASS J04340717+2251227 served as a check star to assess these are non-variable sources. The differential light curve (CI Tau–JH 105) was then calibrated in each filter by adopting JH 105’s magnitudes as provided by the APASS survey, namely, 𝑔′ = 13.235, r′ = 12.238, i′ = 11.666mag. The measurement uncertainty was calculated as RMS/sqrt(2) of the differential light curves of the non-variable comparison and check stars. For both epochs, the measurement uncertainty amounts to 0.020, 0.024, and 0.030 mag, in the 𝑔′r′i′ filters, respectively. Some measurements of lower quality, due to non-photometric conditions and/or the proximity of the bright Moon were discarded. Specifically, we removed images where the comparison and check stars’ differential photometry indicated a measurement discrepant by more than 2σ from the mean. Eventually, for the first epoch, we retained 156 measurements in the 𝑔′-band, and 164 in each of the r′ and i′ -bands, thus sampling the light curves at a sub-day rate nearly continuously over 93 days during semester 20B. For the second epoch, we gathered 165 measurements in the 𝑔′ -band, 172 in r′ -band, and 175 in i′ -band, thus sampling CI Tau’s photometric variability nearly continuously at a daily rate over 183 days during semesters 22B/23A.

2.3 ESPaDOnS data

Spectroscopic observations were conducted using the ESPaDOnS spectropolarimeter at the 3.6 m Canada-France-Hawaii Telescope (CFHT) on Mauna Kea, Hawaii (PI: J. F. Donati). A total of 72 spectra were obtained over a time span of approximately two months, from mid-December 2016 to mid-February 2017. These observations were divided into 18 sequences, with each sequence comprising four consecutive nightly sub-exposures. Each spectrum covers a wavelength range of 3700–10 000 Å at a resolving power of 65 000 (Donati 2003). The data were reduced using the standard ESPaDOnS reduction package (Libre-ESpRIT), which is a newer version of the ESpRIT pipeline (Donati et al. 1997).

A least-squares deconvolution (LSD) was used to combine a large number of photospheric spectral lines into a mean LSD profile with an improved S/N (Donati et al. 1997; Kochukhov et al. 2010). The line mask used for calculating the LSD profiles was extracted from the VALD database (Piskunov et al. 1995; Kupka et al. 1999) for wavelengths between ~3900 Å and 10 500 Å using a template mask derived from a synthetic spectrum corresponding to the effective temperature and gravity of CI Tau (Guilloteau et al. 2014).

2.4 SPIRou data

SPIRou is a fibre-fed, cross-dispersed échelle spectropolarimeter mounted on the Canada-France-Hawaii telescope (CFHT) on Maunakea, Hawaii. The spectrograph provides spectral coverage over a wavelength range from 950 to 2350 nm at a spectral resolving power of 75 000 (Donati et al. 2020b).

In total, 425 spectra (four sub-exposures per night) were obtained for CI Tau during five semesters ranging from 18 December 2018 to 13 January 2023 (PI: Jean-François Donati), covering a total time span of ~1487 days, but with significant gaps in the RV time series and a varying number of RV points among semesters. We measured the mean signal-to-noise ratio (S/N) in three regions covering wavelengths centred at ~1209 nm, ~1557 nm and ~2329 nm, respectively, to be 40, 76 and 115. We note an overall improvement in S/N with MJD, possibly a result of a combination of the instrument optimisation after MJD ~ 58 900 (optical transmission and guiding precision, C. Moutou, priv. comm.) and an overall increase of the total integration time for the four polarimetric sub-exposures per night, ranging from ~1580 s for the earlier observations (around 2018-2020) to ~2000 s for the newer dates (around 2021–2023). Nevertheless, the spectra taken at earlier dates are also valuable.

2.5 ESO ExTrA: Low-resolution near-infrared spectroscopy

We obtained additional data as part of the Exoplanets in Transits and their Atmospheres survey (ExTrA, Bonfils et al. 2015). The ExTrA facility, situated at La Silla Observatory in Chile, comprises three 60 cm telescopes and a single near-infrared (0.88–1.55 µm) fibre-fed spectrograph.

CI Tau was observed on 136 nights between September 28, 2022 and March 21, 2023, using one telescope. Fibre units are located at the focal plane of each telescope, each consisting of two 8″ aperture fibers. One fibre is used to observe a star and the other is used to observe the nearby sky background. We used the higher-resolution mode of the spectrograph (R ~ 200) and 300-s exposures. We obtained between two and 51 exposures per night for a total of 1807 spectra with a median signal-to-noise ratio of 85 at 1.05 µm. The ExTrA data were corrected for dark current, extracted using the flat-field, corrected for sky background emission, and were wavelength calibrated using custom data reduction software. Median spectra of CI Tau were computed for each night and telescope, yielding a total of 136 spectra with a median signal-to-noise ratio of 129 and a standard deviation of 40.

3 Photometric analysis

3.1 K2 light curve

The light curve of CI Tau (Fig. 1: top panel) displays short-term erratic, aperiodic variations on timescales of hours to a few days. These are likely induced by surface inhomogeneities resulting from the combination of cool magnetic and hot accretion spots that appear and disappear on such timescales (Herbst et al. 1994; Grankin et al. 2007; Rucinski et al. 2008; Findeisen et al. 2013). Such variability is commonly identified as red noise in the data and will be explored further in Sect. 3.3. We computed a generalised Lomb-Scargle periodogram (GLS, Zechmeister & Kürster 2009), as implemented in ASTROPY (Astropy Collaboration 2018), to search for periods in the range ~0.1–100 d, using a frequency grid sufficiently dense to provide 50 samples across a given periodogram peak. The GLS periodogram is shown only for periods between 0.1 and 50 d. We also include an L1 periodogram (Hara et al. 2017) of the K2 light curve in Fig. 1, showing the most significant periods with their associated false alarm probabilities (FAPs).

An L1 periodogram is a variant of the traditional Lomb-Scargle periodogram that uses the L1 norm as a measure of goodness-of-fit criterion. The key difference between the traditional Lomb-Scargle periodogram and the L1 periodogram is in the optimisation criteria they use. While the traditional Lomb-Scargle periodogram uses the least squares method to minimise the sum of squared differences between the data and the fitted sinusoidal model, the L1 periodogram analyses multiple sinusoids simultaneously. Furthermore, the L1 periodogram proves particularly useful as a supplement to the Lomb-Scargle peri-odogram, offering a more lucid representation of the number of peaks and their significance. If for example the peaks identified by the L1-periodogram result in a chi-square value of the residuals consistent with periodogram noise, it suggests there may not be many additional signals and therefore helps avoid selecting spurious peaks.

The combined analysis using the GLS (and L1 periodogram), yields periodic variability with significant peaks at approximately 6.6 ± 0.2 d (6.59 d), 9.0 ± 0.5 d (9.24 d), 11.5 ± 0.5 d, 14.2 ± 0.8 d, and 24.2 ± 5 d (24.3 d), as well as (47.2 d). We compute log₁₀ FAPS of −270, −196, −46 and −7 for the 24.3 d, 6.59 d, 9.24 d and 47.2 d periods found in the L1 periodogram, respectively.

In Fig. 2, we present the K2 light curve phased on the dominant period peaks at 6.6d, 9.0d, and 24.2d. Errors on the periods listed for the GLS periodogram were obtained by considering the full-width at half-maximum (FWHM) of the individual periodogram peaks assuming they are Gaussians.

Biddle et al. (2021) investigated the possibility of the 24.2 d period being a synodic period of any combinations of the 6.6 d, 9.0 d and the 14.2 d periods. The authors combined sinusoidal functions with periods of 6.6 and 9.0 days, and 14.2 and 9.0 days respectively, using the same sampling as the K2 data. None of the combinations were able to generate a peak comparable to the one observed in CI Tau with a period of ~24.2 d and therefore, it is unlikely that the 24.2 d peak corresponds to a synodic period. We revisit the beat frequency (ƒ_beat) calculations from combinations of any two significant frequency peaks (ƒ₁ and ƒ₂) that we find in the GLS periodogram, using: ƒ_beat = |ƒ₂ − ƒ₁|. We find that the 14.2 d period is likely a beat period of the 24.2 d and the 9.0 d periods. While the 11.5 d period could be due to beating effects between the 6.6 and 14.2 d peaks.

The nature of the 6.6 and the 9.0 d periods have been debated in the literature. Johns-Krull et al. (2016) and later Biddle et al. (2021) claim that the 6.6 d period corresponds to the stellar rotation period and the 9.0 d is likely induced by a ~11 M_Jup planet. This view was later refuted by Donati et al. (2020a, 2024), where ZDI analysis of ESPaDOnS spectropolarimetric data demonstrate clearly that the ~9.0d period is due to spot rotational modulation. In this paper, we focus on the nature of the longerterm variability, mainly the highest peak seen in the K2 light curve at 24.2 d. Henceforth in the text, for convenience, we refer to this period identified in the K2 light curve as 24 d.

Fig. 1 Light curve, GLS and L1 periodograms for the K2 data. Top: K2 light curve of CI Tau. Middle: GLS periodogram of the K2 light curve, showing three peaks at 6.6, 9.0, and 24.2 d indicated with blue dashed lines. Bottom: L1 periodogram, showing the same periods as found in the GLS periodogram, with an additional significant peak at 47.2 d.

Fig. 2 K2 light curve (in black) phased on the 6.6, 9.0, and 24.2 d periods (from top to bottom). The phase-binned mean light curve is shown in red.

3.2 LCOGT: Multi-wavelength optical photometry

The resulting 𝑔′r′i′ light curves of CI Tau obtained for the first and second epochs are shown in the top two panels of Fig. 3¹. For the first (second) epoch, the peak-to-peak amplitude of variability reaches 0.94 (0.98), 0.67 (0.71), and 0.52 (0.54) mag in the 𝑔′r′i′ filters, respectively.

We compute an L1 periodogram of light curves from both semesters (Fig. 3). We see clear peaks at 11.6 d and 27 d in semester 20B, with log₁₀ FAPS of −22 and −24, respectively. While in semester 22B/23A, peaks are seen at 11.6 d, 25.8 d and 47.1 d, with log₁₀ FAPS of −6, −20 and −23, respectively.. These periods are consistent within their respective errors with the 11.5 ± 0.5 d, 24 ± 5 d and 47.2 d peaks found in the K2 light curve. We choose to conveniently cite the long-term stable period at 25.8/27 d seen in the LCOGT light curve as 26.4 ±0.6 d, which represents the mean period derived from the two semesters, and the error indicating half the difference between the two periods. We highlight that GLS periodograms were generated for the LCOGT light curves as well, producing results similar to those obtained from the L1 periodogram. We therefore opt not to include results from the GLS periodogram in the paper, as done for the K2 light curve.

Interestingly, the 9.0 d, 6.6 d and 14.2 d peaks are not detected in the LCOGT light curve but are present in the K2 data, although the former is well sampled to detect this period. To examine the possibility of a sampling issue, we down-sampled the K2 light curve to mimic a mean sampling of 0.215 d as in the LCOGT light curve. We still recover the same periods in the binned K2 light curve as in the original which implies that this discrepancy in periods found in the two datasets is not due to the time sampling. In addition, the LCOGT data does not sample burst timescales which can be of order of less than a mean sampling of 0.215 d, while the K2 data does.

We note that the K2 data were taken around three to six years before the LCOGT data. A non-detection of the aforementioned periods in the LCOGT data might indicate that these periods are either transient or potentially drowned by other activity-related phenomena operating over a more extended period of several years. Alternatively, the physical processes responsible for these periods might have a lifespan shorter than the time span covered by the LCOGT light curves (3–6 months). These processes may occur randomly in phase within the given time frame, resulting in their lack of detection in the LCOGT periodogram.

Fig. 3 LCOGT light curves and L1 periogogram for each band, 𝑔 (green), r (blue), and i (red). Top panel: first semester (20B) light curve observed between 31 October 2020 and 08 January 2021. Second from top panel: second semester light curve taken between 27 September 2022 and 31 March 2023. Third from top panel: L1 periodogram of the combined dataset for semester 20B. Bottom panel: L1 periodogram of the combined dataset for semester 22B/23A.

3.3 Periodogram-peak significance-test in red noise

A challenging problem in observational astrophysics has been to detect periodic or almost-periodic signals in noisy time series. The task becomes even more challenging when we aim at detecting significant periodic signals hidden in stellar activity (e.g. spot modulation combined with red noise).

As seen in Fig. 4, the power spectrum of CI Tau’s K2 light curve (in black) consists of a red noise component which causes an increase in power at lower frequencies (red continuous line). To test the significance of the periodogram peaks in such red noise, we adopt the recipe of Vaughan (2005). In particular we test the significance of the 24 d peak which is of main interest for this study.

The method involves fitting the underlying power spectral density (PSD) shape of the classical periodogram (Schuster 1898) and identifying any coherent, statistically significant peaks by computing false alarm probability (FAP) levels. The normalised periodogram is given by, (1)

where ΔT is the sampling of the time step for an n evenly sampled time series x_n. For the K2 light curve of CI Tau, ∆T corresponds to ~30 min cadence. The periodogram (Schuster 1898) of the time series is the modulus-squared of the discrete Fourier transform, X_j, at each of the independent Fourier frequencies, with N being a normalisation factor. Only frequency points are considered untill the Nyquist frequency, and frequencies higher than the Nyquist are ignored.

We note that since the classical periodogram of the K2 light curve (Fig. 4) is a variation of the GLS periodogram computed in Sect. 3.1, we do expect slight differences between periodogram peaks obtained from the two methods, but in practice the classical periodogram is accurate enough to give useful insight into important periodic features (e.g. VanderPlas 2018). For example, the 24 ± 5 d period peak obtained using the GLS method (Fig. 1) is fully consistent within errors with the 25.7 d peak found in the classical periodogram (Fig. 4).

At any given frequency f_j, the discrete Fourier transform I(f_j) is scattered around the true power spectrum, P(f_j), following a χ² distribution with two degrees of freedom, the mean and variance: (2)

where, χ² has an exponential probability density (with a mean of 2 and variance of 4) of the form: (3)

We assume that the underlying power spectrum of the time series follows a power law of the form P(f) = Nf^−α. We use a least-squares (LS) fitting method to fit the continuum of the log-periodogram in linear-space and extract and , which are parameters governing the red noise continuum model of the form (see Vaughan 2005 for full description of the process). We use the null hypothesis that the data were generated by a process with the model spectrum and no periodic component.

To infer the statistical significance of the 25.7 d peak, we use an analytic significance test based on the χ² statistics, accounting for the number of frequencies sampled. Since the residual, γ(f_j) = 2I (f_j)/P(f) is χ² distributed, the integrated area under the probability density function of the χ² distribution (Eq. (3)) gives the probability that the power associated with a given peak is significant to a limit (1-∊), where ∊ is the FAP level. When considered in units of the periodogram power, it is given by, (4)

Once ∊ is specified, γ_∊ is calculated to give the significance level used to identify outliers in the periodogram that indicate the presence of statistically significant peaks based on the FAP.

We compute the 0.1 and 0.01% FAP levels (red dotted lines, in Fig. 4). We find that the 25.7 d period peak in the classical periodogram (corresponding to the 24 ± 5 d in the GLS periodogram) is significant beyond the 0.01% level (see inset panel of Fig. 4) and therefore is a significantly true period.

Fig. 4 Power spectrum of the photometric time series (black line) fitted with a power law of the form P(f) = Nf−α (red continuous line). The red dotted lines are at the 0.1 and 0.01% FAP levels from bottom to top. respectively. The 25.7 d period is clearly significant above the 0.01% FAP.

3.4 Gaussian process regression using a quasi-periodic kernel

Rotationally modulated phenomena such as stellar spots, plages and faculae, induce signals in both radial velocities (e.g. Rajpaul et al. 2015) and photometric time series (e.g. Notsu et al. 2013). If these phenomena are regular and stable over a long period of time with a given periodicity, we expect to be able to detect a clear periodicity related to these modulations in both RVs and photometry. However, stellar spots and faculae are not static, but rather evolve over time, emerging, changing shape, and have a finite lifetime. This results into non-periodic/quasi-periodic (correlated-noise-like) variations over each cycle of the stellar rotation, with shape and amplitude of the variability changing from cycle to cycle as seen in CI Tau’s light curve. Inferring rotation periods using straightforward sinusoidal variability models are not suitable in such cases since these only account for white noise. An alternative is to use a model that employs stochastic variability such as a Gaussian process (GP, Angus et al. 2018).

In summary, a GP is a stochastic process that captures the covariance between observational data points and allows for the modelling of correlated (red) noise. A GP is built using a covariance matrix with individual observation variances as diagonal elements and each off-diagonal element describes the covariance between any two observations. A kernel function is used to determine the values of the off-diagonal elements, which encapsulates the nature of the correlated noise. GPs provide a great deal of flexibility that has made them an effective tool to account for stellar activity (Haywood et al. 2014). We employ a GP regression method as prescribed by Rajpaul et al. (2015) to model the effects of stellar activity in the data (spots and random processes) that enables us to gauge the robustness of the 24.2 d period in the light curve. The recommended kernel for stellar activity characterisation is a quasi-periodic (QP) kernel (Aigrain et al. 2012) of the form: (5)

The QP kernel uses four parameters (termed hyperparameters) to model the stellar activity: h, P, λ_p and λ_e. t and t′ are times of observations made at subsequent intervals. P and λ_p represent the period and smoothing parameter of the periodic component of the variations. The parameter λ_e denotes the evolutionary timescale. The parameter h regulates the GP model amplitude. Additionally, we incorporate a white noise term as , where σ_K2 represents the noise amplitude, and δ(t, t′) is the Kronecker delta function.

We employ MCMC sampling to model both the 6.6 d and 9.0 d periods simultaneously, utilising a combination of two QP kernels, each corresponding to these distinct periods. The fitting process is performed using MCMC sampling, implemented through the Python package EMCEE. To explore parameter space, we employed a total of 256 walkers and executed 10 000 iterations during the MCMC sampling process. The initial positions for these walkers were set to correspond to initial parameter guesses. We guess λ_e through a visual examination of the data to identify a characteristic timescale governing variations in the time series – such as changes in shape, amplitude, and phase. Our estimation points towards a value of approximately 20 days, which was fixed in the fitting process. For the remaining hyper-parameters we assume uniform priors, as detailed in Table 1. The best-fit hyperparameters and their corresponding errors are summarised in the same table.

The GP model is shown as a red curve in the top panel of Fig. 5. The L1 periodogram of the residual light curve after removing the modelled periods still shows a clear peak at 24.3 d with a log₁₀ FAP of −58 (Fig. 5, bottom panel). This indicates that the 24 d period is not related to any of the two periods at 6.6 and 9.0 d periods, i.e. not a synodic period, consistent with the findings of Biddle et al. (2021). We also test the possibility of the 24 d period being the beat period between either of the 6.6/14.2 and the 9.0 d. We simulate 1000 light curves for each period combination, with period peak amplitudes and red noise levels that mimic the original K2 light curve. These light curves are then summed by continuously phase-shifting one of them to cover all phases. We then compute FAP levels on the normalised powers using all simulations. The results suggest that at any given phase, adding any of these two periodic models does not give rise to a beat period peak in the periodogram at 24 d as high as the one seen in the original K2 light curve.

Fig. 5 Cleaning of stellar activity signals from the K2 light curve. Top panel: quasi-periodic GP model fit (in red) with two periodic components of 6.6 and 9.0 d. The shaded region in the top panel shows 1, 2, and 3σ from the mean model. The bottom panel shows the L1 periodogram of the residual light curve after removing the modulation due to stellar spots.

4 Spectroscopic analysis

4.1 ESPaDOnS RV analysis

The RV measurements² were obtained by fitting Gaussian functions to each mean LSD profile, using the python LMFIT package (Newville et al. 2016). An L1 periodogram of the full 72 RV data points is shown in Fig. 6. Two significant periods in the whole ESPaDOnS RV dataset are seen at 9.27 d (log₁₀ FAP = −11) and 6.9 d (log₁₀ FAP = −10), both of which are also seen in the K2 light curve.

In addition, we show a 2D periodogram of the LSD profiles, for the 72 ESPaDOnS spectra in Fig. 7. The 2D periodogram is useful to discern periodicities in each velocity channel of the LSD profile. To achieve this, we divided the line profiles into smaller velocity bins. Subsequently, periodograms were independently computed for each velocity bin. The resulting individual periodograms are visualised in a 2D representation as a function of each velocity bin, with the power of the periodogram depicted by colours (e.g. Alencar & Batalha 2002; Sousa et al. 2021). From Fig. 7 it is clear that both the 9 and the 24 d periods seen in the K2 light curve, are also seen to be present in the LSD mean profiles of ESPaDOnS.

4.2 ESPaDOnS bisector analysis

Stellar activity such as spots and faculae induce variability in photospheric lines, which in turn has an effect on RV measurements (e.g. Reiners et al. 2013). A first-order effect is seen when the spots occult either the receding or approaching half of the visible stellar hemisphere, modifying the shape of the spectral lines. This causes a apparent shift in the disk-averaged RV (see e.g. Donati et al. 2020a) that can mimic exoplanet reflexmotion signatures, adding confusions to a true Doppler signal due to a planetary companion that may be present in the system (Queloz et al. 2001; Desidera et al. 2004; Huélamo et al. 2008; Carolo et al. 2014).

The bisector span (BIS span) measures the asymmetry of the mean line profile as was shown in a seminal study by Queloz et al. (2001). This method is a well known proxy to help gauge whether observed RV variations are due to stellar activity (e.g. Queloz et al. 2001; Meunier et al. 2010; Dumusque et al. 2014; Donati et al. 2020a) or otherwise likely induced by Doppler reflex-motion signal by, e.g a planet (Moutou et al. 2005; Naef et al. 2007). If the asymmetry of the mean line profile is purely anti-correlated with the estimated RV, then it is a sign that the line profiles are being modulated by spot activity only. Otherwise, a flat BIS span with respect to the estimated RV indicates that the shift in the lines are purely due to Doppler reflex motion (Queloz et al. 2001). In a mixed scenario when we have both a Doppler reflex by a planet and the effect of activity due to spot, both inducing RV signals that are comparable, we expect both an anti-correlation and a translational horizontal scatter of the BIS span (see Simola et al. 2019, and see also Sect. 4.4, for a detailed analysis).

We carried out a BIS span analysis using the mean LSD profiles of CI Tau. The ESPaDOnS data is of particular interest because at optical wavelengths, the contrast between the spot and the stellar continuum is largest, which makes any stellar activity due to spots more detectable. We followed the standard recipe for BIS span computation in Queloz et al. (2001). We divided each mean LSD profile in 100 bisector vertical intervals and computed the mean of each interval to obtain the bisector. The mean velocity of the top and bottom 10% was computed to obtain the BIS span given by (see Queloz et al. 2001, for details). A plot of the BIS span as a function of the RVs derived from the LSD profiles is shown in Fig. 8.

We find that there is an anti-correlation between the BIS span and RVs at velocities lower than ≲ 17.0 km s⁻¹, after which the BIS span reaches a plateau (Fig. 8). The data points are colour-coded according to the phase using the ephemeris: BJD0 (d) = 2 457 736.76 + 9.01c, where c denotes the spot rotation cycle, starting from an arbitrary initial BJD0 of 2 457 736.76, which corresponds to the first ESPaDOnS observation.

The BIS span is observed to be dependent on the phase of the spot rotation period at 9.0 days, consistent with findings reported by Donati et al. (2020a) and not very different from the newest value of 9.01 d reported in Donati et al. (2024). Between phases ~0.5 to around 1.0, the spot seems to be in a configuration in which it induces a noticeable effect on the overall mean line-shape, causing an apparent RV shift which is intrinsic to the star’s spot activity. During the rest of the phases, the effect of spot activity becomes minimal and the BIS span is flat. A flat BIS span suggests no correlation with RVs and therefore RV shifts of non-intrinsic nature (e.g. Udry et al. 2000). CI Tau’s BIS span therefore indicates a mixed scenario whereby stellar activity due to spot, is seen, but there is also a second component not related to spot modulation, where pure RV shifts are detected in the lines.

In Fig. 9, we display the L1 periodogram of the whole RV dataset and RVs related to BIS span values below and above 0.1 km s⁻¹, which contain 45 and 27 data points, respectively. We arbitrarily chose this threshold because below this limit the BIS span is flat, although we could not exclude some contamination from spot modulation at RVs in the range ~ 16.2–17.0 km s⁻¹. Our selection criteria results in 45 RV data points span over 64 days, which correspond to BIS span values that are essentially in the flatter part of the BIS span plot. We notice that the ~24.6 d period peak is suppressed during spot activity, but rises to a significant FAP level, when RVs corresponding to only the flat part of the BIS span are considered.

Fig. 6 L1 periodogram of RVs related to 72 data points (four nightly sub-exposures) from ESPaDOnS. The spot period is recovered at 9.27 d and another period is recovered at 6.9 d.

Fig. 7 2D periodogram and mean LSD line profiles using 72 ESPaDOnS spectra. Top panel: 2D periodogram with colour range corresponding to the power of the periodogram. The horizontal lines are drawn at 9.0 d and 24.2d. The contours represent FAP levels at 0.01% (black), 0.1 (grey), and 1% (white). Bottom panel: mean LSD profile in black, with one standard deviation of all profiles shaded in blue.

Fig. 8 Bisector span of the mean LSD profiles from ESPaDOnS spectra. A correlation is seen at lower velocities up to around 17 km s−1 and then no correlation where it reaches a plateau. The data points are colour-coded based on phase of the 9.01 d spot rotation period. The dashed line indicates the threshold below which BIS span values were selected to compute the L1 periodogram shown in red in Fig. 9.

Fig. 9 LI Periodogram of RVs for the whole ESPaDOnS RV dataset (black), RVs selected based on BIS span values of >0.1 km s−1 (blue), and RVs selected based on BIS span values of <0.1 km s−1 (red). The level of activity due to stellar spots, corresponding to 6.9 and 9.27 d (vertical green-dotted lines) is high in the whole RV dataset and RVs corresponding to BIS span values of >0.1. For the flatter part of the BIS span with values <0.1, the RVs display a clear period peak at ~24.6 d (vertical red dotted line).

4.3 SPIRou RV analysis

Activity signatures induced by spots are known to be strongly chromatic (e.g. Reiners et al. 2010), since the spot-to-stellar continuum contrast tends to be lower at longer wavelengths in the NIR. In a recent study by Miyakawa et al. (2021), it was shown that chromaticity can actually be even more pronounced if we consider wavelengths from the visible (K_p, Kepler band covering 400–900 nm) to the NIR. Miyakawa et al. (2021) studied RV jitter due to stellar activity on RVs in four different passbands (K_p, J, H and K_s). The authors found that RV ‘jitter’ can vary from 1.4 km s⁻¹ down to 0.16 km s⁻¹ from K_p to K_s bands. This has important implications, especially for planet searches via spectroscopy in very young and active stars like CI Tau, where the effects due to spots may induce apparent RV shifts in the CCFs that are comparable to or even higher than pure Doppler RV shifts that can be induced by a planetary companion (see Sect. 4.4).

We utilised the SPIRou-CCF package, which employs a cross-correlation method to determine RVs³. The input data for our analysis are the telluric-corrected spectrum, which are calculated using the APERO pipeline (Cook et al. 2022). The SPIRou-CCF package includes several processing steps to reduce the significant systematic effects that are commonly found in near-infrared (NIR) data. The cross-correlation process has been described in detail in Martioli et al. (2022).

To test the chromatic effect of the spot on the RVs, we explored two separate cases. In the first case, we constructed a template spectrum (T1), adapted to the spectral type of CI Tau. To generate a template, we formulate a mask by selecting spectral features based on their proximity to the spectral type of the star. Each mask comprises a collection of around 1844 atomic and molecular lines in the whole SPIRou domain, with central wavelengths sourced from the VALD catalog (Piskunov et al. 1995). The line depths are empirically derived from the template spectra of luminous stars observed by SPIRou.

Each individual telluric-corrected spectra were then cross-correlated with T1 to obtain mean cross-correlation function (CCF) profiles of the absorption lines from which the RVs (which we refer to as RVs_full) were extracted. For the second case, to make sure we are least affected by spot activity, we considered a wavelength range as red as possible in the SPIRou domain, making sure we had enough photospheric lines for CCF computation and also making sure that there are not magnetically sensitive lines (e.g. Lavail et al. 2019). Upon close inspection, we found the 2293–2340 nm range in the CO band matched this selection criteria the best. We constructed a sub-template (T2) from T1, over the 2293–2340nm wavelength range. We then cross-correlated all spectra with T2 and the RVs (which we refer to as RVs_SUb) were extracted from this domain separately. Figure 10 shows the L1 periodogram of the CCF RVs for RVs_full and RVs_sub It can be seen that the peak corresponding to the 25.2 d period is significant for RVs_sub compared to RVs_full⁴.

Our analysis aligns well with the results of Reiners et al. (2010); Miyakawa et al. (2021), which outline the chromaticity of the spot signal in the RVs. Therefore, to make sure we are least affected by this effect, we carried out our SPIRou orbital analysis considering RVs extracted over the 2293–2340 nm range. To make sure that the chosen lines in this region are photo-spheric, we determined the veiling value of CI Tau for a given spectrum (S) as described in Sousa et al. (2023). We then computed the residual spectrum, by subtracting S from a template veiled to the veiling value we computed in the K-band. The residual profiles show no feature at the location of the photospheric lines, which indicates that the lines are purely photospheric. A plot of all the CCF CO line profiles and their mean is shown in Fig. A.1.

Fig. 10 L1 periodogram of CCF RVs extracted from the full SPIRou wavelength range (in black) and CCF RVs related to the 2293–2340 nm wavelength range (in red). The dotted vertical lines mark the significant peaks based on their FAP, one of them being the 25.5 d peak.

Fig. 11 Observed (red) and modelled (blue) RVs vs. BIS span for two different cases. Case 1 (left panel): Gspot modulated at 9 d and Gstar kept fixed and case 2 (right panel): spot modulated at 9 d and star Doppler shifted with semi-amplitude of 0.3 km s−1 at 25.2 d period. Individual point error bars are also plotted.

4.4 SPIRou bisector analysis

We use the CCFs obtained by cross-correlating the full SPIRou wavelength range to compute the bisector span (BIS span) values as we described in Sect. 4.2. A Spearman correlation test on the SPIRou RVs against the obtained BIS values, yields an intermediate anti-correlation coefficient of −0.47 ± 0.05. The correlation coefficient uncertainty is obtained using a Monte-Carlo simulation on the individual RVs and BIS values including their errors.

In an attempt to investigate what scenario would give rise to the observed SPIRou bisector span, we have devised a simple model to demonstrate the impact of a spot on the bisector under two circumstances: case (1) when the central star remains stationary and case (2) when it undergoes a Doppler shift caused by the presence of a companion. In our model, we conduct experiments to investigate how the presence of a spot affects the disk-averaged star CCF and its impact on the bisector. We make the assumption that both the star and spot CCF are Gaussians, denoted by G_star and G_spot, respectively. We take a realistic approach by considering G_star to be the mean profile of all available SPIRou CCFs. To achieve that, we fit Gaussians to all of the available 425 SPIRou mean CCF profiles. The Gaussian fit parameters, amplitude (A_star), mean (μ_star) and standard deviation (σ_star), together with the associated error in the respective fit parameters are recorded.

CI Tau’s spot was studied in detail by Donati et al. (2020a) through Zeeman Doppler imaging (ZDI). It was illustrated that CI Tau’s surface contains an accretion region characterised by a bright chromospheric region that aligns spatially with a dark and cool spot. Donati et al. (2020a) were able to constrain the spot’s surface area to be approximately 20% of the stellar surface. Since the impact of this cold spot on the disk-averaged spectrum is less pronounced than that of the unspotted portion of the stellar surface, the presence of a cool spot manifests as a distinctive feature, resembling a bump, which traverses the spectral line profile from the blue to the red end (Rosén & Kochukhov 2012; Hébrard et al. 2014). To simulate this bump, we incorporated the effect of the cold spot on the disk-averaged CCF of the star. We assumed specific parameters for G_spot that are, amplitude (A_spot), mean (μ_spot), and standard deviation (σ_spot). Since the spot’s contribution to the overall CCF is phase dependent from an observer’s point of view, we modulated A_spot in such a way that it had a lesser effect on the wings of G_star.

In the first case, we sinusoidally modulate G_spot with a ~9d period, keeping G_star stationary at a mean value of ~16.39 km s⁻¹, corresponding to the observed mean radial velocity of CI Tau. The BIS span is computed from the combined CCF models (star + spot) and a Spearman coefficient of approximately −0.49 is calculated from the model (Fig. 11, left panel). In this case, since the star is not Doppler shifted, the only effect of the spot is to rotationally modulate the bisector lines about the mean velocity of the star. The RV versus BIS span plot results in a distinctively ordered shape as can be seen (similar to Hébrard et al. 2014), and in this case we are unable to statistically reproduce the observed BIS span versus RVs plot of CI Tau (as discussed below), represented by red circles in Fig. 11.

In the second scenario, we sinusoidally modulate G_spot with a ~9 d period and we incorporate a periodic Doppler shift into G_star, characterised by a 25.2-day period and a semi-amplitude of 0.3 km s⁻¹ as determined from our Keplerian orbital fit (see Sect. 4.5). A translational Doppler shift introduced into G_star due to Doppler reflex induced by the planet, induces a horizontal scatter in the modelled BIS span, which results in a better agreement in the overall shape of the observed BIS span (Fig. 11, right panel). We experiment for different parameter values for the spot’s Gaussian in our model and we find that a value of ~0.7 for the ratio (σ_spot/σ_star) and ~0.1 for the ratio A_spot/A_star, gives the best model match to the observed BIS span in this case. We also experiment with different initial phases for the periods with which G_spot and G_star are modulated. The results are independent of phases.

By employing the relationship between the radial velocity (RV) amplitude of the spot, its latitude (l), and the star’s rotation (∆V_rad = 2υ sin(i) cos(l), as established by Bouvier et al. 2007 and McGinnis et al. 2020), and taking into account a value of 9.5 ± 0.5 for υ sin(i) as reported by Donati et al. (2024), our analysis suggests that the latitude of the cold spot should be close to zero, aligning with the best-fit solution observed in the RV versus BIS span plot. This inference leads us to hypothesise that the observed spectral line-shape distortion attributed to cold spot activity may result from a combination of contributions, not only from the high-latitude spot (Donati et al. 2020a), but potentially from other equatorial cold spots as well, as suggested by Donati et al. (2024).

Our model illustrates the impact of a translational shift in the star’s CCF when accompanied by spot activity, a phenomenon previously examined by Simola et al. (2019). Additionally, as Hébrard et al. (2014) demonstrated, the signature of spots exhibits a distinctive S-type pattern. We can only reproduce this pattern in case 1, where the star remains fixed and only the spot modulates. However, the observed BIS span versus RV plot exhibits significant scatter around the mean, almost resembling an ellipse, even if the respective errors are considered. We are able to replicate this behaviour to a considerable extent when the star’s centre of mass is in motion in tandem with the spot’s motion.

Since in both cases the spearman coefficients of the model and data are not very different from each other, although the shapes of the modelled BIS span versus RV plot differs, we attempt to quantify the shape and the degree of scatter using a simple statistical method that computes a total binned variance (TBV) of the BIS span at binned RVs. For the observed data we compute a TBV 3.54 km s⁻¹. We carry out a Monte-Carlo simulation in which we consider uniformly distributed random values for the depth, mean and standard deviation of both the spot’s and star’s Gaussian profiles. For the latter, the values were uniformly distributed between the errors we obtained from the mean-profile fit of the SPIRou CCFs. We run 10000 iterations in each case (case 1: fixed star and case 2: Doppler shifted star), generating probability density functions (PDFs, Fig. 12) of the TBV in each case. The results strongly indicate that the probability of obtaining the observed shape of the BIS span versus RV plot solely through spot activity is exceedingly low (≤1%). Our analysis suggests the observed BIS span versus RV plot shape agrees best with a scenario in which we have combined effect of spot activity and Doppler reflex motion induced by a companion.

We also experiment on a scenario discussed by Donati et al. (2024), which speculates that the 25.2 d period could be a result of a non-axisymmetric structure in the inner disk at around 0.16 au with a line profile contribution on the CO lines. This disk material with a Keplerian velocity of 57.4 km s⁻¹ could create a line profile contribution of around 3.8 times larger in FWHM in the CO photospheric lines compared to the observed LSD profiles and is expected to modulate RVs at a 25.2 d periodicity. Our veiling measurements unequivocally demonstrate that the CO lines have a photospheric origin. Specifically, when we subtract a template veiled to match the level of veiling observed in CI Tau from the original spectra, the residuals reach the noise level at the location of the CO lines. Assuming that there is still some undetected contribution from a non-axisymmetric material at 0.16 au, we test the scenario by adding this contribution to the total line profile (i.e. spot modulation at 9 d period + line profile contribution due to structure modulated at 25.2 d + fixed star). Our computed the TBV (shown in green distribution in Fig. 12) suggest this is statistically a less probable scenario to account for the observed BIS versus RVs plot presented in Fig. 11.

Fig. 12 TBV vs. probability density after 10 000 MC simulations for three cases. In case 1 (orange) the star is at rest and the spot modulated at 9 d; in case 2 (blue) the star is Doppler shifted with a semi-amplitude of 0.3 km s−1, keeping all other parameters the same as in case 1. In case 3 (green) we test a scenario proposed by Donati et al. (2024), i.e. a contribution from material in the inner disk on the CO lines, that could cause the 25.2 d periodicity. The distributions represent the TBVs of the BIS span plots. The grey dotted line is the observed TBV for the SPIRou BIS span vs. RVs plot.

4.5 Combined SPIRou and ESPaDOnS RVs orbital analysis

We carry out our RV orbital analysis on the combined ESPaDOnS and SPIRou RVs, spanning over approximately eight years. We show the L1 periodogram of this dataset in top panel of Fig. 13. Dominant periodicities are seen at 9.0 d, 11.7 d and 25.2 d⁵.

We conducted a joint Markov chain Monte Carlo (MCMC) fitting process for both the stellar activity spot and Keplerian orbital components using the combined RV data. The MCMC procedure was implemented using the Python package EMCEE and the full fitting process is described in Appendix B. We set uniform priors for all the QP-GP hyperparameters: GP amplitude (h), spot period (P_spot), smoothing parameter , and evolutionary timescale . The Keplerian orbit parameters include the orbital period (P), eccentricity (e), semi-amplitude (K), time at periastron passage (T₀), and the argument of periastron (ω). We incorporate a zero semi-amplitude orbit model into our MCMC priors, specifically to accommodate cases where only a spot GP model is applicable. An observation close to middle of the data were chosen as initial guess for T₀. A Lucy and Sweeney test (Lucy & Sweeney 1971) was performed on the residuals between a circular model and an eccentric model fit to the data. We find a p-value of 0.000851, which suggests an eccentric orbit is preferred over a circular one.

The priors and best-fit parameters resulting from the MCMC fit are listed in Table 2. In Fig. 13, the best-fit Keplerian model is overplotted in phase with the RV data cleaned from the spot GP model.

We use the best-fit Keplerian parameters to determine the mass function given by: (6)

where G is the gravitational constant. Since different values have been reported for the central star mass (0.8 ≤ m₁ ≤ 0.9 M_⊙, Guilloteau et al. 2014; Simon et al. 2019) and disk inclinations (50° ≤ i ≤ 70°, Clarke et al. 2018; GRAVITY Collaboration 2023), we determine the range of possible values for the companion mass, assuming the orbital plane is aligned with one of the disks. The error on the companion mass was obtained considering the taking the entire inclination and mass ranges together with the errors obtained on the Keplerian orbital parameters listed in Table 2. Our result is consistent with a ~3.6 ± 0.3 M_jup planet in a highly eccentric orbit of eccentricity and at a semi-major axis of ~0.17 ± 0.08 au from the central star.

4.6 Planet RV detection limit in spot activity

We conducted experiments to simulate the detection limit of the periodogram peak corresponding to an injected planet with a 25.2-day orbital period, in the presence of varying levels of spot activity. To achieve that we utilised Gaussians constructed for the star (G_star) and spot (G_spot) described in Sect. 4.4.

To mimic the Doppler shift caused by the injected planet onto the star, we introduced a periodic shift of 25.2 days on G_star. This shift was maintained constant at a semi-amplitude value of K = 0.3 km s⁻¹, which was determined through a Keplerian orbital fit (Sect. 4.5). We then modulate G_spot with a period of 9 d and an amplitude (h). At each modulation step, we compute a total profile (G_star + G_spot) and record the modulation steps as time. We fit a Gaussian to the total profile each time to extract the RV (as we would do for a normal observed CCF). The procedure was executed for various values of h, yielding a 1D periodogram for each h/K. The combination of these individual results produces a 2D periodogram, as illustrated in Fig. 14. Subsequently, we determine the threshold at which the peak corresponding to the planet can be consistently detected in the periodogram. Various initial phases between the spot and the injected planet were investigated during this exploration.

Our simulation results indicate that a clear detection of the period peak associated with the injected planet occurs when the ratio h/K is approximately 1.0 or less (see Fig. 14). The planet’s signal is only marginally discernible for h/K values between ~ 1.0 and 2.0 and almost not detected for h/K values higher than ~2.0. For CI Tau, our analysis yields an h value of 0.34 km s⁻¹, determined through GP regression of the spot model in the RVs and K = 0.3 km s⁻¹. This leads to an h/K ratio of ~1.1, where the RV signal arising from spot rotation slightly predominates over the Doppler signal from the injected planet. Nevertheless, the planet’s signal can still be marginally detected (indicated by the white dashed line in Fig. 14). This provides a potential explanation for the limited detection of the 24 d peak in the SPIRou RVs (Fig. 10) that were derived from spectral lines spanning the entire SPIRou spectral range, encompassing bluer wavelengths that are particularly susceptible to the influence of spot activity.

In addition, as shown by Miyakawa et al. (2021), activity levels exhibit chromatic behaviour and can affect optical lines up to a level of ~4 times more than in the red. If we consider spot RV signal to be 0.34 km s⁻¹ from our GP fit to CO RVs and four times more in the optical, i.e. an h/K value of ~4 and will be most likely not detected in lines most affected by spot activity. Donati et al. (2024) demonstrate a approximately 1.5-fold increase in activity levels from the K-band (where the CO line are) to the H-band, within the SPIRou domain. Considering this factor, we anticipate h/K to be 1.7 in the H-band. As depicted in Fig. 14, this is likely why the planet’s RV signal will be hardly detected in the H-band.

Moreover, in the radial velocities obtained from ESPaDOnS, a comparable situation is observed. As depicted in Fig. 8, the BIS values span a range corresponding to an amplitude of 2.5 km s⁻¹, encompassing both spot activity and the planet’s signal. Taking half of that value implies an h/K of 4.2, possibly explaining the limited detection of the planet’s RV signal in the ESPaDOnS data.

Interestingly, our 2D periodogram reveals faint appearance of period peaks at approximately 11.6 days, 40 days and longer periods exceeding approximately 60 days. These peaks become more noticeable when h/K values are below approximately 1.0, where the injected planet’s and spot’s signals are comparable. We suspect that these peaks could possibly be beat periods arising from the combined influence of the 24 d, 9 d peaks and other secondary beat peaks that might arise from any of these periods.

We emphasise that our model solely illustrates the collective impact of both a spot and an introduced planet on the overall CCF and BIS span. One of our key assumptions is that the spot’s influence on the entire CCF profile can be approximated by a Gaussian distribution. However, in reality, it can possess a more complex structure, as demonstrated in Donati et al. (2020a). Furthermore, our model does not account for other quasi-periodic or stochastic effects stemming from magnetic fields or other factors related to stellar activity or accretion, which could potentially introduce additional features in the BIS span and CCF RVs and which could make the planet’s signal detection even harder at bluer wavelengths.

5 Discussion

5.1 Periods in light curve

Our detailed analysis of the K2 and LCOGT light curves show periodicities at 6.6, 9 (spot), 11.5, 14.2 and 24 d/26.4d. We examine whether this 6.6 d period can be understood through an unstable accretion state for CI Tau. Using 3D numerical simulations, Blinova et al. (2016) have shown that the accretion state (stable or unstable) can be determined by computing the ratio of the truncation radius to the corotation radius, r_t/R_cor. They estimated that accretion is unstable if r_t/R_cor ≲ 0.7. Using literature values (e.g. Donati et al. 2020a) for the mass-accretion rate, dipolar magnetic field strength, stellar mass and radius of CI Tau, we find that r_t/R_cor ≈ 0.3−0.57⁶. This value of r_t/R_cor indicates that accretion in CI Tau is unstable. As shown in Blinova et al. (2016), unstable accretion states can produce significant periodogram peaks at values lower than the stellar period. In addition, such peaks were observed, through wavelet analysis, not to be stable over the entire duration of the simulations (see Figs. 3 and 5 of Blinova et al. 2016). Hence, we hypothesise that the 6.6 d period identified in both the K2 light curve and in some parts of the RVs, could be attributed to the influence of this specific accretion regime.

A slightly longer period than the spot modulation is seen to be coherent at 11.5 and 11.6 d in the K2 and LCOGT light curves, respectively. The exact origin of this period is still not clearly understood, but there are clues that it could potentially arise due beating effects between the 6.6 and 14.2 d peaks. In the same way, the 14.2 d period peak could well be a beat period between the 9 and 24 d peaks. One other mechanism giving rise to these peaks could be time-variable accretion. In CTTS, the mass-accretion rate can vary by up to a factor 0.4 dex on timescales of days to weeks (e.g. Costigan et al. 2014; Zsidi et al. 2022). Furthermore, recent studies (e.g. Venuti et al. 2021) have demonstrated that non-steady accretion onto CTTS is expected to produce periodic signals in light curves at periods larger than the stellar period. This particular CTTS feature is further supported by numerical simulations of star-disk interactions (e.g. Romanova et al. 2005; Zanni & Ferreira 2013), which have shown that magnetospheric accretion is not a steady-state process and can occur in quasi-periodic cycles with timescales longer that the stellar period. In summary, periods shorter than the stable long period of 24 d, can be related to magnetospheric accretion processes, that are hard to interpret.

The only stable period that seem to be consistently detected in both the K2 and LCOGT light curves and with high statistical significance, is the longer period at 24 ± 5 d in the K2 data and at 26.4 ± 0.6 d in the LCOGT data, both of which are consistent within their respective errors. Teyssandier & Lai (2020) presented hydrodynamical simulations of hot Jupiters orbiting protoplanetary disks, exploring how a giant planet’s mass, orbital separation, and eccentricity regulate the accretion rate at the inner disk edge. They found that a massive, eccentric planet can drive pulsed accretion at the inner edge of the disk, modulated by the planet’s orbital frequency. The amplitude of accretion variability generally increases with the planet mass and eccentricity in their model. They specifically studied CI Tau, showing that the observed luminosity variability in the star can be explained by accretion signatures modulated by an eccentric planet at 9 d. However, since it has been well established that the 9 d is due to spot modulation and not a planet, it is equally possible that an eccentric planet with a longer period modulates the accretion onto the star. This aligns well with our results, which indicate an eccentric and massive planet is indeed present in CI Tau at an orbital period of 25.2 d.

An alternate mechanism that could account for the 24 d/26.4 d periodicity observed in the light curves may result from the fact that CI Tau is observed at a relatively high inclination, approximately 70 degrees.This may lead to variations in brightness due to circumstellar matter periodically obstructing the star’s photosphere (resulting in variable extinction events). This occulting material might be situated in the innermost region of the accretion disk, at a distance of ~0.16 au from the central star. We note, however, that distinct colour behaviours are expected for a CTTS with variability dominated by circumstellar extinction as opposed to accretion dynamics, and the specific colour properties of CI Tau clearly favour the latter scenario. To test this point, we followed the approach of Venuti et al. (2015), who explored the typical correlation trends between optical variability amplitudes and simultaneous blue-filter colour variability for TTSs dominated respectively by magnetic spots, accretion features, and recurrent extinction by circumstellar material. We used our LCOGT time series photometry, as well as archival U-to-R monitoring data (e.g. Herbst & Shevchenko 1999), to derive colour slopes for CI Tau over days to years timescales and compare them with the range of typical behaviours described in Venuti et al. (2015) for different types of TTS variables (see, in particular, their Fig. 7). This analysis revealed that CI Tau displays a colour behaviour consistent with the properties of accretion-dominated sources (and inconsistent with being extinction-dominated) on all timescales from a few weeks (comparable to the rotational timescales) to several years, as shown in Fig. 15.

Fig. 13 RV orbital analysis. Top: L1 periodogram of the full ESPaDOnS and SPIRou RVs. The 25.5 d period peak, the spot modulation period peak at 9.0 d, and the 11.6 d period peak, with their associated FAPs are indicated. Middle panel: best-fit QP-GP combined model (9 d spot + 25.2 d planet) to the RVs. Bottom panel: best-fit Keplerian orbit model (black curve) overplotted on RVs cleaned from the spot GP model. An arbitrary phase space shift of 0.25 has been introduced to enhance the clarity of the orbital variation. The red dots illustrate a binned median of the data points, with bin sizes chosen to maintain uniform point distribution within each bin. The error bars on the red dots represent the standard deviation within each bin.

Fig. 14 2D periodogram of simulated RVs to show the peak detection limits of an injected planet at 25.2 d period in the presence of spot modulation at 9 d. The value of h/K derived for CI Tau is shown with the white dashed line at 1.1. The cyan line shows the detection limit in the optical where we expect the activity RV signals to be ~ four times larger. The light orange dashed line at 1.7 depicts the h/K expected in the H band.

Fig. 15 Optical variability amplitudes and colour slopes of CI Tau, following the classification of Venuti et al. (2015). It suggests that CI Tau’s colour behaviour is consistent with accretion-dominated sources.

5.2 Periodic variability at 25.2 d seen in the spectroscopic RVs

Our bisector analysis shows that in the presence of solely the spot-induced 9 d modulation period, the CCFs will be modulated periodically in a way so that only a rotational effect on the bisector will be detected about the star’s mean RV. By introducing a periodic component of 25.2 d as a Doppler shift in the CCFs, with an RV semi-amplitude obtained from our Keplerian model, we notice that the observed bisector span is more accurately replicated.

We test detection limits of a 25.2 d planet in CI Tau by varying the amplitude of the spot signal while keeping the Doppler signal due to the injected planet constant. We demonstrate that when the Doppler shift caused by the planet is on the same order of magnitude as the RV amplitude of the spot, the planetary signal can be detected to a significant statistical level in the periodogram of RVs that are less influenced by spot activity. In addition, we note that RV Doppler shift induced by the planet will be of the same RV amplitude in spectral lines at all wavelengths, but the shift will be diluted at bluer wavelengths which are known to be most affected by the star’s activity.

We reduce the effect of spot activity in our RV analysis by using photospheric lines in the SPIRou wavelength domain as red as possible, since the spot contrast is known to decrease at longer wavelengths (e.g. Miyakawa et al. 2021). We then combine all RV data from ESPaDOnS and SPIRou to perform a global MCMC of the spot activity and the Keplerian orbit. Our comprehensive spectroscopic analysis, statistically supports the notion that the star exhibits a Doppler reflex motion with a 25.2 d period due to a massive and eccentric planet. We do not entirely dismiss other potential, although less probable mecha-nism(s). These mechanisms are considered less likely since they should induce a Doppler reflex in the spectral lines, a point we substantiate through our comprehensive bisector analysis.

5.3 Eccentric planet

Our best-fit Keplerian RV model indicates an eccentric planet with an eccentricity of . Given the highly eccentric massive planet located close to the actively accreting star, we expect a pulsed-accretion mechanism operating in CI Tau as proposed by Teyssandier & Lai (2020). Their models predict only an eccentric massive planet in CI Tau can produce a periodic variability in the accretion which will be imprinted on the star’s light curve, making it periodically brighter at phases when the planet is at its closest approach to the star. With the detection of the same coherent periodicity in both K2 and LCOGT light curves, we suspect a pulsed-accretion could be a likely mechanism to create such a modulation in the light curves.

To test this, we analyse the equivalent width (EW) of the Paβ line in the ExTrA data, which is indicative of accretion-induced activity. A distinct peak at approximately 27.0 d is detected (see Fig. 16). We note that the periodicity peak is highest in Paβ EW periodogram, likely because it is the strongest line among the three Paschen lines (Paβ, Paγ and Paδ). The detection of this periodicity in this line, potentially provides observational support for the planet-induced phase-dependent pulsed accretion model proposed by Teyssandier & Lai (2020).

An additional peak at around 18.5 d is detected in the EW of the HeI-10830 line (Fig. 16). This period closely aligns with a periodicity of approximately 19.5 d observed in the blue-shifted absorption component of the HeI-10830 line (see bottom panel of Fig. 16). We posit that this may be tracing periodic outflows with velocities ranging from ~30–200 km s⁻¹ within the system. Outflows such as stellar winds and magnetospheric ejections have shown similar velocities in numerical simulations of star-disk interactions (Romanova et al. 2009; Pantolmos et al. 2020; Bouvier et al. 2023).

Furthermore, the system exhibits a significant excess in the near-infrared (NIR) portion of its spectral energy distribution (SED). Muley & Dong (2021) conducted an SED analysis of CI Tau, concluding the presence of an inner dust-free cavity. This observation is further corroborated by recent interferometric imaging of the very inner region of the CI Tau disk (GRAVITY Collaboration 2023), revealing a resolved inner disk cavity with a radius of 21 ± 2 R_★ (0.18 ± 0.02 au). This discovery aligns with the findings of Muley & Dong (2021) and suggests that a close-in planet might have played a role in carving this inner gap. Recent investigations into the inner CO disk of CI Tau by Kozdon et al. (2023) indicate a CO disk break occurring at approximately 0.14 ± 0.08 au, in close agreement with the expected planet location. Hydrodynamical simulations performed by Debras et al. (2021) demonstrate that a Jupiter-mass planet, after migrating into a low-density gas cavity, can acquire a high eccentricity of 0.3–0.4. Furthermore, Baruteau et al. (2021) propose that a approximately two Jupiter-mass planet can achieve an eccentricity of 0.25 within a few million years (Myr), resulting in pronounced asymmetries in the gas surface density within the cavity.

Since CI Tau is only a 2 Myr old object, we do not expect an inner planet to have circularised in this lifetime. Recent models by Terquem & Martin (2021), aim at assessing the timescale needed for circularisation due to tides raised by a Jupiter-sized planet around a solar mass star. The authors find that tidal dissipation in a Jupiter-sized planet leads to a circularisation timescale of 1 Gyr for an orbital period of 3 days and this relation stays flat after orbital periods of 3 d. Their value is somehow higher than what was predicted in the models by Rodríguez & Ferraz-Mello (2010). The authors used numerical simulations to investigate the long-term tidal evolution of a single-planet system. In their model, they included variations in the semi-major axis and eccentricity of the relative orbit, taking into account the influence of both planetary and stellar tides across various planet types. They determine the critical value of eccentricity at which the stellar tide surpasses the planetary tide in a Sun-Jupiter-type system and they computed a circularisation timescale of at least 0.1 Gyr, which is around two orders of magnitude higher than CI Tau’s age.

6 Conclusions

Our comprehensive examination of the CTTS CI Tau has revealed a statistically significant periodicity of approximately 25.2 d within the system. Through an extensive analysis of bisector span, we have gathered evidence that this 25.2 d period is related to a Doppler shift in the RVs least affected by the effects of spot activity. Our analysis supports the existence of an eccentric planet with a mass of 3.6 ± 0.3 M_jup in orbit around CI Tau at a semi-major axis of 0.17 ± 0.08 au. We investigated alternative scenarios that could potentially induce such a periodicity in the system, but were unable to identify one. Furthermore, we have demonstrated that in an actively accreting system like CI Tau, spot activity signals might attenuate the Doppler reflex signal of a planet in RVs associated with bluer wavelengths, which would make the planetary signal hard to detect in that domain.

CI Tau’s photometric light curves from the K2 mission and LCOGT display other significant periods that range between ~6.6 d and ~47 d. The only period shorter than the 9 d spot rotation period is found to be a 6.6 d period, which is seen only in the K2 light curve and not in the LCOGT light curve that is taken ~ three to six years later. This period is also detected in the ESPaDOnS RVs, but is not very clear in the SPIRou RVs. We hypothesise that this period is not stable over the long term, and may be related to unstable accretion on the stellar surface. Other periods that are slightly longer than the spot rotation period are 11.5 d and 14.2 d, which we think could be related to either beating effects or non-steady accretion variability on the stellar surface, although we cannot reach any clear conclusions about these periodicities.

Our study outlines the difficulty of detecting inner embedded planet(s) in the inner ~0.1 au of young CTTS systems which are dominated by activity. We demonstrate that when working with light curves and RVs, one of the main limitations to the detection of exoplanets in active stars such as CI Tau is the different sources of variability induced by the accreting star itself rather than the instrument’s RV precision. Only planets that induce Doppler signals in the RVs that are of the same order as or higher than the spot activity signals in highly accreting systems like CI Tau are likely to be detected in lines that are least affected by the star’s activity. Our study indicates that, unless future research convincingly establishes an alternative explanation for the 25.2 d periodicity, our findings, in conjunction with recent independent studies of the inner regions of CI Tau, strongly support the planet hypothesis as the most plausible explanation. Importantly, this marks the first identification of an inner embedded planet in a CTTS, highlighting that even in systems with pronounced activity signals the detection of inner embedded planets is achievable.

Fig. 16 L1 periodogram analysis of LCOGT light curves and trailed plot of He I line. Top: time series of EW measured for He I, Paβ, Paγ, and Paδ lines. Middle: L1 periodogram of EW. A clear peak is detected in the strongest line tracing accretion (Paβ) at ~27.0 d. Additional peaks that are relatively significant are seen at 7.2 d, 9.1 d (spot), 11.6 d, 18.5 d (see text). Bottom: trailed plot of the He I line phased at a period of 19.5 d showing high-velocity outflows traced by the blue absorption component in the line. The colour bar represents the normalised flux.

Acknowledgements

R.M. would like to respectfully dedicate this study to his dearly departed mother (Uma Devi Manick), who sadly passed away around the time of detection. We express our gratitude to the anonymous referee for providing valuable suggestions that significantly contributed to the improvement of the paper. We acknowledge support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 742095; SPIDI: Star-Planets-Inner Disk-Interactions, http://www.spidi-eu.org). R.M. thanks Dr. Evelyne Alecian, Dr. Oscar Barragan, Dr. Pia Cortes-Zuleta, Dr. Anthony Soulain, Dr. Konstantin Grankin and Dr. Catherine Dougados for fruitful discussions on the analyses presented in the paper. R.M. also thanks Hugo Nowacki for providing the code for creating the 2D periodogram of the ESPaDOnS RVs. This research made use of NASA’s Astrophysics Data System; SciKit GStat (Mälicke et al. 2021), SciPy (Virtanen et al. 2020); NumPy (Harris et al. 2020); matplotlib (Hunter 2007); and Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2018). We thank the director of LCOGT for a generous allocation of discretionary time. E.M. acknowledges funding from FAPEMIG under project number APQ-02493-22 and a research productivity grant number 309829/2022-4 awarded by CNPq, Brazil. This work was also supported by the NKFIH excellence grant TKP2021-NKTA-64. This research was made possible through the use of the AAVSO Photometric All-Sky Survey (APASS), funded by the Robert Martin Ayers Sciences Fund and NSF AST-1412587. R.M. acknowledges routines used from the IvS repository of KU Leuven to create trailed plot shown in Fig. 16.

Appendix A CCF profiles of CO lines from SPIRou observations and their mean profile.

Each SPIRou observation consists of four subexposures to measure Stokes V, taken at different orientations of the polarimeter and used to compute the non-polarised and the circularly polarised profiles. We treat each sub-exposure as an individual spectrum and compute a mean CCF of photospheric lines in the CO region for each of them.

In Figure A.1, we show the average of all the CCF profiles in black dotted points. We note that barycentric velocity correction and local normalisation to the continuum level were applied to all the spectra, achieved through a polynomial function fitting the continuum.

Fig. A.1 All CCF mean profiles of the photospheric CO lines (in red) and their mean profiles (in black).

Appendix B MCMC sampling plot for best-fit QP-GP hyperparameters and RV orbital parameters.

We perform an MCMC sampling to simultaneously fit a combined QP-GP model for the spot and a Keplerian orbit model. Our analysis begins by defining a log-likelihood function which quantifies the agreement between the combined model and observed data through a chi-squared fit. Subsequently, we computed the posterior probability distribution, which integrates information from both the log-prior and log-likelihood functions. To explore parameter space, we employed a total of 256 walkers and executed 10,000 iterations during the MCMC sampling process. The initial positions for these walkers were set to correspond to initial parameter guesses.

We set uniform priors for all the QP-GP hyperparameters and Keplerian orbit parameters with values listed in Table 2. The resulting MCMC chain was recorded, and posterior distributions were generated based on the entire 10,000 iterations with 256 walkers. For improved sampling efficiency and to reduce autocorrelation within the chain, we applied a thinning factor of 10. Consequently, this led to a final distribution of 256,000 points for each sampled parameter (Figure B.1).

Fig. B.1 Corner plot of QP-GP hyperparameters and best-fit Keplerian orbit model parameters to the combined ESPaDOnS and SPIRou RV dataset, after 10,000 MCMC iterations with 256 walkers thinned to 256,000 points per sampled parameter.

References

All Tables

All Figures

	Fig. 1 Light curve, GLS and L1 periodograms for the K2 data. Top: K2 light curve of CI Tau. Middle: GLS periodogram of the K2 light curve, showing three peaks at 6.6, 9.0, and 24.2 d indicated with blue dashed lines. Bottom: L1 periodogram, showing the same periods as found in the GLS periodogram, with an additional significant peak at 47.2 d.
In the text

	Fig. 2 K2 light curve (in black) phased on the 6.6, 9.0, and 24.2 d periods (from top to bottom). The phase-binned mean light curve is shown in red.
In the text

Fig. 3 LCOGT light curves and L1 periogogram for each band, 𝑔 (green), r (blue), and i (red). Top panel: first semester (20B) light curve observed between 31 October 2020 and 08 January 2021. Second from top panel: second semester light curve taken between 27 September 2022 and 31 March 2023. Third from top panel: L1 periodogram of the combined dataset for semester 20B. Bottom panel: L1 periodogram of the combined dataset for semester 22B/23A.

In the text

	Fig. 4 Power spectrum of the photometric time series (black line) fitted with a power law of the form P(f) = Nf−α (red continuous line). The red dotted lines are at the 0.1 and 0.01% FAP levels from bottom to top. respectively. The 25.7 d period is clearly significant above the 0.01% FAP.
In the text

	Fig. 5 Cleaning of stellar activity signals from the K2 light curve. Top panel: quasi-periodic GP model fit (in red) with two periodic components of 6.6 and 9.0 d. The shaded region in the top panel shows 1, 2, and 3σ from the mean model. The bottom panel shows the L1 periodogram of the residual light curve after removing the modulation due to stellar spots.
In the text

	Fig. 6 L1 periodogram of RVs related to 72 data points (four nightly sub-exposures) from ESPaDOnS. The spot period is recovered at 9.27 d and another period is recovered at 6.9 d.
In the text

	Fig. 7 2D periodogram and mean LSD line profiles using 72 ESPaDOnS spectra. Top panel: 2D periodogram with colour range corresponding to the power of the periodogram. The horizontal lines are drawn at 9.0 d and 24.2d. The contours represent FAP levels at 0.01% (black), 0.1 (grey), and 1% (white). Bottom panel: mean LSD profile in black, with one standard deviation of all profiles shaded in blue.
In the text

Fig. 8 Bisector span of the mean LSD profiles from ESPaDOnS spectra. A correlation is seen at lower velocities up to around 17 km s−1 and then no correlation where it reaches a plateau. The data points are colour-coded based on phase of the 9.01 d spot rotation period. The dashed line indicates the threshold below which BIS span values were selected to compute the L1 periodogram shown in red in Fig. 9.

In the text

Fig. 9 LI Periodogram of RVs for the whole ESPaDOnS RV dataset (black), RVs selected based on BIS span values of >0.1 km s−1 (blue), and RVs selected based on BIS span values of <0.1 km s−1 (red). The level of activity due to stellar spots, corresponding to 6.9 and 9.27 d (vertical green-dotted lines) is high in the whole RV dataset and RVs corresponding to BIS span values of >0.1. For the flatter part of the BIS span with values <0.1, the RVs display a clear period peak at ~24.6 d (vertical red dotted line).

In the text

	Fig. 10 L1 periodogram of CCF RVs extracted from the full SPIRou wavelength range (in black) and CCF RVs related to the 2293–2340 nm wavelength range (in red). The dotted vertical lines mark the significant peaks based on their FAP, one of them being the 25.5 d peak.
In the text

	Fig. 11 Observed (red) and modelled (blue) RVs vs. BIS span for two different cases. Case 1 (left panel): Gspot modulated at 9 d and Gstar kept fixed and case 2 (right panel): spot modulated at 9 d and star Doppler shifted with semi-amplitude of 0.3 km s−1 at 25.2 d period. Individual point error bars are also plotted.
In the text

Fig. 12 TBV vs. probability density after 10 000 MC simulations for three cases. In case 1 (orange) the star is at rest and the spot modulated at 9 d; in case 2 (blue) the star is Doppler shifted with a semi-amplitude of 0.3 km s−1, keeping all other parameters the same as in case 1. In case 3 (green) we test a scenario proposed by Donati et al. (2024), i.e. a contribution from material in the inner disk on the CO lines, that could cause the 25.2 d periodicity. The distributions represent the TBVs of the BIS span plots. The grey dotted line is the observed TBV for the SPIRou BIS span vs. RVs plot.

In the text

Fig. 13 RV orbital analysis. Top: L1 periodogram of the full ESPaDOnS and SPIRou RVs. The 25.5 d period peak, the spot modulation period peak at 9.0 d, and the 11.6 d period peak, with their associated FAPs are indicated. Middle panel: best-fit QP-GP combined model (9 d spot + 25.2 d planet) to the RVs. Bottom panel: best-fit Keplerian orbit model (black curve) overplotted on RVs cleaned from the spot GP model. An arbitrary phase space shift of 0.25 has been introduced to enhance the clarity of the orbital variation. The red dots illustrate a binned median of the data points, with bin sizes chosen to maintain uniform point distribution within each bin. The error bars on the red dots represent the standard deviation within each bin.

In the text

Fig. 14 2D periodogram of simulated RVs to show the peak detection limits of an injected planet at 25.2 d period in the presence of spot modulation at 9 d. The value of h/K derived for CI Tau is shown with the white dashed line at 1.1. The cyan line shows the detection limit in the optical where we expect the activity RV signals to be ~ four times larger. The light orange dashed line at 1.7 depicts the h/K expected in the H band.

In the text

	Fig. 15 Optical variability amplitudes and colour slopes of CI Tau, following the classification of Venuti et al. (2015). It suggests that CI Tau’s colour behaviour is consistent with accretion-dominated sources.
In the text

Fig. 16 L1 periodogram analysis of LCOGT light curves and trailed plot of He I line. Top: time series of EW measured for He I, Paβ, Paγ, and Paδ lines. Middle: L1 periodogram of EW. A clear peak is detected in the strongest line tracing accretion (Paβ) at ~27.0 d. Additional peaks that are relatively significant are seen at 7.2 d, 9.1 d (spot), 11.6 d, 18.5 d (see text). Bottom: trailed plot of the He I line phased at a period of 19.5 d showing high-velocity outflows traced by the blue absorption component in the line. The colour bar represents the normalised flux.

In the text

	Fig. A.1 All CCF mean profiles of the photospheric CO lines (in red) and their mean profiles (in black).
In the text

	Fig. B.1 Corner plot of QP-GP hyperparameters and best-fit Keplerian orbit model parameters to the combined ESPaDOnS and SPIRou RV dataset, after 10,000 MCMC iterations with 256 walkers thinned to 256,000 points per sampled parameter.
In the text