© ESO 2020

1 Introduction

The exoplanetary systems known to date show a large variety of architectures that results from the diverse outcomes of planet formation and evolution processes. Planetary migration mechanisms are acknowledged to be the main factor that shaped the observed systems, and that might be the origin of the hot Jupiters (HJs, e.g. Dawson & Johnson 2018). Planet–disc interaction (Baruteau et al. 2014), high-eccentricity migration (Rasio & Ford 1996) produced by secular interaction among bodies in the system or planet–planet scattering, or in situ formation (Batygin et al. 2016) are expected to produce observable trends in the planet population that can be used to gauge their respective effectiveness. Theoretical works partially describe the observed distribution of the HJ population (Ford & Rasio 2008; Matsumura et al. 2010; Hamers et al. 2017), which seems to be mainly produced by a high-eccentricity migration process associated with tidal interactions (e.g. Bonomo et al. 2017). However, this scenario cannot fully explain the observational evidence, and a clear view of the conditions that favour one mechanism over the others is still lacking. Information on the HJ formation path can be obtained by determining their orbital parameters, in particular, eccentricity and obliquity, but also from understanding the migration timescales and age-dependent frequency of different types of systems. However, these clues cannot be easily provided by the available and well-known distribution of mature systems. Instead, observation of HJs around young stars allows directly spotting the ongoing planetary evolution and provides crucial indications to this open question.

In recent years, first detections of exoplanets in young open clusters and stellar associations have been claimed (e.g. Quinn et al. 2012, 2014; Malavolta et al. 2016; Mann et al. 2017, 2018). One of the most intriguing results is the apparent high frequency of HJs around stars younger than a few dozen million years (Donati et al. 2016; Yu et al. 2017, and recently Rizzuto et al. 2020) relative to their older counterparts. This finding places strong constraints on HJ migration timescales, showing that planet-disc interactions may play a significant role in the genesis of such planets. In this respect, the role of dynamical interaction of planets with perturbing stars within a cluster as well as planet-planet interactions in a multiple planetary system cannot be neglected (Cai et al. 2017; Flammini Dotti et al. 2019; van Elteren et al. 2019). Moreover, a short migration timescale would imply that HJs could undergo strong X and UV irradiation from their hosts, sufficient to remove their outer envelope and modify their physical properties with time (e.g. Locci et al. 2019). This knowledge still relies on a small number of discoveries because a robust confirmation of the presence of such young planets is generally difficult. The typically very high levels of activity of the host stars hamper detections, in particular, for blind searches using the radial velocity (RV) method.

Even when evidence of a planetary companion is found, for instance through transits observed in the light curves of the space telescopes Kepler/K2 and Transiting Exoplanet Survey Satellite (TESS), the amplitude of the RV signal generated by the stellar activity could be up to several hundred m s⁻¹, and the planetary signal could go undetected even when sophisticated modelling is used to account for the activity of the host star. All this makes the measurement of planetary mass and bulk density very challenging. Exemplary cases are represented by two of the handful of exoplanetary systems younger than ~20 Myr. The first is represented by the super-Neptune-sized companion to the low-mass star K2-33 (R_p = 0.451 ± 0.033 R_Jup), a M3V star in Upper Scorpius (age 11 Myr; David et al. 2016; Mann et al. 2016). The expected RV semi-amplitude due to the planet is about 20 m s⁻¹, which is dwarfed by an activity variability of about a few hundred m s⁻¹. A reliable mass determination is still lacking for this object, preventing the understanding of the planet bulk structure and further studies of its evolution at early stages based on solid observational results. The second case, still more complicated, is represented by the multi-planetary system V1298 Tau (age ~20 Myr) detected by Kepler/K2 (David et al. 2019a), with a Jupiter-sized planet cohabiting with three more companions, all between the size of Neptune and Saturn. The high-amplitude variability in the RVs caused by stellar activity (~200 m s⁻¹ over nearly five days, as measured from Keck/HIRES VIS RVs by David et al. 2019b), and dynamical effects due to mean-motion resonances, makes the characterisation of the V1298 Tau system very challenging with RV follow-up. More recently, the detection with TESS and Spitzer of a 0.4 R_Jup transiting the bright, pre-main-sequence M dwarf AU Mic every ~8.5 d (Plavchan et al. 2020, age ~20 Myr) made headlines, in that the planet AU Mic b co-exists with a debris disc. The RV follow-up of the star revealed a variability due to stellar activity with amplitudes of ~150 and 80 m s⁻¹ in the visible and near-infrared, respectively, that allowed only for a measurement of the planet mass upper limit (<0.18 M_Jup, or K < 28 m s⁻¹, at 3σ confidence).

In 2012, the Global Architecture of Planetary Systems (GAPS) project (Covino et al. 2013) started a large and diversified RV campaign with the HARPS-N spectrograph (Cosentino et al. 2014) at the Telescopio Nazionale Galileo (TNG) focused on exoplanetary science. One main goal pursued by GAPS is assessing the planet occurrence rates around different types of stars (e.g. Barbato et al. 2019), and understanding the origin of planetary-system diversity. Since 2017, the characterisation of exoplanetary atmospheres (Borsa et al. 2019; Pino et al. 2020; Guilluy et al. 2020) and the RV search for planetary companions around young stars (Carleo et al. 2018, 2020) became main scientific themes. The RV survey was specifically designed to confirm the emerging evidence of a higher frequency of planets around young T Tauri stars than around more evolved stars (e.g. Yu et al. 2017, underlining, however, that the sample is still too small for any reliable statistics), and to determine their orbital and physical parameters for a comparison with the older population.

Within this framework, we monitored a sample of targets in young associations (e.g. Taurus, Cepheus, AB Doradus, Coma Berenices, and Ursa Major) to search for planetary companions, and a small sample of targets with confirmed or candidate planets from other surveys (e.g. Carleo et al. 2020). In these, we observed the weak-line T Tauri star V830 Tau (age ~2 Myr), known to host a HJ (minimum mass m_b sin i_b = 0.57 ± 0.10 M_Jup; P_b = 4.927 ± 0.008 d) announced by Donati et al. (2016), and further characterised by the same team (Donati et al. 2017, hereafter DO17). This detection came as a breakthrough and was followed by the transiting HJ HIP 67522 b (Rizzuto et al. 2020) because it showed that giant planets can form at very early stages of star formation and migrate within a gaseous protoplanetary disc. However, the host star V830 Tau shows a very high level of activity, with an RV scatter of ~660 m s⁻¹ as measured by DO17, which is at least an order of magnitude higher than the semi-amplitude of the detected planetary signal. We therefore became interested in observing this star with a different spectrograph, aiming to confirm the presence of the planet and refine the planetary orbital and physical parameters with a careful treatment of the RV activity signal, whilst also monitoring simultaneously the star with a dedicated photometric follow-up.

In this paper we present the results of our independent 2.5-yr-long follow-up campaign of V830 Tau. It is structured as follows. We describe the original datasets used in our analysis in Sect. 2, and present updated stellar fundamental parameters derived from HARPS-N spectra in Sect. 3. A characterisation of stellar activity observed during our campaign is given in Sect. 4, where we also present the results of a transit search for V830 Tau b in the Kepler/K2 light curve. In Sects. 5 and 6 we discuss several methods and techniques that we used to extract and model the RVs to characterise the planetary signal. We support our conclusions from our search for planet b by quantifying the detection limits through injection-recovery simulations (Sect. 7), and present a summary and final discussion in Sect. 8.

2 Description of the datasets

2.1 HARPS-N spectra

We collected 146 spectra of V830 Tau with the HARPS-N spectrograph (Cosentino et al. 2014) between 17 October 2017 and 15 March 2020 (time span 878 d), almost two years after the observations of DO17, with a median signal-to-noise ratio (S/N) of 19.6 measured over all the échelle orders. We excluded from further analysis the spectra collected at epochs BJD 2 458 052.702 and 2 458 098.589, which have S∕N = 4.1 and 0.5, respectively. HARPS-N is a cross-dispersed high-resolution (R = 115 000) and high-stability échelle spectrograph, covering the wavelength range 3830–6930 Å. The spectrograph is fed by two fibres, one on the target and the second, used as a reference, illuminated by the sky in the case of V830 Tau.

Following the method described in Malavolta et al. (2017), we did not find evidence for spectra contaminated by moonlight, which in general could affect the RV of the target in a measurable way. During the last season, from 16 November 2019 to 28 February 2020, we adopted a denser sampling by scheduling the target twice per night whenever possible, with the pair of observations separated by at least three hours. This change in the observing strategy was intended to improve the fit of the short-period component of the activity signal, modulated over the known ~2.7 d rotation period, especially during consecutive observation nights.

2.2 Photometric light curves

V830 Tau was observed by the Kepler extended K2 mission in long-cadence mode (one point every 30 min) during campaign 13, from 8 March to 27 May 2017 (K2 target ID EPIC 247822311). Almost 1.5 yr after K2 observations, we followed V830 Tau up from October 2018 to the end of January 2020 with the STELLA facility in Tenerife (Strassmeier et al. 2004)and its wide-field imager WiFSIP. This time span corresponds to the last two seasons of our monitoring with HARPS-N. Blocks of five exposures per filter were collected, with exposure times per single image of 60 s (V band) and 25 s (I band). Standard data reduction, including bias subtraction and flat field correction, was performed, and aperture photometry was used to extract the differential light curve, as described in Mallonn et al. (2018). We then analysed the averaged values per observing block and filter, resulting in 125 and 122 data points for the V and I filter, respectively. The ensemble of comparison stars (UCAC4 573-011610, 573-011607, 574-011176, and 573-011630) was chosen in automatic mode by the pipeline to minimise the scatter of the differential photometry.However, due to the strong variability of the target, a specific choice of the reference stars has minor effects on the differential light curve. The data are listed in Tables A.1 and A.2.

3 Stellar parameters and lithium abundance

Spectroscopic determination of stellar parameters through standard techniques based on line equivalent widths (EWs) and imposing excitation or ionisation equilibria is not possible for our star because of its low effective temperature (T_eff) and relatively rapid rotation. We therefore performed a spectral synthesis analysis of the co-added spectrum of the target in four spectral regions around λλ5400, 5800, 6200, and 6700 Å, and using the 2017 version of the MOOG code (Sneden 1973; 2017 version) with the driver synth. We considered the line list kindly provided by Chris Sneden (priv. comm.), the Castelli & Kurucz (2004) grid of model atmospheres, with solar-scaled chemical composition and new opacities (ODFNEW), and linear limb-darkening coefficient taken from Claret (2019). Assuming [Fe/H] = 0.0 dex as iron abundance, in compliance with the metallicity of the Taurus-Auriga star-forming region as derived by D’Orazi et al. (2011), microturbulence velocity of 1.0 km s⁻¹ and macroturbulence of 1.4 km s⁻¹ by Brewer et al. (2016), we varied the effective temperature in the range of 3800–4500 K and surface gravity in the range 3.6–4.5 dex at steps of 100 K and 0.1 dex, respectively, until the best match between target and syntheticspectra was found for each spectral region. In the end, the mean values of the effective temperature obtained from the four spectral regions were adopted. The final values of the derived parameters were T_eff = 4050 ± 100 K and log g = 3.95 ± 0.20 dex, where the errors take into account both the best-fit determination and the standard deviation of the mean obtained from the four spectral regions considered. Our estimate for T_eff is slightly lower than that reported by DO17, who adopted T_eff = 4250 ± 50 K from Donati et al. (2015), but consistent with values found by other authors (see e.g. Sestito et al. 2008, and references therein). As a byproduct, we also derived a projected rotational velocity of v sin i = 30 ± 1 km s⁻¹, consistent with the value 32.0 ± 1.5 km s⁻¹ reported by Nguyen et al. (2012).

After the stellar parameters were derived, we also measured the EW for the lithium doublet at 6708 Å, finding EW(Li I) = 658 ± 5 mÅ, which corresponds to an abundance A(Li I) = 3.19 ± 0.16 dex (Lind et al. 2009). The error on T_eff is the dominant source of uncertainty, implying an error of 0.15 dex in A(Li I). This means that our measured Li abundance agrees within the uncertainties with values measured in T Tauri stars and meteoritic abundance. Furthermore, our value is consistent within the errors with previous findings using the same method (see Sestito et al. 2008, and references therein). Applying non-local thermal equilibrium (NLTE) corrections and following the prescription by Lind et al. (2009), we obtain A(Li I)_NLTE = 3.16 dex.

Fig. 1 Kepler/K2 light curve. The time series collected during the first half of 2017, and the same data phase-folded to the stellar rotation period are shown in the upper and lower panel, respectively. The most powerful flares occurred at phases close to maximum brightness.

4 Stellar activity analysis

4.1 Light-curve analysis

Figure 1 shows the time series of Kepler/K2 photometry and the data phase-folded to the known rotation period P_rot = 2.74 d. We corrected the K2 light curve following the procedure described in detail in Nardiello et al. (2016) and Libralato et al. (2016). Briefly, we performed a two-step correction: first, we corrected the systematic trends associated with the conditions of the spacecraft, of the detector, and the environment. In order to perform this correction, we fitted and applied to the light curve of the target a linear combination of orthonormal bases, the co-trending basis vectors, released by the Kepler team¹. In a second step, we corrected the position-dependent systematics due to the large jitter of the spacecraft. We refer to Sect. 3.1 of Nardiello et al. (2016) for a detailed description of the corrections. We flattened the light curve as done in Nardiello et al. (2019, 2020): we defined anumber of knots spaced 6.5 h on the light curve, and interpolated these knots with a fifth-order spline to obtain a model of the stellar variability. We used this model to flatten the light curve. In this process, we also excluded all the points 3σ above the average value of the flattened light curve. The corrected photometric data, phase-folded at the orbital period P = 4.927 d found by DO17, are shown in Fig. 2. We did not detect any transit signal at the period of the claimed planet, as further confirmed by the analysis with the box least-squares algorithm (BLS; Kovács et al. 2002). For comparison, assuming R_⋆ = 2 R_⊙ (after DO17), and a planetary radius between 1.5 and 2 R_Jup, estimated from the giant planet thermal evolution models by Fortney et al. (2007) for a planet of mass, age, and orbital distance (scaled to take into account the irradiation of an M star, according to its luminosity) as those found by DO17, we would expect a transit with depth between ~0.6 and 1%, that is, 6–10 times the RMS of the flattened light curve.

Figure 3 shows the differential photometry of STELLA and their generalized Lomb-Scargle periodograms (GLS; Zechmeister & Kürster 2009), which show clear peaks at the rotation period and 1d aliases. We note that the full amplitude of the STELLA V band photometry is the same (ΔV = 0.28 mag) as that measured by DO17 between 2015 Oct 30 and 2016 Mar 15.

Fig. 2 Detrended and flattened K2 light curve of V830 Tau, phase-folded at the orbital period of the planet announced by DO17. Red dots represent five-point binned data.

Fig. 3 STELLA light curves in V and I bands, spanning the last two seasons of HARPS-N observations, and their respective GLS periodograms. The peak corresponding to the rotation period is labelled.

4.2 Spectroscopic diagnostics

We extracted the chromospheric activity indexes based on the CaII H&K and Hα spectral lines using the method described in Gomes da Silva et al. (2011) and the code ACTIN v1.2.2 (Gomes da Silva et al. 2018). Their time series and GLS periodograms are shown in Fig. 4. Peaks at the rotation period and 1d aliases are clearly visible for both indicators, which show an increase in the level of the stellar magnetic activity over the time span of our spectroscopic follow-up.

As we show in Sect. 5, the cross-correlation function (CCF) of the HARPS-N spectra appears highly distorted by the high level of stellar activity. Therefore, activity indicators based on a measure of the CCF asymmetry, such as the full width at half maximum (FWHM) and bisector inverse slope (BIS) might not be suitable to correct for the activity signal in the RVs. This isevident from Fig. 5, showing the correlation between the BIS index and one of the RV datasets considered in this study, and where seven BIS outliers are visible. With this caveat, some other indicators of the CCF shape were further considered in our RV modelling attempts, as discussed in Sect. 6.5.

We also analysed the stellar activity by phasing the CCF residuals to the rotational period, that is, the CCFs divided by the average CCF calculated over the whole dataset (Fig. 6). While the shape of the map qualitatively confirms the rotational period of the star of 2.741 d, it is interesting to note that active regions (positive and negative deviations from the average CCF) are moderately stable in position in the three observation seasons. The time series of the spectroscopic activity diagnostics are listed in Table B.1.

Fig. 4 Time series of the activity indicators extracted from the HARPS-N spectra based on the CaII H&K and H-α spectral lines (left panels), and their corresponding GLS periodograms (right panels). The peak corresponding to the rotation period is labelled.

Fig. 5 Correlation plot between the CCF BIS index and the RVs extracted with the DRS pipeline (see Sect. 5). Left panel: all data. Right panel: data without outliers to better visualise the RV-BIS correlation.

5 Extraction of the radial velocities with different methods

To complete the very challenging task of detecting a signal induced by the Keplerian motion of V830 Tau b, whose expected semi-amplitude is more than an order of magnitude smaller than the RV scatter due to magnetic activity, we extracted RVs using three independent methods, and different from the least-squares deconvolution (LSD) used by DO17. Each extraction method could be affected by the very high stellar activity in a different way, therefore we decided to analyse different datasets in the attempt to detect and characterise the signal induced by V830 Tau b.

We used the standard DRS pipeline version 3.7.0 to extract the RVs through the CCF technique (Pepe et al. 2002). To calculate the CCF, we adopted the template mask for a K5V dwarf and a half-window of the CCF of 200 km s⁻¹, to account for the line broadening through the fast stellar rotation and to include a good portion of the continuum for a proper fitting of the CCF profile. Figure 7 shows the CCF of one HARPS-N spectrum with average S/N, which clearly illustrates the strong deformation in the core of the average line profile due to stellar activity.

We also used the template-matching TERRA pipeline (Anglada-Escudé & Butler 2012) to extract an independent dataset. Our default TERRA dataset is that obtained by considering orders corresponding to a reference wavelength higher than λ = 4530Å, as recommended for stars with a spectral type like that of V830 Tau. The computation of the RVs includes a correction for the secular perspective acceleration.

Finally, we extracted the RVs from the observed spectra using the Gaussian-process-based, template-free approach recently proposed byRajpaul et al. (2020) (hereafter identified as dataset R20). In brief, Gaussian processes (GP) are used to model and align all pairs of spectra with each other; the pairwise RVs thus obtained are combined to produce accurate differential stellar RVs, without having to construct a template. Such differential RVs can be extracted on a localised basis, for example to yield an independentset of RVs for each échelle order, or indeed for much smaller subdivisions of orders. The rationale behind the latter approach is that regions of spectra affected by, for instance, stellar activity or telluric contribution may in principle be identified and excluded (effectively a data-driven masking of the spectrum, without any knowledge of line locations or properties) from the calculation of the final RVs, which are obtained by an inverse variance-weighted average of the localised RVs. The RVs used in this work were obtained by combining RVs from each échelle order, and without anymasking. This approach was found to yield the highest signal- to- (estimated) noise ratio. We did, however, explore several alternative schemes for RV extraction, where we divided each order into anything from 16 to 128 chunks, each of which might have contained between zero and several lines, and then selectively recombined these localised RVs trying to produce a final set of RVs that minimised correlations with the FWHM or BIS time series, minimised periodogram power near the stellar rotation period of 2.74 d, and/or maximised the power near the putative planetary orbital period of 4.93 d, for example. We explored both iterative optimisation schemes and more sophisticated machine-learning approaches (e.g. the HDBSCAN algorithm; Campello et al. 2013) to optimise the masking procedure. Unfortunately, we found that virtually all the useful Doppler information was contained in spectral regions that are strongly contaminated by rotational activity: all attempts to suppress this stellar signal while trying to boost periodicity at 4.93 d led to RV error bars that were at least an order of magnitude larger (>250 m s⁻¹) than in the mask-free case, thus thwarting our attempts to tease out the putative planetary signal.

We list in Table C.1 the datasets used in this work, and show the time series in Fig. 8. We note that the third season is characterised by a higher RV dispersion, likely due to the increasing stellar activity visible in the spectroscopic diagnostics, as discussed in Sect. 4. We summarise in Table 1 the main properties of each dataset compared to that of the MaTYSSE large programme RVs analysed by DO17, which were collected over 91 days with the ESPaDOnS and Narval spectropolarimeters linked to the 3.6 m Canada-France-Hawaii, the 2 m Bernard Lyot, and the 8 m Gemini-North telescopes. ESPaDOnS and Narval collect spectra covering the wavelength range 3700–10000 Å, which overlaps with that of HARPS-N and extends to the near-IR (NIR) region, with a resolving power of 65 000.

We note that the median internal error σ_RV of the HARPS-N RVs (TERRA and DRS extraction) is nearly half that of the MaTYSSE data, and that the TERRA dataset has a scatter reduced by ~28% and 9% with respect to that of the DRS and R20 extractions, respectively.

Fig. 6 Residuals of the CCF, i.e. each individual CCF is divided by the average CCF calculated over the whole dataset, phase-foldedto the rotational period of V830 Tau (2.741 d).

Fig. 7 Cross-correlation function of V830 Tau for a HARPS-N spectrum with average S/N. The CCF was calculated by the DRS pipeline using a template mask for a KV5 star.

6 Radial velocity analysis

We describe below the results obtained from the analysis of our HARPS-N RVs. We start with showing the GLS periodograms, to illustrate that the time series are clearly dominated by signals produced by stellar activity and ± 1 d⁻¹ aliases.

6.1 Frequency content analysis

We calculated the GLS periodograms for the original data and residuals after recursive pre-whitening, as shown in Fig. 9 for all the datasets. The periodogram of the original data (first panel) shows a very sharp maximum at the stellar rotation frequency, and signals related to stellar activity (at the rotation frequency or its harmonics, and 1d aliases) clearly dominate the periodograms even after five pre-whitening iterations. The periodogram of the RV residuals after the last pre-whitening is still characterised by high dispersion (228 and 327 m s⁻¹ for TERRA and DRS data, respectively), around ~3–5 times the semi-amplitude of the claimed signal induced by the planet. In general, pre-whitening using only sinusoids is not an optimal wayto account for stellar activity and search for planetary signals, and more sophisticated functions should be used to model the complex activity-related signals; however, these results do at least illustrate well that unearthing a small planetary signal is not a trivial task for V830 Tau. Taking advantage of the fact that the STELLA data span the last two seasons of the RV follow-up, it is interesting to compare the photometric and spectroscopic datasets by phase-folding the data to a common period and phase-zero epoch. This can provide some insights into the nature of stellar surface patterns responsible for the periodic modulation observed in the RVs. We present this comparison in Fig. 10, using the TERRA dataset. The light curves and the spectroscopic activity indicators are anti-correlated, and this can be explained by spot-dominated activity. This evidence is also confirmed by the smaller amplitude of the I band light curve. This indicates that active regions are dominated by cooler features (star spots) rather than hotter, facular-like features. The ~0.2 phase shift between RVs and light curves is also typical of the effect related to the flux deficit due to spots that affect the CCF. We note some differences in the RV-folded curves distinguished by observing season.

Fig. 8 Radial velocities of V830 Tau extracted from the HARPS-N spectra with different methods (average subtracted). From top to bottom: DRS, TERRA, and template-free algorithm (R20) by Rajpaul et al. (2020). The error bars are nearly two orders of magnitude smaller than the major ticks on the y-axis, and are not visible. For each season we indicate the RMS of the corresponding RV subsample.

6.2 Gaussian process regression

We turned to more sophisticated analysis techniques to mitigate the stellar activity contribution to the variability observed in the RVs. Gaussian process regression, which has often been applied to detect and characterise planetary signals in RV time series (e.g. Haywood et al. 2014; Dumusque et al. 2017), was used by DO17 to derive their planet parameters from the raw ESPaDOnS, Narval, and GRACES RVs, and proved to be a very efficient way to model the activity over the shorter time span of their observations (about three months). We apply here the same technique and model, using a quasi-periodic covariance matrix for the correlated signal due to stellar activity, $$\begin{array}{lll}k(t,t^{\prime })\hfill & =\hfill & h^2\cdot \exp \left[-\frac{{(t-t^{\prime })}^2}{2{\lambda }^2}-\frac{{\text{sin}}^2(\pi (t-t^{\prime })/\theta )}{2w^2}\right]\hfill \\ \hfill & \hfill & +({\sigma }_{\text{RV}}^2(t)+{\sigma }_{\text{jit}}^2)\cdot {\delta }_{t,t^{\prime }},\hfill \end{array}$$(1)

where t and t′ represent two different epochs, σ_RV is the radial velocity uncertainty, and ${\delta }_{t,t^{\prime }}$ is the Kronecker delta. Our analysis takes into account other sources of uncorrelated noise (instrumental and/or astrophysical) by including a constant jitter term σ_jit that is added in quadrature to the formal uncertainties σ_RV. h, λ, θ, and w are the GP hyper-parameters: θ represents the periodic timescale of the modelled signal, and corresponds to the stellar rotation period; h denotes the scale amplitude of the correlated signal; w describes the level of high-frequency variation within a complete stellar rotation; and λ represents the decay timescale in days of the correlations, which relates to the temporal evolution of the magnetically active regions responsible for the correlated signal observed in the RVs.

We performed a Monte Carlo analysis with the open-source Bayesian inference tool MULTINEST v3.10 (e.g. Feroz et al. 2013), through the PYMULTINEST PYTHON wrapper (Buchner et al. 2014), including the publicly available GP PYTHON module GEORGE v0.2.1 (Ambikasaran et al. 2015). Our setup included 500 live points and a sampling efficiency of 0.3. The use of a nested sampler allows for a robust Bayesian model comparison through the calculation of marginal likelihood (or evidence) $\mathcal{Z}$ for each model with good accuracy, which is a crucial point for our analysis. We then compare the Bayesian evidence of a model containing only the correlated stellar activity term (that includes 5 free (hyper-)parameters), with that of a model including a planetary signal (that includes 9 or 11 free (hyper-)parameters, for a circular and eccentric orbit respectively) for a robust statistical analysis of our dataset.

Because our goal is an independent confirmation of the presence in our data of the planetary signal detected by DO17, we adopted uninformative priors for all the free parameters except for θ to guaranteean unbiased analysis. The stellar rotation period can be reliably constrained using a Gaussian prior based on the result of a GP quasi-periodic regression of the H-α activity diagnostic time series (θ = 2.7417 ± 0.0007 d), which is extracted from the same spectra used to derive the RVs. However, we adopted a more conservative value for the σ of the prior, that is, one order of magnitude larger that the uncertainty associated with the rotation period derived from the H-α time series. The orbital period of the planet was uniformly sampled up to ten days.

The results of the analysis for each of the different RV datasets are summarised in Table 2. As an example of posterior distributions of model parameters, we show the corner plot for the GP+1 planet model (TERRA RVs) in Fig. D.1. The planetary signal detected by DO17 is not recovered in any of our datasets, and the marginal likelihoods always favor the zero-planet model. The GP regression is able to model the stellar activity signal effectively, as can be seen by comparing the RMS of the original data to the RMS of the residuals. However, the latter are still above 100 m s⁻¹, which is nearly three times larger than the RMS of the residuals of DO17 (35 m s⁻¹). We do not have an explanation for this observed difference, which may be partly due to a higher level of activity of the star during our follow-up, or could be partly explained with the different wavelength ranges covered by HARPS and ESPaDOnS/Narval, with the latter reaching the NIR region where the RVs are expected to be less contaminated by stellar activity.

We performed one more test by taking the RVs extracted with TERRA using échelle orders starting from no. 45, with the central wavelength λ = 5463Å, that is, we used a narrower region corresponding to a redder part of the spectra. The RVs have median uncertainties σ_RV = 10.2 m s⁻¹ and RMS of 828 m s⁻¹, which are slightly smaller than that of the default dataset, in agreement with the evidence from STELLA data that the photometric variability is lower in I band than in V band (Fig. 10). Despite the lower scatter due to a reduced contribution from stellar activity, we did not find evidence for the planetary sign.

Leaving the eccentricity unconstrained does not improve the fit (e.g. we obtain $e_{\text{b}}={0.46}_{-0.31}^{+0.36}$ and $\ln \mathcal{Z}=-992.3$ for the TERRA dataset). We also used a looser prior for P_b, increasing the upper limit to 100 d, in order to explore the possible presence of longer period planets. Even so, we did not find evidence for any significant signal in the data.

Fig. 9 GLS periodograms of the original TERRA (a), DRS (b), and R20 (c) RVs and their residuals after recursive pre-whithening. For each figure: the vertical and dashed blue lines indicate the stellar rotation frequency and its first and second harmonic; the red line marks the orbital frequency of the planet announced by Donati et al. (2016, 2017); each panel reports the RMS of the dataset used for calculating a specific periodogram; the window function is shown in green in the bottom panel.

6.3 Gaussian process RV modelling jointly with ancillary activity indicators

To investigate the interplay between RV variations and stellar activity in more detail, we applied a more sophisticated GP-based approach. We analysed all the different RV datasets in Table 1 using the framework described by Rajpaul et al. (2015). The RVs were fitted jointly with the DRS CCF asymmetry indicators BIS and FWHM. In short, this GP framework assumes that all observed stellar activity signals are generated by some underlying latent function G(t) and its derivatives; this function, which is not observed directly, is modelled with a GP with a quasi-periodic covariance function. G(t) and its derivative are then allowed to manifest (using physically motivated relationships) in all observable, activity-sensitive time series, while Keplerian terms for one or more planets are incorporated into the RVs alone. This GP-based approach to model RVs jointly with activity indicators might enable a more reliable planet characterisation than traditional approaches that assume simple parametric forms for the stellar signals, or that try to exploit simple correlations between RVs and activity indicators.

In an effort to detect V830 Tau b in our RVs, we combined the GP-based activity framework with models including either one or no planet(s). We placed uninformative priors on all standard parameters in the Rajpaul et al. (2015) framework, as well as on the planet parameters, except for the period, which we constrained to 4.93 ± 0.25 d (i.e. we now searched for an expected signal at this period, rather than searching blindly in all the periods). To compute model posteriors and Bayesian evidence, we employed POLYCHORD (Handley et al. 2015), which is a state-of-the-art nested sampling algorithm, and an efficient alternative to MULTINEST, designed to work especially with very high dimensional parameter spaces.

In brief, we found that the Bayes factor for the one-planet models versus the zero-planet models ranged from1.26 to 5.52, depending on the RV extraction algorithm used (e.g. DRS vs. TERRA): in no instance, then, was a one-planet model strongly favoured. More decisively, the RV semi-amplitude for the 4.93 d period Keplerian was in all cases consistent with zero within 1σ, indicating the non-detection of V830 Tau b.

In various other tests where we used the same GP framework but replaced the narrow planet period prior with an uninformative one, one-planet models were always rejected outright compared to the zero-planet case. The period posteriors had a low probability density around ~ 4.9 d, and the semi-amplitudes associated with ~4.9 d period samples were clustered tightly around zero. These results indicate even more strongly than in the case of the narrow period prior a non-detection of V830 Tau b.

Fig. 10 STELLA light curves (upper panel), CaII H&K and H-α activity indicators (middle panel), and TERRA RVs (lower panel) phase folded using the same period P = 2.741 d and reference epoch of the STELLA V band photometry andKepler/K2 light curve in Fig. 1. Subplots for the light curves and RVs show the residuals O–C of the best-fit sine function obtained with the GLS software.

6.4 Modelling the RV variations from the observed photometric modulation

Wide-band photometry can be used to map the longitudinal distribution of active regions on the surface of an active star and to predict the activity-induced RV variations to some extent (e.g. Lanza et al. 2011). In Fig. 11 we plot the two seasons of V band optical photometry of V830 Tau versus the rotation phase together with continuous interpolations obtained for the individual seasons as well as for the dataset as a whole. The rotation phase is computed assuming a constant rotation period of P_rot = 2.7409 d. To computethe interpolations, we performed a kernel regression (KR), that is, a locally linear regression of the RV versus the phase giving decreasing weights to the data points that are more distant in phase from the given data point for which the regression value is to be computed (see Sect. 6.5 for details). Only minor changes occur between the two seasonal light curves, although they are separated by more than 300 d. This indicates that the photospheric brightness inhomogeneities are very stable in V830 Tau. Therefore, we consider our V band photometricdataset as a whole, thus obtaining a more continuous phase coverage for our subsequent analysis. An analogous approach can be applied to the I band light curves, but we focus on the V band light curves because they show a flux modulation of greater amplitude because the star spot contrast is higher in the optical, which permits a more precise RV reconstruction.

To compute the activity-induced RV variations, we applied the so-called FF′ method introduced by Aigrain et al. (2012). It accounts for both the variation ΔRV_rot induced by the spectral line distortions produced by surface brightness inhomogeneities, the visibility of which is modulated by stellar rotation, and the variation ΔRV_c produced bythe inhibition of surface convection in the regions where photospheric magnetic fields are more intense, which reduces the local convective blueshifts of spectral lines. We adopted the following expressions for the two components: $$\Delta {\text{RV}}_{\text{rot}}(\varphi )=\frac{A}{F_0}\frac{\text{d}F(\varphi )}{\text{d}\varphi }\left[\frac{F(\varphi )}{F_0}-1\right]$$(2)

and $$\Delta {\text{RV}}_{\text{c}}(\varphi )=B{\left[\frac{F(\varphi )}{F_0}-1\right]}^2,$$(3)

where ϕ is the rotation phase, F(ϕ) the interpolated V band flux at phase ϕ, F₀ the flux in the absence of spots (that we take equal to the maximum flux along the interpolated light curve), and A and B are two coefficients to be determined together with the RV offset RV₀ between the model and the observations by minimising the χ². This is computed as the sum of the squares of the residuals between the model RVs and the observations, normalised by the respective standard deviations of the RV measurements.

The χ² minimisation with respect to the RV extracted with the TERRA procedure gives the model plotted as a solid red line in the top panel of Fig. 12, the residuals of which are plotted in the bottom panel of the same figure and have a standard deviation of 509.33 m s⁻¹, while the original RV time series in the top panel has a standard deviation of 874.99 m s⁻¹. The model RV variations are dominated by the effect of the brightness inhomogeneities with a mean value of |ΔRV_rot| equal to 4.93 times the mean value of |ΔRV_c|, as might beexpected given the rather large rotational broadening of the spectral lines (v sini ~ 30 km s⁻¹) and the large area occupied by star spots on V830 Tau. Similar results are obtained with the RV extracted by the HARPS-N pipeline DRS, although with a larger RV dispersion and greater residuals for the FF′ model.

In conclusion, our RV model based on wide-band photometry does not adequately reproduce the observed RV variations of V830 Tau, probably because the pattern of surface brightness inhomogeneities is much more complex than the simple spot distribution assumed by the FF′ model. This is clearly indicated by the Doppler imaging maps of DO17, who cautioned about the limitations of any reconstruction of the RV variations based solely on the photometry for this very active star.

Fig. 11 V band optical light curves of V830 Tau with the flux plotted vs. rotation phase. The data points collected between BJD 58 393.742 and 58 607.363 are plotted as open green diamonds and their kernel regression vs. the rotational phase is given by the dashed green line, while the data points collected between BJD 58 765.645 and 58 916.387 are plotted as open orange triangles and their regression is given by the dot-dashed orange line. The solid red line is the regression to the whole V band photometric dataset vs. the rotation phase.

Fig. 12 Top panel: RVs obtained with the TERRA procedure vs. time (filled dots). The model RV variation as derived from the solid red interpolation in Fig. 11 by means of the FF′ method is superposed (solid red line). Bottom panel: residuals between the observed TERRA RV and the RV predicted by the FF′ model vs. time. In both panels the size of the RV error bars is comparable with that of the data points in most of the cases.

6.5 Kernel regression analysis of the RV time series

In another, complementary analysis of the V830 Tau RVs, we first tried to remove the rotational modulation produced by stellar activity. The evolution timescale of the surface features produced by magnetic activity is comparable with the time span of the observations in individual seasons and is thus an effective method to reduce the activity-induced RV modulation is to place the RV data in phase for each season and make a regression versus the rotation phase. By subtracting this regression, we remarkably reduced the activity-induced RV variations. We considered our three seasons of RV data, the first between BJD 58 044.6623 and 58 192.3573, the second between BJD 58 341.7344 and 58 566.3847, and the third between BJD 58 804.4827 and 58 924.4063, and placed each of them in phase assuming a constant rotation period P_rot = 2.7409 d.

To compute the regression of the RV versus the rotation phase, we performed a KR. More precisely, to compute the regression value for the ith data point observed at time t_i, corresponding to rotation phase ϕ_i, we performed a linear regression of the RV versus the phase over all the data points, giving them a weight $$w_j\propto \exp \left\{-\left[{\left(\frac{{\varphi }_i-{\varphi }_j}{h_{\varphi }}\right)}^2+{\left(\frac{t_i-t_j}{h_{\text{t}}}\right)}^2\right]\right\},$$(4)

where ϕ_j is the rotation phase of the generic jth data point observed at time t_j, while h_ϕ and h_t are the so-called bandwidths; they govern the decrease of the weight w_j of the generic jth data point as it becomes increasingly distant from the considered ith data point. More details on the KR implementation can be found in Lanza et al. (2018, 2019), and references therein.

The results of the application of KR to our seasonal RV datasets extracted with the TERRA procedure are illustrated in Fig. 13. The whole TERRA RV time series consists of 144 data points with a standard deviation of 874.99 m s⁻¹ and is plotted in the top panel with different colours indicating data collected in different seasons. The same colour code is used to plot the corresponding seasonal KRs. The residual time series obtained by subtracting the seasonal KRs has a standard deviation of 123.91 m s⁻¹ and is plotted in the bottom panel. We refer to this residual RV time series as the cleaned RV time series. The mean bandwidths over the three seasons are h_ϕ = 0.095 and h_t = 19.73 d.

To further reduce the RV scatter, we considered the indicators of the shape of the CCF and the chromospheric index $\log R_{\text{HK}}^{\prime }$ that measures the excess flux in the core of the Ca II H&K lines produced by the non-radiative heating controlled by magnetic activity. In addition to the commonly used BIS index, our suite of CCF shape indicators included the contrast of the CCF, its FWHM, ΔV, and V_{asy (mod)} introduced in Lanza et al. (2018). We performed KRs of the cleaned RV time series versus each of these indicators and the time. Additionally, we performed a further KR with respect to the rotational phase and time. In all the cases, as in Lanza et al. (2018), a 3σ clipping wasapplied by performing a preliminary KR to exclude possible outliers.

None of these KRs gave a significant reduction of the standard deviation of the data points as measured by the Fisher-Snedecor F statistics (see Lanza et al. 2019); therefore we simply considered the one giving the smallest standard deviation of the residuals. This was the KR with respect to stellar rotation phase and time, probably because the CCF indicators lose most of their power when the CCF is strongly distorted, as in the case of a very active star such as V830 Tau, while the chromospheric index $\log R_{\text{HK}}^{\prime }$ is not strongly correlated with the photospheric activity that is mainly responsible for the RV variations. This is supported by the analysis of Hα emission performed by DO17. The second KR with respect to phase and time is different from the first KR applied to obtain the cleaned RV time series because it has a longer time bandwidth h_t = 68.15 d, although the phase bandwidth is the same h_ϕ = 0.095. The standard deviation of the residuals after this second KR is 65.33 m s⁻¹ for a total of 140 data points because the 3σ clipping excluded four outliers. The KR applied to the cleaned time series and the obtained residuals are shown in Fig. 14.

In Fig. 15 we plot a GLS periodogram of the residual RV time series in the bottom panel of Fig. 14. The false-alarm probability corresponding to the highest peak is 0.642 as given by the analytical formula of Zechmeister & Kürster (2009); therefore, there is no indication of significant periodicities in the explored period range. The peak closest to the orbital period of DO17 falls at 4.9545 d. By fitting a sinusoid with this period to the RV time series, we find a semi-amplitude of only 18.63 m s⁻¹, much smaller than the orbital RV semi-amplitude of 68 ± 11 m s⁻¹ reported by DO17.

The possibility that our two successive KRs with respect to stellar rotation phase and time might have removed a signal at the period of the putative planet appears to be very low because the two h_t of their time bandwidth are significantly longer than the period of 4.927 d. We acknowledge that the approach we used for our KR analysis is not completely appropriate from a statistical point of view because the activity and the sinusoidal fit should be performed simultaneously rather than applying the GLS to the KR residuals (see e.g. Anglada-Escudé & Tuomi 2015). Nevertheless, it is much simpler and can be adequate for an exploratory analysis such as that presented here. We therefore conclude that even with the alternative and complementary KR technique we cannot detect a significant signal at the period of V830 Tau b.

Fig. 13 Top panel: RV time series of V830 Tau as extracted with the TERRA procedure vs. the rotation phase (filled dots). Different colours indicate data collected in different seasons: black dots indicate data points collected between BJD 58 044.6623 and 58 192.3573 (first season), red dots between BJD 58 341.7344 and 58 566.3847 (second season), and green dots between 58 804.4827 and 58 924.4063 (third season). The KR performed vs. the rotation phase and time in the first season is indicated by the solid black line, in the second season by the solid red line, and in the third season by the solid green line. The RV error bars are smaller than the size of the plotted dots. Bottom panel: RV residuals obtained by subtracting the seasonal KRs from the corresponding data points.

Fig. 14 Top panel: cleaned RV time series of V830 Tau vs. time as obtained after the first KR with respect tothe rotation phase and time (filled dots). The solid red line indicates a further KR with respect to the phase and time that gives the largest reduction of the standard deviation of the residuals in comparison with the KRs computed with respect to the CCF indicators or the chromospheric index $\log R_{\text{HK}}^{\prime }$. The dashed green lines indicate the 1σ confidence range of the KR (see Lanza et al. 2018, for details). Bottom panel: residualsof the KR in the top panel vs. time, with a standard deviation of 65.13 m s−1.

Fig. 15 Generalized Lomb-Scargle periodogram of the residual RV time series of V830 Tau in the bottom panel of Fig. 14. The power is normalised to its maximum in the considered period interval and is plotted vs. the period itself.The vertical dashed green line marks the rotation period, and the dotted red line indicates the orbital period of the planet proposed by Donati et al. (2017).

7 Planet detection through injection-retrieval simulations

The lack of the exoplanetary signal reported by DO17 in our time series needs to be investigated in terms of effective detectability in presence of such a high level of stellar activity. For this purpose, we devised GP-based simulations to test the feasibility of retrieving the planetary signal claimed by DO17 after this was injected into our data (adopting the TERRA dataset), following two different approaches.

Fig. 16 Results for the simulations with direct injection of the planetary signal into the real data.

7.1 Detection sensitivity by direct injection of the planetary signal into the data

The first set of simulations was built by direct injection of the planetary signal into our data after we randomly drew the parameters from normal distributions determined from the DO17 results (K_b = 68 ± 11 m s⁻¹, P_b = 4.927 ± 0.008 d, T _{conj, b} = 2457360.523 ±0.124 BJD). We produced 50 mock datasets, which were analysed with a GP regression including a sinusoid to fit the planetary signal, and with the same set-up used for fitting the real data (e.g. adopting a uniform prior $\mathcal{U}$(0,10) days for P_b). As a figure of merit, we then inspected the distribution of the P_{b, retrieved}/P_{b, inj} ratio between the 50th percentile (P_{b, retrieved}) of each P_b posterior and the corresponding injected orbital period P_{b, injected}. We derived a similar distribution using the maximum a posteriori probability (MAP) values P_{b, MAP} in place of P_{b, retrieved}. Both histograms are shown in Fig. 16. For more than 50% of the simulated datasets we retrieved an accurate value for P_b within the range 4.9 < P_b < 5.0 d, with a median significance level of 9.5σ (20% of the whole dataset having a significance higher than 100σ). For all the datasets of this sub-sample, P_{b, MAP} corresponds to P_{b, inj} with high accuracy, except in one dataset for which P_{b, MAP} = 3.24 d. In this sub-sample, we were not able to recover the semi-amplitude of the planetary signal with the same degree of accuracy. We obtained 0.98 and 0.28 for the median and RMS of the K_{b, retrieved}/K_{b, inj} ratio, respectively. When we assume the MAP values as an estimate for the orbital period, the percentage of the cases for which we can claim an accurate recovered P_b increases to ~60%. Within this sub-sample, we found five datasets for which P_{b, retrieved} is not in the range 4.9–5.0 d, while the median K_{b, retrieved} is overestimated by 20%.

We tested our ability to recover the planetary signal also by setting the semi-amplitude to a couple of illustrative higher values K_b = 100 and 130 m s⁻¹, and increasing the upper prior bound to 200 m s⁻¹. For both cases, we considered only one realisation. For the first case, we retrieved K_b =${89}_{-36}^{+18}$ m s⁻¹ (MAP value K_b = 103 m s⁻¹) and P_b = 4.915 ± 0.004 d (P_{b, injected} = 4.921 d); the model including the planetary signal was only moderately favoured over the model with just the correlated stellar activity signal, with a Bayes factor of about 5, which is not enough to claim a statistically significant detection. When K_b = 130 m s⁻¹ is used, the planetary signal was recovered far more precisely ($K_{\text{b}}={136}_{-16}^{+14}$ m s⁻¹ and P_b = 4.932 ± 0.002 d, with P_{b, injected} = 4.933 d), and with a high significance (Bayes factor of about 9 × 10⁶). This simple test demonstrates that we can reliably detect the planet when the semi-amplitude of the injected signal is greater than the RMS of the RV residuals of the real data when the quasi-periodic activity signal is removed (see Table 2).

7.2 Injection-recovery simulations under a more general scheme

For the second set of simulations, we kept the same time stamps as for the real data, and generated 100 mock datasets as follows. We added the planetary signal of DO17 to the best-fit stellar activity signal that we determined through a GP regression including one planet. We used the error bars of each GP hyper-parameter and of the uncorrelated jitter σ_jit to randomly draw arrays of parameters from a multi-dimensional normal distribution. The arrays of hyper-parameters were used to generate the quasi-periodic stellar activity term, to which we added the planetary signal (with K_b, P_b, and T_conj,b drawn from normal distributions, as done before). Finally, we added a randomly generated white-noise term with an RMS equal to that of the residuals of the real data (105 ms⁻¹), and the RV values obtained in this way were randomly shifted within the internal errors given as $\sqrt{{\sigma }_{\text{RV}}^2+{\sigma }_{\text{jit}}^2}$ (with ${\sigma }_{\text{jit}}^2$ = 115 m s⁻¹), still adopting a normal distribution.

Statistical properties of the simulated datasets are shown in Fig. 17. They should be compared with those of the real TERRA RVs from which they were derived (see periodograms and RMS values in Fig. 9). The mean GLS periodograms (left column) and the distributions of the RMS of the simulated data (original data and residuals determined through iterative pre-whitening) are on average well consistent with those of the real data.

We analysed each mock dataset with a GP regression including a sinusoid, and with the same setup used for the real RVs. The real σ_RV uncertainties were used as error bars. Figure 18 summarises some main outcomes of the analysis. Panel a shows examples of posterior distributions for the orbital period P_b, corresponding to datasets for which P_b was recovered with high and good accuracy, and one case for which P_b was not recovered at all. Panel b shows the distribution of the P_{b, retrieved}/P_{b, inj} ratio, and panel c shows the distribution of the P_{b, MAP}/P_{b, inj} ratio. For 61% of the samples, the best-fit median P_{b, retrieved} falls within ± 0.5 d of the injected orbital period, which we assume as the interval corresponding to an accurate and potentially precise (andsignificant) detection. The MAP value for nearly half of this sub-sample (corresponding to 32% of the total mock datasets, which we call the S₃₂ sample for convenience) lies within the same interval². The percentage of the total samples for which P_{b, MAP} falls within ± 0.5 d of the injected orbital period is 40%, which means that we failed to recover precise best-fit values of P_b for 8% of the samples with quite accurate MAP values. The median and mean significance of the recovered orbital period for the S₃₂ sample is 1.8σ and 164σ, respectively (see panel d of Fig. 18). Panel e shows the distribution of the ratio between the recovered and injected semi-amplitude K_b for sample S₃₂. For more than half of the samples, we recovered inaccurate values for K_b, which are less than half the injected value.

We also analysed each simulated dataset without including a sinusoid to model the signal due to the injected planet. We compared the Bayesian evidence $\ln {\mathcal{Z}}_{1pl}$ and $\ln {\mathcal{Z}}_{0pl}$ derived by MultiNest in order to assess how much the correct model (i.e. the model with the planetary signal) is statistically favoured. This in turn provided information on how effective our methods are at retrieving the injected signal. The result is shown in the last panel of Fig. 18. The model that includes one circular planetary signal is never more probable, except for two datasets only with $\ln {\mathcal{Z}}_{1pl}-\ln {\mathcal{Z}}_{0pl}>$2.

These simulations indicate that even with a large number (well over one hundred) of high-quality RVs from an instrument such as HARPS-N, it is almost invariably going to be extremely difficult to reliably detect the planet that was claimed by D017 to orbit V830 Tau.

Fig. 17 Statistical properties of the simulated dataset described in Sect. 7.2. The average GLS periodograms and the distributions of the RMS of the data (m/s; original data and residuals) are shown in the left and right columns, respectively. Each row, starting form the second, refers to residuals determined through iterative pre-whitening.

8 Discussion and conclusions

After the announcement of the discovery of a HJ with the radial velocity method, the ~2 Myr old star V830 Tau has become a milestone for our understanding of the formation and evolution timescales of extrasolar planets (Donati et al. 2016, 2017). The detection of a planet at an early stage of formation in a close-in orbit around its host (a = 0.057 au) revealed that Jupiter-like planets can migrate inwards in less than 2 Myr. This discovery motivated an RV follow-up campaign of V830 Tau within the GAPS programme, with the main goal of improving the planetary parameters using the high-resolution HARPS-N spectrograph.

With a variability of about 1 km s⁻¹ observed in the RV time series that is almost entirely due to magnetic activity, V830 Tau is one of the most active young stars monitored for blind planet searches. It therefore represents an a priori very challenging target even when the best currently available spectrographs are used, and claiming the detection of even a massive HJ with high statistical significance can be very difficult.

The conclusion from our analysis of HARPS-N RVs is that we cannot confirm the existence of the planetary signal attributed to V830 Tau b, which has been claimed with a very strong statistical significance by DO17 (Bayes’ factor of 10⁸ −10⁹). To investigate the presence of the planet as carefully as possible, we analysed RVs extracted with three different pipelines, and used different methods and tools to account for the dominant activity signals and perform robust Bayesian model comparison. Our analysis alsotook advantage of the information embedded in simultaneous activity diagnostics from photometry and spectroscopy. The internal errors of two of the HARPS-N RV datasets are almost half those of DO17, but the scatter in our measurements is higher. This could be due to an increase in the level of the stellar activity since 2016, and we found evidence foran increasing activity within the time span of our spectroscopic follow-up (Fig. 4). This may represent a further obstacle for recovering the planetary signal, despite the quality and sampling of our data. Figure 19 shows that our HARPS-N observations are distributed quite uniformly over the orbit of V830 Tau b. This means that our non-detection of the planet signal cannot be attributed to poor sampling.

We also devised detailed injection-retrieval simulations based on our data (TERRA dataset) and analysis set-up to meticulously investigate our sensitivity to the presence of a planetary companion. The main results are summarised below.

• After injecting the putative planetary signal into our real RV time series (Sect. 7.1), we were able to accurately recover the correct orbital period for almost 50% of the cases (4.9 < P_b < 5.0 days), but the semi-amplitude was not accurately retrieved. We were able to recover the planet with the same high statistical significance as was claimed by DO17 by only injecting a signal with twice the semi-amplitude of that reported in their work. However, it must be highlighted that our results depend on the adoption of quite broad, uninformative priors for the planetary model parameters in all our analyses, as expected when a blind search is conducted.

• After injecting the planetary signal of DO17 into randomly generated RV datasets with average properties similar to the real RVs (Sect. 7.2), we were able to retrieve accurate values for P_b (relying on the maximum a posteriori probability) for 32% of the realisations. For just eight of these favourable cases, however, were we able to recover a precise P_b (significance of the detection >5σ), while the retrieved semi-amplitude K_b was in general not accurate for this sub-sample. For all the mock datasets, except for one, the model comparison based on Bayesian marginal likelihoods showed that the model that includes the planetary signal is never statistically favoured.

While these results do not rule the existence of V830 Tau b out with high confidence, they nonetheless clearly show that retrieving the DO17 signal (or a signal very similar to it) with high significance is far from a simple task, even given the high quality of the HARPS-N data and the state-of-the-art tools and methods we used for the analysis.

We further explored the potential of our data by calculating the detection limits provided by the HARPS-N RVs. For this purpose, we derived a diagram showing the lowest minimum planetary mass we are sensitive to as a function of the planet orbital period. The calculation assumes the TERRA RV residuals of the GP model with N = 0 planets, and is based on the following frequentist approach. First, through a bootstrap analysis we derived the power pow_0.1% of the GLS periodogram corresponding to a level of false-alarm probability of 0.1%. We then defined arrays of velocity semi-amplitudes K, orbital periods P, and orbital phases ϕ to generate sinusoids that simulate signals induced by a planet on a circular orbit. We adopted 10–1000 m s⁻¹ and 2.2–440 d as the variability ranges for K and P, respectively (the upper limit on the orbital period is equal to half the time span of our observations), and 100 linearly spaced values between 0 and 1 were generated for ϕ. Each simulated sinusoid was injected into the original RV residuals, and we calculated the power pow_trial of the GLS periodogram at the planet orbital frequency. When given a pair (K,P), pow_trial > pow_0.1% for all the orbital phases, we consider the planet detected. We used 1 M_⊙ for the massof V830 Tau (from DO17) to transform the velocity semi-amplitudes into values of minimum mass for the injected planet. According to the results of this simplified but still illustrative calculation (Fig. 20), we straddle the detection limit for the m_p sin i_p = 0.57 M_Jup planet claimed by DO17. This agrees with the difficulties encountered in retrieving the planetary signal through the more complex and rigorous statistical analysis described in Sect. 7.1.

Our work was intended as an independent investigation of the V830 Tau system using HARPS-N, and therefore we do not present hereany reanalysis of the data from DO17. We cannot fully reject the reality of the 4.9 d signal claimed by DO17, but the HARPS-N observations and our analyses cast doubt on a planetary origin of the DO17 signal. Further work, new observations (even with NIR spectrographs, and possibly during epochs of lower stellar activity), more sophisticated analysis techniques, and/or perhaps a better understanding of nuisance signals in existing RVs are clearly required to definitively confirm or refute the existence of V830 Tau b. This point is of crucial relevance for assessing the occurrence rate of HJs around young stars, and for understanding the formation paths and migration mechanisms that apparently might cause them to move close to their newborn star on a short timescale, before the dissipation of the protoplanetary disc.

One main conclusion of our work is that any detection based on the RVs alone should be taken with extreme caution, and that independent reanalysis and follow-up are strongly encouraged on a case-by-case basis. Good examples of debated RV-detectedHJs around young stars are TW Hya and Cl Tau. The first is the closest T Tauri star to the Sun (~10 Myr old), with a candidate giant planet (~10 M_Jup) detected at a separation of 0.04 au (P ~ 3.5 d) by Setiawan et al. (2008). The existence of this close-in companion was debated by Huélamo et al. (2008), who concluded that the RV signal could best be explained by a long-lasting cool stellar spot on the stellar surface. The star Cl Tau is coeval to V830 Tau with the first HJ candidate (P ~ 9 d) detected within the very young protoplanetary disc (Johns-Krull et al. 2016). This detection has recently been questioned and the signal attributed instead to stellar activity (Donati et al. 2020). Interestingly, for Cl Tau there is evidence for ongoing giant planet formation at larger separations (10–100 au), as revealed by Clarke et al. (2018) using high-resolution imaging with the Atacama Large Millimeter/submillimeter Array (ALMA). Moreover, the existence of a HJ orbiting the older ~ 150 Myr star BD+20 1790 was ruled out by the GAPS collaboration using NIR RVs (Carleo et al. 2018), and still GAPS observations, using the combination of HARPS-N (VIS) and GIANO (NIR) RVs, enabled excluding the existence of an HJ orbiting AD Leo (age between 25 and 300 Myr) (Carleo et al. 2020). Detecting young planets in close-in orbits with the photometric transit method remains the more secure way today to ascertain their existence. Precisely characterising them with spectroscopic follow-up observations is still challenging, however.

Fig. 18 Panel a: posterior distributions of the orbital period Pb of the injected planet for some selected mock datasets. The histogram in green corresponds to a dataset for which Pb has been recovered with high precision and accuracy. The posterior in blue corresponds to a case with a very accurate but less precise retrieved Pb, and the posterior in red corresponds to a mock dataset for which the orbital period was not recovered. Panel b: distribution of the Pb, retrieved/Pb, inj ratio between the median best-fit orbital period and the injected value for each simulated dataset. Panel c: same as for panel b, but using the values of Pb corresponding to the MAP estimates. Panel d: distribution of the significance of the retrieved best-fit values Pb for the S32 sub-sample described in the text. The significance is expressed as the ratio between the median of the posterior and the corresponding lower uncertainty for each simulated dataset. Panel e: distribution of the ratio between the recovered and injected semi-amplitude Kb of the planetary RV signal for the 32% of the simulated datasets for which an accurate and precise estimate of Pb was recovered (S32 sample). Panel f: results of a model comparison analysis. The plot shows the distribution of the difference $\ln {\mathcal{Z}}_{1pl}$-$\ln {\mathcal{Z}}_{0pl}$, calculated for the 100 mock RV datasets, between the Bayesian evidence of the model including one sinusoid (i.e. one planet on a circular orbit) and the model with only a GP-term. The model including one planet is marginally more significant ($\ln {\mathcal{Z}}_{1pl}$-$\ln {\mathcal{Z}}_{0pl}>$2) only in two cases, one being that corresponding to the posterior of Pb in green shown in panel a.

Fig. 19 Distribution of the number of HARPS-N observations of V830 Tau according to the orbital phase of the planet announced by Donati et al. (2017).

Fig. 20 Detection limits for planets on circular orbits based on the TERRA RV residuals for the GP model with N = 0 planets. For each trial orbital period, we calculated the minimum detectable value of the planet minimum mass. The red star identifies the position on the diagram of the planet announced by Donati et al. (2017).

Acknowledgements

We thank the anonymous referee for her/his consideration and useful suggestions. We thank Daniele Locci (INAF-OAPa) for useful comments. We acknowledge the support by INAF/Frontiera through the “Progetti Premiali” funding scheme of the Italian Ministry of Education, University, and Research. We acknowledge the computing centres of INAF – Osservatorio Astronomico di Trieste / Osservatorio Astrofisico di Catania, under the coordination of the CHIPP project, for the availability of computing resources and support. We thank Chris Sneden for providing us with the line list useful for measuring stellar parameters. V.M.R. thanks the Royal Astronomical Society and Emmanuel College, Cambridge, for financial support. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier).

Appendix A Light-curve measurements

Appendix B Spectroscopic activity indexes

Appendix C Radial velocity measurements

Appendix D Radial velocity analysis

Fig. D.1 Posterior distributions of the free (hyper-)parameters of the GP+1 circular planet regression (Table 2). The corresponding RV dataset is that extracted with TERRA pipeline.

References

All Tables

All Figures

	Fig. 1 Kepler/K2 light curve. The time series collected during the first half of 2017, and the same data phase-folded to the stellar rotation period are shown in the upper and lower panel, respectively. The most powerful flares occurred at phases close to maximum brightness.
In the text

	Fig. 2 Detrended and flattened K2 light curve of V830 Tau, phase-folded at the orbital period of the planet announced by DO17. Red dots represent five-point binned data.
In the text

	Fig. 3 STELLA light curves in V and I bands, spanning the last two seasons of HARPS-N observations, and their respective GLS periodograms. The peak corresponding to the rotation period is labelled.
In the text

	Fig. 4 Time series of the activity indicators extracted from the HARPS-N spectra based on the CaII H&K and H-α spectral lines (left panels), and their corresponding GLS periodograms (right panels). The peak corresponding to the rotation period is labelled.
In the text

	Fig. 5 Correlation plot between the CCF BIS index and the RVs extracted with the DRS pipeline (see Sect. 5). Left panel: all data. Right panel: data without outliers to better visualise the RV-BIS correlation.
In the text

	Fig. 6 Residuals of the CCF, i.e. each individual CCF is divided by the average CCF calculated over the whole dataset, phase-foldedto the rotational period of V830 Tau (2.741 d).
In the text

	Fig. 7 Cross-correlation function of V830 Tau for a HARPS-N spectrum with average S/N. The CCF was calculated by the DRS pipeline using a template mask for a KV5 star.
In the text

	Fig. 8 Radial velocities of V830 Tau extracted from the HARPS-N spectra with different methods (average subtracted). From top to bottom: DRS, TERRA, and template-free algorithm (R20) by Rajpaul et al. (2020). The error bars are nearly two orders of magnitude smaller than the major ticks on the y-axis, and are not visible. For each season we indicate the RMS of the corresponding RV subsample.
In the text

Fig. 9 GLS periodograms of the original TERRA (a), DRS (b), and R20 (c) RVs and their residuals after recursive pre-whithening. For each figure: the vertical and dashed blue lines indicate the stellar rotation frequency and its first and second harmonic; the red line marks the orbital frequency of the planet announced by Donati et al. (2016, 2017); each panel reports the RMS of the dataset used for calculating a specific periodogram; the window function is shown in green in the bottom panel.

In the text

	Fig. 10 STELLA light curves (upper panel), CaII H&K and H-α activity indicators (middle panel), and TERRA RVs (lower panel) phase folded using the same period P = 2.741 d and reference epoch of the STELLA V band photometry andKepler/K2 light curve in Fig. 1. Subplots for the light curves and RVs show the residuals O–C of the best-fit sine function obtained with the GLS software.
In the text

Fig. 11 V band optical light curves of V830 Tau with the flux plotted vs. rotation phase. The data points collected between BJD 58 393.742 and 58 607.363 are plotted as open green diamonds and their kernel regression vs. the rotational phase is given by the dashed green line, while the data points collected between BJD 58 765.645 and 58 916.387 are plotted as open orange triangles and their regression is given by the dot-dashed orange line. The solid red line is the regression to the whole V band photometric dataset vs. the rotation phase.

In the text

Fig. 12 Top panel: RVs obtained with the TERRA procedure vs. time (filled dots). The model RV variation as derived from the solid red interpolation in Fig. 11 by means of the FF′ method is superposed (solid red line). Bottom panel: residuals between the observed TERRA RV and the RV predicted by the FF′ model vs. time. In both panels the size of the RV error bars is comparable with that of the data points in most of the cases.

In the text

Fig. 13 Top panel: RV time series of V830 Tau as extracted with the TERRA procedure vs. the rotation phase (filled dots). Different colours indicate data collected in different seasons: black dots indicate data points collected between BJD 58 044.6623 and 58 192.3573 (first season), red dots between BJD 58 341.7344 and 58 566.3847 (second season), and green dots between 58 804.4827 and 58 924.4063 (third season). The KR performed vs. the rotation phase and time in the first season is indicated by the solid black line, in the second season by the solid red line, and in the third season by the solid green line. The RV error bars are smaller than the size of the plotted dots. Bottom panel: RV residuals obtained by subtracting the seasonal KRs from the corresponding data points.

In the text

Fig. 14 Top panel: cleaned RV time series of V830 Tau vs. time as obtained after the first KR with respect tothe rotation phase and time (filled dots). The solid red line indicates a further KR with respect to the phase and time that gives the largest reduction of the standard deviation of the residuals in comparison with the KRs computed with respect to the CCF indicators or the chromospheric index $\log R_{\text{HK}}^{\prime }$. The dashed green lines indicate the 1σ confidence range of the KR (see Lanza et al. 2018, for details). Bottom panel: residualsof the KR in the top panel vs. time, with a standard deviation of 65.13 m s−1.

In the text

	Fig. 15 Generalized Lomb-Scargle periodogram of the residual RV time series of V830 Tau in the bottom panel of Fig. 14. The power is normalised to its maximum in the considered period interval and is plotted vs. the period itself.The vertical dashed green line marks the rotation period, and the dotted red line indicates the orbital period of the planet proposed by Donati et al. (2017).
In the text

	Fig. 16 Results for the simulations with direct injection of the planetary signal into the real data.
In the text

	Fig. 17 Statistical properties of the simulated dataset described in Sect. 7.2. The average GLS periodograms and the distributions of the RMS of the data (m/s; original data and residuals) are shown in the left and right columns, respectively. Each row, starting form the second, refers to residuals determined through iterative pre-whitening.
In the text

Fig. 18 Panel a: posterior distributions of the orbital period Pb of the injected planet for some selected mock datasets. The histogram in green corresponds to a dataset for which Pb has been recovered with high precision and accuracy. The posterior in blue corresponds to a case with a very accurate but less precise retrieved Pb, and the posterior in red corresponds to a mock dataset for which the orbital period was not recovered. Panel b: distribution of the Pb, retrieved/Pb, inj ratio between the median best-fit orbital period and the injected value for each simulated dataset. Panel c: same as for panel b, but using the values of Pb corresponding to the MAP estimates. Panel d: distribution of the significance of the retrieved best-fit values Pb for the S32 sub-sample described in the text. The significance is expressed as the ratio between the median of the posterior and the corresponding lower uncertainty for each simulated dataset. Panel e: distribution of the ratio between the recovered and injected semi-amplitude Kb of the planetary RV signal for the 32% of the simulated datasets for which an accurate and precise estimate of Pb was recovered (S32 sample). Panel f: results of a model comparison analysis. The plot shows the distribution of the difference $\ln {\mathcal{Z}}_{1pl}$-$\ln {\mathcal{Z}}_{0pl}$, calculated for the 100 mock RV datasets, between the Bayesian evidence of the model including one sinusoid (i.e. one planet on a circular orbit) and the model with only a GP-term. The model including one planet is marginally more significant ($\ln {\mathcal{Z}}_{1pl}$-$\ln {\mathcal{Z}}_{0pl}>$2) only in two cases, one being that corresponding to the posterior of Pb in green shown in panel a.

In the text

	Fig. 19 Distribution of the number of HARPS-N observations of V830 Tau according to the orbital phase of the planet announced by Donati et al. (2017).
In the text

	Fig. 20 Detection limits for planets on circular orbits based on the TERRA RV residuals for the GP model with N = 0 planets. For each trial orbital period, we calculated the minimum detectable value of the planet minimum mass. The red star identifies the position on the diagram of the planet announced by Donati et al. (2017).
In the text

	Fig. D.1 Posterior distributions of the free (hyper-)parameters of the GP+1 circular planet regression (Table 2). The corresponding RV dataset is that extracted with TERRA pipeline.
In the text