The T Tauri star V410 Tau in the eyes of SPIRou and TESS

We performed spectropolarimetric observations with SPIRou between 2019 October 31 and 2019 December 13 (UTC). SPIRou works in the near-infrared (NIR) domain between 950 nm and 2500 nm with a spectral resolving power of $\mathord{\sim}$70 000 (Donati et al., 2020). Each polarimetric observation is composed of a sequence of 4 subexposures of 300 seconds each taken with different configurations of the polarimeter (i.e., different azimuths of the polarimeter retarders) in order to get rid of potential spurious signals in the polarisation and systematic errors at first order (Donati et al., 1997). Twenty sequences were collected, yielding spectra in both unpolarized (Stokes $I$) and circularly polarized (Stokes $V$) light. Data reduction was performed with a pipeline based on the ESPaDOnS pipeline Libre-ESpRIT (Donati et al., 1997) adapted for SPIRou observations. We then obtained telluric-corrected spectra using a PCA approach as mentionned by Artigau et al. (2014). The signal-to-noise ratio (SNR) per pixel of these spectra peaks in the $H$ band with a median value of 179 (ranging from 140 to 200). We applied Least-Square Deconvolution (LSD; Donati et al. 1997) on all spectra in order to add up information from all lines. We only used 18 out of the 20 recorded spectra, the two remaining ones (collected on 2019 Nov 10 and Dec 13) being either much noiser than the average or suffering from a technical issue.

Three series of LSD profiles were computed with three different masks generated with the VALD-3 database (Ryabchikova et al., 2015) and including lines ranging from 950 to 2600 nm. The first one (hereafter M1) contains $\mathord{\sim}$10 000 atomic and molecular lines. It yields Stokes $I$ LSD profiles only, as the magnetic sensitivity (i.e. Landé factor) of many molecular lines is unknown, with SNRs ranging from 1890 to 3060 (median of 2770). The second mask (hereafter M2) contains $\mathord{\sim}2000$ atomic lines with known Landé factors and relative depths with respect to the continuum >10%. It yields Stokes $I$ LSD profiles with SNRs ranging from 1270 ro 1930 (median of 1740) and Stokes $V$ LSD profiles with SNRs ranging from 3380 to 4950 (median of 4440). The last mask (M3 hereafter), containing $\mathord{\sim}900$ molecular lines with a depth relative to the continuum down to 5%, yields Stokes $I$ LSD profiles with SNRs ranging from 600 to 1020 (median of 835).

To phase our observations on the rotation cycle of the star, we used the same reference date as in Skelly et al. (2010) and Yu et al. (2019), namely the Barycentric Julian Date BJD₀=2,454,832.58033, along with the well defined stellar rotation period obtained by Stelzer et al. (2003): $P_{\mathrm{rot}}~{}=~{}1.871970~{}\pm~{}0.000010$ d. Despite the very well constrained rotation period, the reconstructed maps presented in Sec. 3 cannot be directly compared to the previously published ones, as the brightness and magnetic surface maps of V410 Tau evolved since then as a result of both differential rotation and intrinsic variability.

A summary of our observations is given in Table 1. Since our data are spread over 23 rotation cycles only (i.e. 43 d), we do not expect major changes in the surface brightness map nor in the large-scale magnetic topology of V410 Tau, given the typical timescale on which both quantities are found to evolve, either as a result of intrinsic variability or differential rotation (Yu et al., 2019).

Contemporaneously with SPIRou observations, V410 Tau (TIC 58231482) was observed by TESS (Ricker et al., 2014) from 2019 November 28 to 2019 December 23 during its monitoring of sector 19. TESS being mostly sensitive to wavelengths in the range 600 to 1000 nm (centred on the Cousins $I_{\rm c}$ photometric band), it thus probes a different surface brightness distribution than that observed with SPIRou.

The target was observed with a cadence of 2 min over a total time span of 25 d. These observations have been re-processed by the Science Processing Operations Center (SPOC; Jenkins et al. 2016) data pipeline (version 4.0), with light curves available from the Mikulski Archive for Space Telescopes (MAST)¹¹1https://archive.stsci.edu/missions-and-data/tess. In particular, we used the Pre-search Data Conditioning Single Aperture Photometry (PDCSAP) flux that provides a better estimate of the intrinsic variability of the star since instrumental variations have been corrected in this light curve, as well as contamination from some nearby stars (Smith et al., 2012; Stumpe et al., 2012; Stumpe et al., 2014). We only kept the data that were not flagged from the SPOC pipeline, thus those with no known problems. We also rejected the observations carried out between BJD 2,458,826.5 and 2,458,827.8, for which a high background level from the Earth led to spurious photometric variations²²2Details can be found in the Data Release notes of sector 19 (DR26 and DR30) available at https://archive.stsci.edu/tess/tess_drn.html (see Fig. 1).

Since TESS is sensitive to small flux variations caused by flares, we first filtered them through a $3\sigma$-clipping process, which also allowed us to mitigate the potential effects of activity from stellar companions (V410 Tau B and C) on the photometry. This process consisted in an iterative fit of the light curve with a Gaussian Process Regression (GPR; Rasmussen & Williams 2006). At each iteration, we rejected the points having a residual larger than 3$\sigma$ until all the residuals were lower than this threshold. As we expected that the variations in the light curve are mainly due to stellar activity, we chose a quasi-periodic kernel to model the TESS light curve, which is well adapted to model signals induced by active regions rotating with the stellar disk (Haywood et al., 2014):

where $t_{i},t_{j}$ are the times associated to the observations $i$ and $j$, respectively. $\theta_{1}$ is the amplitude of the Gaussian Process (GP), $\theta_{2}$ is the decay timescale (i.e. the exponential timescale on which the modeled photometry departs from pure periodicity) or equivalently typical spots lifetime, $\theta_{3}$ is the period of the GP (expected to be close to $P_{\mathrm{rot}}$) and $\theta_{4}$ is the smoothing parameter that controls the amount of short-term variations that the fit can include (within a rotation cycle). For our purpose, we imposed a large decay timescale $\theta_{2}=300$ d to avoid the GP to fit rapidly evolving patterns. The resulting light curve is shown in Fig 1. We then modelled this filtered light curve with the same GPR but letting all 4 parameters free in order to derive an estimate of the typical timescale on which the main surface features evolve (found to be equal to $162^{+30}_{-25}$ d). This value is lower than that obtained by Yu et al. (2019) from ground-based $V$ measurements ($314^{+31}_{-29}$ d) most likely because TESS photometry is much more accurate and sensitive to small structures evolving rapidly. We also found $\theta_{3}=1.873\pm 0.001$ d, consistent with the rotation period obtained by Stelzer et al. (2003) and the smoothing parameter $\theta_{4}=1.02\pm 0.06$. Considering only the data obtained before the end of SPIRou observations (2019 December 12), we computed median time and relative photometry every ten points, resulting in 757 photometric values, in order to get a smoother curve, to reduce computational time and to balance the relative weights of photometry and spectroscopy when applying ZDI (see Sec. 3.2).

Despite the very-high-precision photometry provided by TESS, no periodic signal beyond that due to V410 Tau A is detected in the light curve, which indicate that the temporal variations are attributable to the main star. In addition, the amplitude ratio between the peaks associated with component A and the noise in the TESS light curve periodogram (equal to $\mathord{\sim}32$) provides an independant confirmation that the magnitude contrasts with the two other components in the $I_{\rm c}$ band is about 4 or more, in agreement with Ghez et al. (1997).

Multicolour photometry of V410 Tau was collected with the ground-based 1.25 m AZT-11 telescope at the CrAO from 2019 September 02 to December 18. Using a ProLine PL23042 CCD camera, 40 brightness measurements were collected in the $V$, $R_{\rm c}$ and $I_{\rm c}$ bands. For these estimates, the wTTS V1023 Tau was used as a comparison star as this star shows small variability amplitude. We note that V410 Tau has an average magnitude in the $V$ band equal to 10.85 mag (see Fig. 2), thus consistent with observations of this star at previous epochs (Grankin et al., 2008). A full log of the CrAO photometric observations is given in Table 5.

We fitted the light curves in the $V$, $R_{\rm c}$ and $I_{\rm c}$ bands separately, using a periodic fit including the fundamental frequency and the first two harmonics. As the formal uncertainties on the measured magnitudes (7 mmag for $V$ and $R_{\rm c}$, 5 mmag for $I_{\rm c}$) are underestimated, missing some sources of noise like the intrinsic variability of the observed star, we derived empirical estimates of these error bars by scaling them up to the values that ensure a unit reduced chi square for the fit in each band (using a sine wave plus 2 harmonics, see Fig. 2). We find these empirical error bars to be equal to $20$, $14$, and $13$ mmag for $V$, $R_{\rm c}$ and $I_{\rm c}$, respectively. The models for $V-R_{\rm c}$ and $V-I_{\rm c}$ were then obtained by subtracting the models for the associated magnitudes while the error bars on these colours were estimated by propagating the uncertainties (Fig. 2).

As expected we find that the light curve amplitude decreases with wavelength: $0.231\pm 0.012$ mag in the $V$ band, $0.212\pm 0.008$ mag in the $R_{\rm c}$ band and $0.174\pm 0.008$ mag in the $I_{\rm c}$ band (with the error bars on the amplitudes derived from the empirical uncertainties on the measured magnitudes). As expected, TESS photometry (amplitude of $0.17656\pm 0.0009$ mag) is consistent with that obtained from the ground in the $I_{\rm c}$ band. As all light curves are more or less consistent in shape, we will only use TESS data in the following as they are much more accurate than ground-based photometry.

To recover maps of the surface brightness and/or of the magnetic field topology of the star, we applied ZDI (Semel, 1989; Brown et al., 1991; Donati & Brown, 1997; Donati et al., 2006) on our time series of Stokes $I$ and/or Stokes $V$³³3Stokes $V$ profiles are obviously different from $V$ magnitudes despite sharing similar notations. As there is no real ambiguities between the 2 quantities, we kept the notations unchanged. LSD profiles. ZDI aims at constraining surface (brightness or magnetic) maps of rotating stars from time series of (Stokes $I$ and $V$) LSD profiles collected as the star rotates. To achieve this we proceed as follows. ZDI uses a conjugate gradient algorithm to deduce iteratively, from an initially non-spotted (resp. non-magnetic) distribution, the maps of relative brightness with respect to the quiet photosphere at SPIRou wavelengths (resp. of the magnetic field vector) at the surface of the star, until the corresponding synthetic LSD profiles fit the observed ones down to a unit reduced chi-square ($\chi^{2}_{r}$). This allows to look for the maximum-entropy solution of this optimization problem, i.e. the map containing the smallest amount of information capable of fitting the data down to the noise level.

In a second step, we applied ZDI to both SPIRou and TESS data simultaneously. We proceeded as in Yu et al. (2019) but this time including photometric data as part of the fit (instead of simply comparing the light curve predicted by ZDI with photometric observations). This is achieved by deriving, as part of the iterative imaging process and using Planck’s law, the brightness contrast in the TESS bandpass that we expect from the one in the SPIRou bandpass (which we reconstruct through ZDI).

In practice, we divide the surface of the star into a grid of a few thousand cells; to estimate the spectral contribution of each cell to the measured Stokes $I$ and $V$ LSD profiles, we use Unno-Rachkovsky’s analytical solution to the polarized radiative transfer equations in a Milne-Eddington atmosphere (see e.g. Landi Degl’Innocenti & Landolfi 2004) with appropriate values for relevant parameters such as the limb-darkening coefficient known to strongly depend on wavelength ($0.3\pm 0.1$ in the $H$ band for $T_{\rm eff}=4500$ K and $\log g=4$; Claret & Bloemen 2011). The synthetic Stokes $I$ (resp. Stokes $V$) LSD profiles are then computed by integrating all the local Stokes $I$ (resp. Stokes $V$) LSD profiles over the visible stellar hemisphere, while the photometric values are computed by summing up the value of the continuum over all grid cells.

The poloidal and toroidal components of the magnetic field are decomposed into spherical harmonics (Donati et al., 2006) while the photospheric relative brightness is computed independently for each cell of the grid.

In a first approach, we will assume that the star rotates as a solid body. We will then take into account differential rotation at the surface of V410 Tau in Sec. 3.4.

In order to fit our LSD profiles, we chose a line model with mean wavelength, Doppler width and Landé factor of 1650 nm, 1.8 $\mathrm{km}\,\mathrm{s}^{-1}$ and 1.2, respectively. As the depth of the LSD profiles varies depending on the mask and as we kept constant the Doppler width, our line models features an equivalent width (EW) of 1.5 $\mathrm{km}\,\mathrm{s}^{-1}$, 1.9 $\mathrm{km}\,\mathrm{s}^{-1}$ and 1.4 $\mathrm{km}\,\mathrm{s}^{-1}$ for the M1, M2 and M3 masks, respectively.

Using ZDI on our Stokes $I$ LSD profiles obtained from SPIRou data, we found $v\sin{i}~{}=~{}72.8~{}\pm~{}0.5$  $\mathrm{km}\,\mathrm{s}^{-1}$ for the line-of-sight-projected equatorial rotation velocity and $i=45\pm 10^{\circ}$ for the inclination of the rotation axis to the line of sight. These values being consistent with those of Yu et al. (2019) within the error bars, we decided to follow Yu et al. (2019) and set $v\sin{i}=~{}73.2$ $\mathrm{km}\,\mathrm{s}^{-1}$ and $i=50^{\circ}$. Using ZDI, we found that the bulk RV of the star is equal to $17.4\pm 0.3$ $\mathrm{km}\,\mathrm{s}^{-1}$ when considering M1, $18.1\pm 0.3$ $\mathrm{km}\,\mathrm{s}^{-1}$ when considering M2 and $16.0\pm 0.6$ when considering M3. Although these values remain compatible within 3$\sigma$, we suspect that these differences may partly come from inaccuracies in the empirical mask line wavelengths, known to be less reliable for molecules than for atoms. It may also reflect a systematic RV blueshift of molecular lines with respect to atomic lines, that would suggest that atomic lines are more affected than molecular lines by the inhibition of convective blueshift by stellar magnetic activity for a reason yet to be clarified.

Finally, as in Yu et al. (2019), we set the maximum number of spherical harmonics to $\ell=15$ to describe the large-scale magnetic topology.

A summary of the information provided by our ZDI analyses with the three different masks is provided in Table 2.

We focused on the Stokes $I$ LSD profiles obtained with M1 to deduce the surface brightness map of the star. We performed two analyses, one with SPIRou data only and one considering simultaneously SPIRou and TESS data.

For the first analysis, we just applied ZDI on our time series of Stokes $I$ LSD profiles that we fitted down to $\chi^{2}_{r}=1$ (see Fig. 3). The reconstructed surface brightness distribution is shown in the left panel of Fig. 4. We note the presence of large features, such as dark spots at phase 0.70 or 0.85, along with smaller ones. However, rather surprisingly, no polar spot is visible in this map from SPIRou data, whereas previous maps from optical data always showed conspicuous dark features covering the polar regions (e.g. Joncour et al. 1994; Hatzes 1995; Yu et al. 2019).

For our brightness reconstruction using both SPIRou and TESS data, we provided ZDI with the 757 photometric values (with error bars set to 1.6 mmag per data point) from the TESS light curve before 2019 December 12 (see Sec. 2.2) in addition to the Stokes $I$ LSD profiles. This yielded the reconstructed map and the fitted light curve shown in the right panel of Fig. 4 and in Fig. 5, respectively. Both spectroscopic and photometric data were fitted down to $\chi^{2}_{r}=~{}1$, with the fit of Stokes $I$ profiles being almost identical to that obtained using SPIRou data only. We find that including photometry in addition to SPIRou data yields brightness maps with enhanced contrasts, especially at low latitudes, that were needed to fit the light curve at a RMS level of 1.6 mmag.

When fitting Stokes $I$ profiles only, the spot coverage is close to 8.4% while it increased to 10.7% when TESS data were taken into account.This 2.3% increase is significant, firstly because this difference is larger than the typical error we can expect on this parameter (of order 0.4% in the context of this particular data set), and secondly since maximum entropy is profiled to provide the image containing the smallest amount of information. We come back on the origin of this difference in Sec. 5.1. We note that, in the fitted light curve, the residuals still exhibit correlated noise (at a level of 1.6 mmag RMS) that does not repeat from one rotation cycle to the next, and that likely reflects small-scale structures at the surface that evolve with time. (This explains a posteriori why the error bar on the TESS data points was set to this value).

As no polar spot is visible whether we take into account the TESS data or not, we reconstructed the brightness distribution from Stokes $I$ LSD profiles using the M3 mask, expected to be more sensitive to cooler regions (Fig. 6 and 19). Once again, the reconstructed map does not show a polar spot but displays some low-level differences with the map reconstructed with the M1 mask, e.g., the cool spot at phase 0.70 (see Fig 6).

Obviously, we do expect the brightness distribution to evolve over timescales of several years. However, we expect histograms of brightness contrasts to remain more or less the same in a given spectral range. In this context we can compare maps from NIR and optical data by rescaling, with Planck’s law, brightness maps reconstructed from optical data (see left panel of Fig. 7 in the particular case of epoch 2013 December; from Yu et al. 2019) to the image one would have reconstructed in the NIR (right panel of Fig. 7). We find that, as expected, the NIR map is less contrasted than the original optical map. Moreover, the histogram of rescaled optical data is comparable with that directly obtained from SPIRou data (Fig. 8).

To reconstruct the large-scale magnetic topology, we focused on LSD profiles provided by mask M2, as it contains only lines with well known magnetic sensitivity.

We fitted simultaneously Stokes $I$ and Stokes $V$ LSD profiles with ZDI. Once again, we achieved a fit of the data down to $\chi^{2}_{r}=1$. The profiles and the associated maps are shown in Fig. 9 and 10, respectively. We found a map of the logarithmic surface brightness similar to those reconstructed with mask M1 (with and without the TESS data). We see that the magnetic field is complex, with a topology similar to those derived by Yu et al. (2019). We note that the magnetic field average strength is close to 410 G. The poloidal component of the field, which encloses about 60% of the overall reconstructed magnetic energy, is essentially non-axisymmetric ($\mathord{\sim}50\%$) and weakly dipolar ($\mathord{\sim}35\%$) while the toroidal one presents the inverse properties ($\mathord{\sim}70\%$ axisymmetric and $\mathord{\sim}55\%$ dipolar). We find that the dipole component has a polar strength of 390 G and that its axis is tilted at 15^∘ to the rotation axis, towards phase 0.70.

We see no clear correlation between brightness and magnetic field maps obtained with ZDI as was already the case in maps derived from optical data (Yu et al., 2019).

We determined the longitudinal field⁴⁴4The longitudinal field is defined as the (algebraic) line-of-sight projected magnetic field component averaged over the visible hemisphere and weighted by brightness inhomogeneities. $B_{\ell}$ at each epoch by computing the first moment of Stokes $V$ profiles (Donati et al., 1997). The longitudinal field varies between about +50 and -250 G. The associated uncertainties range between 29 G and 42 G with a median of 32 G. We clearly see a periodic pattern in this index that can be fitted down to the noise level with a sine curve with one harmonic (period of $1.873\pm 0.002$ d) as shown in Fig. 11.

As our spectropolarimetric observations are spread over more than one month, the data can potentially exhibit some variability due to differential rotation. By computing the amount of shearing by latitudinal differential rotation that the surface (brightness or magnetic) maps experienced with time, ZDI allows one to estimate the surface differential rotation assuming a solar-like shear, given by:

where $\theta$ is the colatitude, $\Omega_{\mathrm{eq}}$ and $\mathrm{d}\Omega$ are the rotation rates at the equator and the difference of rotation rate between the pole and the equator, respectively. This differential rotation law was found to be successful at modeling the surface shear of low-mass stars, including those of rapidly rotating, fully convective dwarfs (Morin et al., 2008), including in particular V410 Tau (Yu et al., 2019). One can measure both $\Omega_{\mathrm{eq}}$ and $\mathrm{d}\Omega$ by finding out the values that minimize the $\chi^{2}_{r}$ for a given amount of reconstructed information.

The $\chi^{2}_{r}$ maps derived from Stokes $I$ and $V$ data respectively are shown in Fig. 12, where contours of 68% (1$\sigma$) and 99.7% (3$\sigma$) confidence levels are depicted. We find from Stokes $I$ LSD profiles that $\Omega_{\mathrm{eq}}=3358.8\pm 0.5$ mrad d^-1 and $\mathrm{d}\Omega=6.4\pm 2.2$ mrad d^-1, while Stokes $V$ LSD profiles yield $\Omega_{\mathrm{eq}}=3358.7\pm 0.4$ mrad d^-1 and $\mathrm{d}\Omega=9.0\pm 1.9$ mrad d^-1, both estimates being mutually consistent within 1.5 $\sigma$. This implies that the rotation period ranges from about 1.8707 d ($\pm 0.0003$ from Stokes $I$ and $\pm 0.0002$ from Stokes $V$) at the equator to $1.8742\pm 0.0013$ d (from Stokes $I$) or $1.8757\pm 0.0011$ d (from Stokes $V$) at the pole. Compared to previous shear detections in the optical (Yu et al., 2019), we find that the uncertainty on $\Omega_{\rm eq}$ and d$\Omega$ obtained from Stokes $I$ are about 5 and 3 times larger, respectively, which is likely related to the lower amount of data in our set as well as to the lower contrast of the reconstructed brightness features. On the contrary, uncertainties derived from Stokes $V$ profiles are comparable to that estimated from optical data at previous epochs.

We also computed the colatitude corresponding to the barycentre of the brightness and magnetic field distributions from the slope of the major axis of the confidence ellipse. In particular, our estimate from Stokes $I$ LSD profiles is slightly larger than the ones derived by Yu et al. (2019) in the optical at previous epochs, suggesting that large surface features are indeed migrating poleward as speculated by these authors.

Table 3 gathers our results about differential rotation.

We computed the RV of V410 Tau at each observed phase as the first moment of Stokes $I$ LSD profiles (Donati et al., 2017) for both observed (obtained with M1) and ZDI synthetic ones. With the set of synthetic profiles, we simulated noisy profiles with the same SNRs as the observed ones. For several realisations of the noise we computed the RV and then the dispersion of these measurements to estimate the error bars on our actual radial velocity data. We compared raw observed RVs with synthetic RV curves obtained from the maps in Fig. 4 and we computed activity filtered RVs as the difference between the observations and the model (Fig. 13). We note that both maps lead to models that fit reasonably well the data.

In a second step, we used GPR to model the impact of stellar activity on the observed RVs. For that, we chose the same kernel as that presented in Sec. 2.2, in Eq. (1). Moreover, we added an additional term representing a potential excess of uncorrelated noise (in case our RV error bars are underestimated). The log likelihood function $\log\mathcal{L}$ we wanted to maximize becomes:

where $\mathrm{K}$, $\mathrm{\Sigma}$ and $\mathrm{S}$ denote the covariance matrix with a quasi-periodic kernel, the diagonal matrix containing the variance of the observed RVs, and the diagonal matrix containing the additional noise $s$ to the square. $N$ refers to the number of points (i.e. number of observed RVs) and $y$ corresponds to the observed raw RVs.

Given the low number of RV points, we chose to fix 2 of the 4 hyperparameters, namely the decay timescale $\theta_{2}$ (exponential timescale on which modeled RVs depart from pure periodicity) and the smoothing parameter $\theta_{4}$ (controling the amount of short-term variations in the fit) at 160 d (as derived from TESS light curve) and 0.35, respectively, following Yu et al. (2019)⁵⁵5As shown in Klein & Donati (2020), surface features induce more complex modulation on RVs than on light curves, hence the smaller value of the smoothing parameter $\theta_{4}$ used here compared to that derived when fitting the TESS data (see Sec. 2.2).. Through a Markov Chain Monte Carlo (MCMC) approach using the EMCEE PYTHON module (Foreman-Mackey et al., 2013), we sampled the posterior distribution of the other parameters given the priors listed in Table 4. We ran our MCMC on 5000 iterations of 100 walkers, and removed a burn-in period of 250 iterations, that is about 5 times larger than the autocorrelation time of the chain ($\mathord{\sim}50$ iterations). We then chose the median of these posterior distributions as best values for the free parameters. The resulting phase plot is illustrated in Fig. 14. The amount of excess uncorrelated noise in the data (modeled with $s$) is found to be compatible with 0. From the best set of parameters, we obtained the GPR fit that is shown in Fig. 13.

Each model yields a synthetic RV curve that we compared to our RV measurements; the corresponding $\chi^{2}_{r}$ values are equal to 0.53 when applying ZDI to SPIRou data alone, 1.97 when ZDI is applied to combined SPIRou and TESS data and 0.14 when applying GPR. The corresponding dispersion of the activity-filtered RVs is about twice lower when we use GPR (0.08 $\mathrm{km}\,\mathrm{s}^{-1}$) rather than ZDI (0.13 $\mathrm{km}\,\mathrm{s}^{-1}$ and 0.23 $\mathrm{km}\,\mathrm{s}^{-1}$, using SPIRou data alone or both SPIRou and TESS data). We come back on the potential origin of this difference in Sec. 5.5. We note that these dispersions are consistent with the typical uncertainty on our RV measurements ($\mathord{\sim}180$ $\mathrm{m}\,\mathrm{s}^{-1}$), demonstrating that our models are successful at reproducing the activity-induced RV variations, both for atomic and molecular lines (and despite the difference in bulk RVs for both sets of lines).

Applying the FF’ method (Aigrain et al., 2012) to the light curve predicted with ZDI, we can investigate the precision level at which this technique can mitigate activity in RV curves. We find that RV residuals exhibit a dispersion of 810 $\mathrm{m}\,\mathrm{s}^{-1}$ RMS, i.e., 6-7$\times$ larger than those predicted with ZDI (130 $\mathrm{m}\,\mathrm{s}^{-1}$ RMS), confirming that the FF’ technique is not adequate for filtering out RV curves of moderately to rapidly rotating active stars, whose brightness distributions are often rather complex.

We focused on three main activity indicators: the He i triplet at 1083 nm, the Paschen $\beta$ (Pa$\beta$) and Brackett $\gamma$ (Br$\gamma$) lines at 1282 nm and 2165 nm, respectively (Fig. 20).

We identified a flare on December 09 (cycle 21.752) with the corresponding spectrum being blueshifted by 100 $\mathrm{km}\,\mathrm{s}^{-1}$ with respect to the stellar rest frame, and the flux in the three lines being stronger than the typical one, especially in the Pa$\beta$ line (Fig. 20). This flare occurs in the main gap of the TESS data (BJD 2,458,826.5 to BJD 2,458,829.3, see Fig. 1) and therefore does not show up in the light curve. In addition, two observations collected just after the flare (i.e. on December 10 and 11) were also affected by the flaring episode and we thus decided not to use these three observations for our analyses based on the He i triplet and Pa$\beta$ and Br$\gamma$ lines. We also identified a feature in the red wing of Pa$\beta$ (at $+120$ $\mathrm{km}\,\mathrm{s}^{-1}$), likely tracing Ti, Ca and Fe lines blending with Pa$\beta$, that does not vary more than the continuum around this line and is thereby expected not to affect significantly our analyses. We then compute the amount by which the equivalent width of these lines vary as a result of activity, which we call activity ‘equivalent width variations’ (EWVs). In the stellar rest frame, we divided each telluric-corrected Stokes $I$ spectrum by the median spectrum shown in Fig. 21, yielding the median-divided spectra in Fig. 22. The activity EWVs are then defined as the EW of these median-divided spectra counted as negative when absorption is larger than average. The values of the EWVs (and the corresponding error bars) were obtained through a gaussian fit to the median-divided spectra, using a gaussian of full-width-at-half maximum equal to 130 $\mathrm{km}\,\mathrm{s}^{-1}$ centred on the stellar rest frame (consistent with the median He i, Pa$\beta$ and Br$\gamma$ median profiles). We note that an activity indicator equal to 0 at a specific epoch does not indicate a lack of detection but rather that the corresponding profile is identical to the median one. The EWVs are provided in Table 1.

For each line, we assumed equal error bars for all spectral points of all observations, which we set to the dispersion between spectra in the continuum about each line (equal to 0.014 for He i and Br$\gamma$, and 0.008 for Pa$\beta$, and tracing mostly photon noise). The corresponding error bars we derive for the EWVs are equal to 2, 1 and 4 pm, respectively. As the integrated flux in the He i and Pa$\beta$ lines is variable at a higher level than that expected from photon noise, we empirically derived the error bars on the EWVs of the three lines using the same method as in Sec. 2.3 (likely overestimating the uncertainties) to account for the intrinsic variability and other main sources of noise that cannot be easily quantified. We achieved this by fitting the He i and Pa$\beta$ EWVs with a sine curve (including 2 harmonics for He i), whereas the Br$\gamma$ EWVs (showing essentially no variation with time) were fitted with a constant, yielding error bars of 25, 3 and 6 pm for the He i, Pa$\beta$ and Br$\gamma$ lines, respectively. To assess the significance of our models, we computed the $\chi^{2}_{r}$ when fitting a constant instead of a periodic curve for the He i ($\chi^{2}_{r}$$=13.6$) and Pa$\beta$ ($\chi^{2}_{r}$$=3.5$) lines, yielding probabilities of 0 and $3.6\,10^{-6}$, respectively, for the detected modulation to be spurious by chance. Even with these simple models and pessimistic estimates of the error bars, we detected a significant modulation of both the He i and Pa$\beta$ EWVs. We however caution that the false alarm probabilies (FAPs) we quote, assuming white noise, may be underestimated if correlated noise dominates, even though pessimistic error bars were used.

The activity EWVs reveal enhanced absorption in phase range 0.4-0.6 for both the He i triplet and the Pa$\beta$ line. This feature is also seen in the dynamic spectrum of the He i triplet (Fig. 23).

From the median-divided spectra, we computed autocorrelation matrices for each of the three lines, considering the line relative intensities within an interval of $\pm 200$ $\mathrm{km}\,\mathrm{s}^{-1}$. The coefficient $C_{ij}$ between velocity bins $i$ and $j$ is defined as:

where $r_{ij}$ is the Pearson linear coefficient between the two velocity bins, $\sigma_{i}$ and $\sigma_{j}$ are the standard deviation in the velocity bins $i$ and $j$, respectively. This definition of the unnormalised coefficient allows us to estimate the relative importance of the correlations, as a high value of $r_{ij}$ associated with a high level of variability is better emphasized than a high value of $r_{ij}$ associated with a low level of variability.

The autocorrelation matrix of the He i triplet reveals that the entire profile correlates well with itself (left panel of Fig. 16). As a strong correlation indicates a common origin for the components, the observed correlation suggests that the entire He i triplet emerges from a single region, likely the stellar chromosphere. The autocorrelation matrix of Pa$\beta$ (middle panel of Fig. 16) shows a correlation / anticorrelation chessboard pattern above the noise level suggesting that the line width is slightly changing with time, getting narrower and deeper at times, then broader and shallower at some other times. In addition, both matrices show that the variability is asymmetric, being larger in the blue wing of these lines, possibly indicating the presence of a stellar wind. No particular pattern is apparent in the autocorrelation of the Br$\gamma$ line (right panel of Fig. 16). We also show the normalized autocorrelation matrices (i.e. the $r_{ij}$ coefficients) in Fig. 24.

For each velocity bin of the median-divided spectra, we computed a Generalized Lomb-Scargle periodogram as introduced by Zechmeister & Kürster (2009) thanks to the PyAstronomy PYTHON module (Czesla et al., 2019). We show them as 2D maps in Fig. 17 and 25. Using the typical error bars derived in Sec. 4.2 for all spectral points, we found that the He i and Pa$\beta$ profiles exhibit rotational modulation (with aliases associated with the observing window), with a stronger variability in the blue wing (as seen in Sec. 4.3). These results are consistent with those obtained from the EWVs in Sec. 4.2, though less obvious as information is not integrated over the line profile. The Br$\gamma$ periodogram does not show any clear period, consistent with EWV$\mathord{\sim}0$.

For the first time, we jointly used high-resolution spectropolarimetry and high-precision photometry in ZDI to reconstruct the brightness distribution at the surface of the star. Both data are complementary as spectropolarimetry mostly constrains the location of the spots (from the distorsions of profiles) while photometry mainly informs on their contrast relative to the quiet photosphere.

Taking into account photometry yields a higher number of brightness features and higher contrasts in the ZDI image, especially at low latitudes (Fig. 4). These features are needed so that both SPIRou and TESS data can be fitted at unit $\chi^{2}_{r}$  with the light curve fitted down to 1.6 mmag RMS (with residual correlated noise likely attributable to small-scale rapidly-evolving surface brightness features that cannot be properly reproduced by ZDI). These latitudes are usually not well reconstructed when considering spectropolarimetry only, as ZDI is best sensitive to features located at higher latitudes, in the visible hemisphere (Vogt et al., 1987; Brown et al., 1991). We suspect that this increase in spot coverage reflects that TESS and SPIRou do not see the same spot distributions because of the difference of spectral domains. To further improve the accuracy of the brightness modeling at the surface of the star, one would need to secure photometric data in the $JHK$ bands that would provide an ideal match to the SPIRou spectropolarimetric data. This would ensure in particular that all brightness features present at the surface of the star (including those in the polar regions) affect spectroscopic and photometric data in the same way (which is not the case for, e.g., the cool polar spot detected in the optical but not in the NIR).

Ground-based photometry allowed us to perform a similar analysis as Yu et al. (2019) (Fig. B7 of their paper) but with $V-R_{\rm c}$ and $V-I_{\rm c}$ colour indexes. We fitted these indexes as a function of the magnitude in the $V$ band using a simplistic two-temperature model based on colour indexes from Bessell et al. (1998). Our model features a fixed temperature of 4500 K for the photosphere, a surface gravity $\log{g}=4.0$, and a fixed temperature for spots with varying filling factor. We found an optimal spot temperature of 3750 K with a typical surface spot coverage of about 70% (Fig. 26), consistent with previous photometric measurements (Yu et al., 2019). This rather high spottedness level (consistent with that found for the similar wTTS LkCa 4, Gully-Santiago et al. 2017) suggests in particular that a large fraction of the stellar surface is more or less evenly covered with small features that are not accessible to (and thus not reconstructed by) ZDI.

V410 Tau is a wTTS that has been extensively studied in the past, mainly in the optical. Our study is innovative as we used NIR observations to constrain the brightness and the magnetic field of the star.

As expected, NIR leads to a less contrasted surface brightness map. We find a reasonable statistical agreement for low- and mid-latitude spots (but less so in the polar regions) between optical and NIR maps, even though secured at different epochs. However, the absence of a polar spot in the map reconstructed from NIR data is surprising since such a feature consistently showed up in images derived from optical data up to now (Joncour et al., 1994; Hatzes, 1995; Skelly et al., 2010; Rice et al., 2011; Carroll et al., 2012; Yu et al., 2019). Our reconstructed map obtained with a mask containing only molecular lines, more sensitive to cooler regions, further confirms that no polar spot is detected at NIR wavelengths. Although we cannot entirely dismiss it, the option that the polar spot disappeared at the time of our SPIRou observations seems unlikely given the persistent presence of this feature in all previously published studies. The fact that the TESS light curve yields an average rotation period that is consistent with the trend derived from previous photometric data by Yu et al. (2019) suggests that the spot configuration at the surface of V410 Tau did not drastically evolve since 2016 and in particular that the cool spot reconstructed near the pole from optical data was likely still present in 2019.

If the dark polar feature systematically seen at optical wavelengths is indeed not visible in the NIR, it suggests that continuum opacity above polar regions of V410 Tau is much larger in the optical domain than in the NIR for some reason. A speculative option, to be investigated further, may be that dust grains, such as those present in the upper solar atmosphere though in larger concentrations, tend to cluster in polar regions of the upper atmosphere of V410 Tau, making them appear much darker at optical than at NIR wavelengths.

Applying ZDI to our Stokes $I$ and $V$ LSD profiles simultaneously allowed us to reconstruct the large scale magnetic topology of the star. Our results are consistent with previous studies (Skelly et al., 2010; Yu et al., 2019). We find that the large-scale magnetic field has an average surface strength of about 410 G and that the radial field can be more intense locally, reaching up to 1.1 kG. Although V410 Tau is still fully convective, the magnetic field presents a strong toroidal component of unclear origin, as for the other fully-convective wTTS LkCa4 (Donati et al., 2014). More observations of fully-convective wTTSs are thus clearly needed to further constrain the origin of this strong toroidal field.

We also found that the poloidal component contributes to nearly 60% of the overall magnetic energy, compatible with the recent measurements derived from NARVAL (at the Telescope Bernard Lyot) optical data in 2016 (Yu et al., 2019). In addition, the polar strength of the dipole component (of the poloidal field) is close to 400 G, which again supports the reported increase in the intensity of the dipole from 2008 (Yu et al., 2019). These properties are compatible with those obtained by Yu et al. (2019) from their 2016 data set, and more generally with the long-term evolution they pointed out. These results suggest that, if a magnetic cycle exists, it is likely longer than 11 years. More observations of V410 Tau would be needed to confirm whether the observed tendency reflects part of a magnetic cycle as suggested by other studies (Stelzer et al., 2003; Hambálek et al., 2019) or rather intrinsic variability of a stochastic nature.

The longitudinal field as derived from SPIRou data shows similar fluctuations than that from optical studies (of period $\mathord{\sim}P_{\rm rot}$), but with error bars that are about 1.7 times smaller (typically 30 G) in half the exposure time, clearly demonstrating the benefits of studying magnetic fields of young stars in the NIR thanks to the enhanced Zeeman effect.

We constrained the surface differential rotation of V410 Tau with ZDI from our Stokes $I$ and $V$ LSD profiles separately, both results being compatible within 1.5$\sigma$. Our estimates are also consistent with those provided by Yu et al. (2019) within $\mathord{\sim}3\sigma$, although our Stokes $I$ LSD profiles yielded slightly lower value. We note that the error bars obtained from Stokes $I$ LSD profiles are larger in the NIR than in the optical, which is likely due to the lower number of observations but also to the lower contrast of the brightness features. For Stokes V data, the differential rotation parameters we derived are similar to the optical measurements of Yu et al (2019), with error bars of the same magnitude despite the sparser data thanks to the enhanced Zeeman effect in the NIR.

We note that our estimates of differential rotation are larger than those derived by Siwak et al. 2011 from photometric data collected with the MOST space-telescope in 2009. This photometric measurement of differential rotation is also inconsistent with the estimates of Yu et al. (2019), despite having been collected at a close-by epoch. We thus suspect that this difference is related to the two-spot model used by Siwak et al. (2011) known to be inappropriate for stars like V410 Tau given the complex spot distributions reconstructed with ZDI (featuring both bright and dark spots).

The He i triplet at 1083 nm, the Pa$\beta$ and Br$\gamma$ lines are used as proxies to study the chromospheric activity of V410 Tau. A flare was detected at cycle 21.752 which also affected the two subsequent observations. Our analyses reveal that both the He i and the Pa$\beta$ lines are rotationally modulated, while no significant variations are observed in Br$\gamma$.

To obtain a rough description of the large-scale stellar magnetosphere, we extrapolated our magnetic image at the surface of V410 Tau into 3D maps, assuming that the magnetic field is potential (following the method described by Jardine et al. 1999) and that the source surface at which field lines open is located at 2.1 R_⋆, following Yu et al. (2019) (Fig. 18). We see that enhanced absorption in chromospheric lines, occuring in phase range 0.4-0.6, takes place slightly before the magnetic pole crosses the line of sight (at phase 0.7, see Fig. 18), i.e., when one may have expected it to occur by analogy with solar coronal holes (darker in regions of open field lines). This phase lag may relate to the potential field assumption being no more than a rough approximation in our case. The reconstructed large-scale magnetic field indeed features a strong toroidal component (with intense azimuthal fields located close to the open field line region at phase 0.7, see Fig. 10) that may suggest that the large-scale surface field is significantly stressed at these phases. Another option is that this enhanced absorption episode is due to the presence of massive prominences trapped in closed coronal loops (such as those reported in Yu et al. 2019) and crossing the stellar disc at phases 0.4-0.6.

Obviously, the way the He i triplet and the Pa$\beta$ line behave in wTTSs, and in particular how the He i and Pa$\beta$ fluxes respond to the topology of the large-scale field remains to be investigated in more details. This will be the subject of forthcoming papers.

RVs of V410 Tau derived from Stokes $I$ LSD profiles exhibit a full amplitude of about 4.5 $\mathrm{km}\,\mathrm{s}^{-1}$ and a dispersion of 1.40 $\mathrm{km}\,\mathrm{s}^{-1}$ RMS. These values are smaller than those generally observed in the optical at roughly the same SNRs, with amplitudes ranging from 4 to 8.5 $\mathrm{km}\,\mathrm{s}^{-1}$ and a typical dispersion of 1.8 $\mathrm{km}\,\mathrm{s}^{-1}$ RMS, respectively (Yu et al., 2019). This confirms the gain in using NIR observations to reduce the activity jitter in RV measurements. We find that the amplitude of RV jitter is reduced by up to a factor of 2, consistent with results of previous optical and NIR RV studies of TTSs (Prato et al., 2008; Mahmud et al., 2011; Crockett et al., 2012).

For each of our ZDI reconstructed brightness maps (Fig. 4) we computed the RV curve that results from the brightness features at the surface of V410 Tau. The ZDI image taking into account SPIRou data only, corresponding to a static brightness distribution, yields filtered RVs with a dispersion of 0.13 $\mathrm{km}\,\mathrm{s}^{-1}$ RMS, i.e. about 25% lower than in the optical for this star (typically 0.17 $\mathrm{km}\,\mathrm{s}^{-1}$; Yu et al. 2019), which suggests that the evolution of spots is not significant over our observations. We note that adding TESS data to the SPIRou data in the ZDI modeling does not improve, and actually even degrades, the accuracy of the filtering process (even after explicitly taking into account the difference in brightness contrasts at SPIRou and TESS wavelengths in the imaging process). This result demonstrates that the brightness distributions as seen by SPIRou and TESS are genuinely different and cannot be simply scaled up from one another, e.g., using Planck’s law, with some features showing up in one spectral domain but not in the other (like the prominent polar spot detected ine the optical but not seen at NIR wavelengths). Filtering the activity jitter would thus likely be more efficient with ZDI applied to data sets combining SPIRou data with high-precision NIR photometry. We also modeled the RV activity jitter using GPR yielding a dispersion of filtered RVs about twice smaller than with ZDI models thanks to the higher flexibility of GPR to model intrinsic variability in the periodic modulation of the RV curve, that results from the evolution of the spot configuration at the surface of the star.

The periodograms of RVs (Fig. 27) do not show any periodic signature beyond that from V410 Tau A, which further confirms that our spectropolarimetric data mainly probe the primary star, and not (or no more than very marginally) its 2 companions. In addition, our filtered RVs show no evidence for a RV signal from a potential giant planet on a close-in orbit (Fig. 27), consistent with previous observations that did not suggest the presence of a hJ (Yu et al., 2019). To derive an upper limit on the mass of a potential planet from our data, we proceeded as in Yu et al. (2019) and applied GPR on simulated datasets (with the same temporal sampling as that achieved for our 2019 observations) featuring both a RV activity jitter (computed from the results of Sec. 4.1) and a RV signal from a planet on a circular orbit with a white noise identical to that of our measurement (of 181 $\mathrm{m}\,\mathrm{s}^{-1}$ RMS). For each simulation, we compared models including both the planet and the activity jitter, with those including only the activity jitter, to assess the significance level at which a close-in giant planet of given mass could be detected from our data. From the difference of logarithmic marginal likelihood between both models (detection threshold set at $\Delta\mathcal{L}=10$), we found that, for a planet-star separation lower than 0.09 au, only planets with a mass larger than $\mathord{\sim}5$ M_jup can be reliably detected (at a >3$\sigma$ level), consistent with the (more stringent) upper limit (of $\mathord{\sim}1$ M_jup) derived by Yu et al. (2019) for the same planet-star separation.

Detecting close-in massive planets typically requires carrying out monitorings over several months during which the surface of the star can evolve significantly. This intrinsic variability cannot be modeled with the current version of ZDI that assumes a static distribution of features at the surface of the star (except for differential rotation), often forcing one to split data sets into smaller subsets that can be modelled independently from one another (e.g. Donati et al. 2017; Yu et al. 2017, 2019). In order to get a more global and consistent description of the stellar surface activity over several months, one needs to be able to model at the same time both the distribution of surface (brightness and magnetic) features and its evolution with time using all data at once. In this aim, we started to modify the original ZDI code to couple it with GPR, in order to simultaneously benefit from the physical modeling provided by ZDI (to detect and characterize stellar surface features), and from the flexibility provided by GPR (to describe the temporal evolution of these features). This new version of ZDI is currently under development and will be the object of forthcoming publications.

Our study illustrates the benefits of NIR (vs optical) observations with instruments like SPIRou, to investigate the magnetic topologies of young stars and look for the potential presence of hJ on close-in orbits through RV measurements. New monitorings of V410 Tau will provide strong constraints on the existence of a magnetic cycle (and the underlying dynamo processes), will bring further clues on the enigmatic strong toroidal field that the star is able to trigger despite being fully convective, and will allow us identifying the main differences between images reconstructed from optical and NIR data, especially in the polar regions. More generally, SPIRou observations of PMS stars, including those carried out within the SLS, will offer the opportunity to investigate in more details the impact of magnetic fields on star / planet formation, and in particular to accurately characterize young planetary systems hosting transiting planets such as AU Mic and V1298 Tau (David et al., 2019a, b; Plavchan et al., 2020; Klein et al., 2021), allowing one to refine the mass-radius relation of planets at an early stage of evolution.