© The Authors 2024

1. Introduction

The accretion processes of pre-main sequence (PMS) stars induce several forms of variability in stellar spectra. The classical T Tauri stars (CTTS) are low-mass PMS objects surrounded by an accretion disc and present a magnetically driven accretion process. The strong magnetic field of these objects exerts magnetic pressure on the disc, forcing the material to leave the disc plane and fall onto the star along the magnetic field lines, a process known as magnetospheric accretion (see, e.g., Bouvier et al. 2007; Hartmann et al. 2016 as reviews). The material thus reaches the stellar surface in a very localised region near to the dipole magnetic pole. The material’s kinetic energy is dissipated in a shock, heating the region and producing a bright spot, also called a hot or accretion spot. The signature of this spot on the stellar photospheric lines is the same as the cold spots, that is, the darker magnetic regions seen on the Sun, which are also present on CTTSs: a bump in the line profile, blueshifted and redshifted with the stellar rotation, and producing an apparent radial velocity (RV) modulation with the same period as the rotation period of the object (Vogt & Penrod 1983). This feature is a well-known signature of stellar activity and may hide other modulations of the RV induced by a planetary or stellar companion.

The aim of the present work is to study the RV of HQ Tau (right ascension (RA): 4h35m, declination (Dec): +22°50′), an intermediate-mass T Tauri star (IMTTS) of which the magnetospheric accretion process was reported and studied in detail by Pouilly et al. (2020) (hereafter Paper I). These latter authors reported a modulation of the RV on the star’s rotation period and ascribed this to a large polar spot. However, the mean RV (⟨V_r⟩ = 7.22 ± 0.27 km s⁻¹) was far below the median velocity of the Taurus region (about 15.5 km s⁻¹, Galli et al. 2019), and was found to be inconsistent with the work of Nguyen et al. (2012) and Pascucci et al. (2015) (16.3 ± 0.02 and 17.2 ± 0.2 km s⁻¹, respectively). In addition to the erroneous renormalised unit weight error (RUWE) measured by Gaia (RUWE = 13.39; Gaia Collaboration 2023), we thus suspected the presence of a companion, and the additional RV measurements obtained in late 2019 and early 2020 confirmed a larger variation of the RV than that produced by the stellar activity only (see Table 3 of Paper I).

HQ Tau was first suspected to be a tight binary by Simon et al. (1987) based on lunar occultation, with a separation of 4.9 ± 0.4 mas. This estimation was later revised by Chen et al. (1990) to 9.0 ± 0.2 mas, who noted that the binary nature of this system is not obvious. All subsequent studies (Richichi et al. 1994; Simon et al. 1995, 1996; Mason 1996) led to the same conclusion, that HQ Tau is a single object, until our work in Paper I. Indeed, the large RV variation, in addition to the inner cavity in HQ Tau’s disc suspected by Long et al. (2019) and Akeson et al. (2019) from ALMA observations, are in favour of the presence of a companion, motivating further investigations of this aspect of this system.

The new dataset presented in this work covers a six-month period, which should allow us to probe the presence of a companion and derive the orbital elements of the system given the large RV variation that occurs over about 2 months (between our measurements in 2019 and 2020), meaning an orbital period of about 4 months if this variation corresponds to the peak-to-peak amplitude of the velocity modulation. The spectroscopic dataset is described in Sect. 2, and the entire analysis is presented in Sect. 3 and discussed in Sect. 4. We present our conclusions from this work in Sect. 5.

2. Observations

In this section, we describe the dataset used to perform our analysis. We used four different instruments: the Spectropolarimètre Infrarouge (SPIRou, Donati et al. 2020) and the Echelle SpectroPolarimetric Device for the Observation of Stars (ESPaDOnS, (Donati 2003)) at Canada-France-Hawaii Telescope (CFHT), Neo-Narval at Telescope Bernard Lyot (TBL), and SOPHIE at Observatoire de Haute Provence (OHP). The 39 observations were performed in late 2020 and early 2021, covering almost 6 months in total, with a sampling cadence that varies between less than 1 day and 28 days. The complete log of observations is provided in Appendix A.

2.1. SPIRou

Most of the observations taken at the CFHT were made using SPIRou, a near-infrared spectropolarimeter covering the YJHK band in a single shot at a resolution of about 75 000. The ten observations were taken between 2020 September 4 and 2021 January 5 using the polarimetric mode, meaning that they are composed of four subexposures in different polarimeter configurations, allowing us to derive the intensity (Stokes I), and the circularly polarised (Stokes V) and null spectra. The exposure time of each subexposure is about 150 s. The observations were reduced and telluric corrected using the APERO data reduction system (Cook et al. 2022), and reached a signal-to-noise ratio (S/N) ranging between 73 and 90 for the subexposure combination at the order #47.

2.2. ESPaDOnS

Four observations were performed using ESPaDOnS mounted at the CFHT between 2020 November 30 and 2020 December 8. ESPaDOnS covers the 300–1000 nm wavelength range and reaches a resolving power of 68 000. Here again, observations were made using the polarimetric mode, with a ∼150 s exposure time for each subexposure. These observations were reduced using the Libre-ESpRIT package (Donati et al. 1997), and reach a S/N of between 73 and 90 for the subexposure combination on Stokes I at 730 nm. Unfortunately, the observation on 2020 November 30, which has the highest S/N, seems to be contaminated by the Moon and is not used in this work.

2.3. Neo-Narval

We obtained 11 observations at the TBL using the ESPaDOnS’s twin Neo-Narval (López Ariste et al. 2022) between 2020 September 5 and 2021 February 23. These data were reduced using the latest DRS version of the instrument. The polarimetric mode was used here again, with approximately 450 s of exposure time for each subexposure. The resolution of this instrument is about 65 000 and it covers the 380–1000 nm wavelength range. The S/N of the intensity spectrum ranges from 7 to 20 at 600 nm.

2.4. SOPHIE

The largest individual data set used in this work consists of 14 observations carried out with the SOPHIE spectrograph mounted at the OHP between 2020 September 3 and 2021 February 17. The observations were taken in a single exposure of ∼1800 s. This instrument has a resolution of 40 000 and covers the 390–700 nm wavelength range. The raw data were reduced using the SOPHIE real-time pipeline (Bouchy et al. 2009), and the resulting spectra reach a S/N of between 40 and 64 at 600 nm. One of the observations (2020 November 27) was affected by Moon contamination and is therefore not used in this work.

3. Results

3.1. Radial velocity

The crucial first step of this analysis is to derive the RV. As described in Paper I, the amplitude of the RV modulation seems to be of around 14 km s⁻¹, including a spot modulation of about 6 km s⁻¹. We are thus looking for an orbital modulation of the same order as the spot modulation. We used the cross-correlation method described in Paper I to derive the RV using the same photospheric template (i.e. Melotte25-151/V1362 Ori) for the observations in the visible frame, meaning the SOPHIE, Neo-Narval, and ESPaDOnS’ spectra, and V819 Tau for the SPIRou observations. We computed the RV on ten wavelength windows of about 10 nm in width between 480 and 880 nm for each visible observation and on five windows, also of about 10 nm in width, between 1550 and 2150 nm for the SPIRou observations. Then, for each observation, we averaged the results of each window and computed a standard deviation to derive the RV and its uncertainty. All the computed values are provided in Table A.5 and the V_r curve is shown in Fig. 1.

Fig. 1. Radial velocity curve (top) and corresponding Lomb-Scargle periodogram analysis (bottom) of HQ Tau over the 20B semester. The colors of the radial velocity curve correpond to the different instruments used: NeoNarval at TBL (purple), SOPHIE at OHP (pink), ESPaDOnS (orange), and SPIRou (brown) at CFHT. On the bottom panel the Lomb-Scargle periodogram is shown in black and the FAP in grey.

Although the values are relatively widely scattered, they follow a trend in that they oscillate around the mean of V_r ∼ 17 km s⁻¹ on a period of about 100 days. Furthermore, the amplitude of the variation (about 15 km s⁻¹) and the mean velocity are far larger than those measured in Paper I (∼4 km s⁻¹ and 7.22 km s⁻¹, respectively), while the photometric variations were within 0.14 mag in V band at the time this new dataset was obtained according to the American Association of Variable Star Observers (AAVSO) survey. This indicates that another phenomenon modulates the RV in addition to a surface spot.

The Lomb-Scargle periodogram (Scargle 1982) shown in Fig. 1 yields a detection at f ∼ 0.4 d⁻¹, with a false-alarm probability (FAP) of 10⁻³ (computed following the approximation of Baluev (2008)), meaning P ∼ 2.5 d, which is consistent with a stellar rotation period, as expected. Furthermore, we note a secondary peak at f ∼ 0.009 d⁻¹, meaning P ∼ 111 d, which is consistent with the apparent modulation seen on the RV curve; however, the FAP of 0.48 for this signal is too high for it to be considered as an independent detection of orbital motion.

3.2. Least-squares deconvolution profiles and bisector inverse slope

The analysis of the bisector is commonly used to identify the origin of a RV variation. Indeed, when a companion shifts the line velocity, the stellar activity alters the line shape, producing an asymmetric line profile, which can be studied using the bisector inverse slope (BIS; Queloz et al. 2001).

Here, we study the BIS of the least-squares deconvolution (LSD; Donati et al. 1997) profiles of the unpolarised signal (Stokes I). The S/N increase provided by this method, as well as the switch to the velocity space, allow us to take into account the observations with low S/N, as well as their different wavelength regions. In order to normalise the LSD weights, we used an intrinsic depth of 0.2, a mean Landé factor of 1.2, and a mean wavelength of 520 (1200) nm for the optical (infrared) spectra. For the computation, we used a mask derived from a line list provided by the VALD (Ryabchikova et al. 2015) database covering the wavelength range of each instrument and from which we excluded the regions contaminated by emission lines or telluric absorption. The resulting Stokes I profiles used in this study are shown in Fig. 2.

Fig. 2. Measurement of the BIS of the Stokes I profiles. The HJD is indicated at the top of each plot. The black lines are the Stokes I profiles and the grey dotted line delimits the top and bottom regions selected for the BIS computation. The blue points are the bisector in these regions and the red dots are the mean bisector of each region used to compute the BIS.

To compute the BIS, we selected two regions at the top and the bottom of the line, in which we computed the mean velocity of the bisector. The BIS is then defined by the difference between the top velocity and the bottom one. The two regions were selected as follows: we ignored the first 15% of the profile containing the continuum and the wings of the line, as well as the last 15% containing the bottom of the profile, where the bisector cannot be computed because of the noise or because of an excessively strong activity signature splitting the profile into two parts in this region. The top region is thus defined as 25% of the profile from the upper limit, and the bottom region is defined as 25% of the line from the lower limit. The error on the BIS is defined as the quadratic mean of the standard deviation of the bisector of the top and bottom regions. We removed five observations¹ from this analysis because the S/N did not allow us to properly derive the bisector of the profile. The computation of the BIS is illustrated in Fig. 2, and the values are provided in Table A.5.

The BIS varies with a period of 2.42 days (f = 0.41 d⁻¹, FAP = 0.09; see Fig. 3), which is consistent with the rotation period of the star and thus confirms that it is tracing the stellar activity. We do not observed any signal around 111 d. Figure 4 shows the BIS versus the RV. In the case of a RV that is fully modulated by the stellar activity, a linear relation with a slope of −1 is expected, and a BIS of 0 km s⁻¹ is expected at the star’s real velocity as well. For HQ Tau, a linear downward trend is observed, as well as a clustering of measurements around BIS = 0 km s⁻¹ and V_r ∼ 17 km s⁻¹ (the mean velocity of our data set). However, the fitted slope is −0.50 ± 0.09, and the measurements show a larger dispersion at the lower and higher velocities. This indicates that, in addition to the stellar activity, another phenomenon is affecting the RV modulation.

Fig. 3. Lomb-Scargle periodogram of the BIS. The peak at f = 0.41 d−1 corresponds to a period of 2.42 d and has a FAP of 0.09.

Fig. 4. BIS versus RV. The colour scale corresponds to the HJD, the black line is the linear fit, and the grey-shaded area is the 1σ uncertainty on the slope of the fit.

3.3. Multi-dimensional Gaussian process regression

Although visual inspection of the RV curve and the Lomb-Scargle periodogram analysis seems to indicate the presence of a second modulation with a lower frequency and lower amplitude than the stellar activity, the signal remains insufficiently strong to constrain any orbital parameters. We therefore used Gaussian process regression (GPR) to model the stellar activity signal with pyaneti, a tool developed and fully described by Barragán et al. (2019); see also Barragán et al. (2022) for the latest version. Here we simply describe its basic principles. This method, which is often used in searches for exoplanets around active hosts, considers that the modulation can be described as a finite set of random variables related by a multi-variate normal distribution, defining a Gaussian process (GP). According to Tracey & Wolpert (2018), this distribution is given by:

$$\begin{array}{c}\hfill P(\mathbf{t})=\frac{1}{\sqrt{{(2\pi )}^N|\mathbf{K}|}}\exp [-\frac{1}{2}{(\mathbf{t}-\mathit{\mu })}^T{K}^{-1}(\mathbf{t}-\mathit{\mu })],\end{array}$$(1)

where t = t_{i, (i = 1, ..., N)} are the locations of the N variables, μ are the mean values, and K is the covariance matrix (or kernel), giving the relation between the variables.

The aim of GPR is, assuming that the time series is a sample of GP, to optimise the parameters ϕ of the mean function μ(t; ϕ), and the hyper-parameters Φ of the kernel function γ(t_i, t_j : Φ). In our case, μ is physically motivated by a Keplerian, and γ is a composition of the ‘white-noise’ (WN) and ‘quasi-periodic’ (QP) kernels, to describe the rotational modulation and the error of our data. These kernels are given by

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{WN}}(t_i,t_j)={\sigma }_i^2{\delta }_{\mathit{ij}},\end{array}$$(2)

where σ_i is the error of the value i, and δ_ij is the Kroenecker delta, and

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{QP}}(t_i,t_j)=A^2\exp \{-\frac{{sin}^2[\pi (t_i-t_j)/P_{\mathrm{GP}}]}{2{\lambda }_p^2}-\frac{{(t_i-t_j)}^2}{2{\lambda }_e^2}\},\end{array}$$(3)

where A is the amplitude scaling the deviation from the mean function, P_GP is the period of the GP (in our case, the rotation period), λ_p is the inverse of the harmonic complexity (i.e. the complexity of the variation within a period), and λ_e is the long-term evolution timescale.

In this work, we used the multi-dimensional GPR implemented in pyaneti, meaning that we used two time series simultaneously: the RVs and an activity indicator, in our case the BIS of the LSD profiles. These two sets of variables are described by a composition of the function drawn by the GP, G(t), and its derivative $\dot{G}(t)$ as follows:

$$\begin{array}{c}\hfill \{\begin{array}{c}\mathrm{\Delta }\mathrm{RV}(t)=A_{1}G(t)+B_1\dot{G}(t)\hfill \\ \mathrm{BIS}(t)=A_{2}G(t)+B_2\dot{G}(t),\hfill \end{array}\end{array}$$(4)

where ΔRV is the RV modulation due to both the stellar activity and the companion. The RV is described by RV(t) = RV_syst + ΔRV, where RV_syst is the systemic RV.

Physically, the function G(t) represents the projected area of the active regions on the visible stellar surface that induces the RV (and other activity indicator) variations, as a function of time. The RVs and the BIS are affected by these regions, and by their time evolution (Dumusque et al. 2014); the best description of these data is thus given as a function of G(t) and its time derivative $\dot{G}(t)$.

The optimisation of the complete set of parameters is done through a Markov Chain Monte Carlo (MCMC) procedure using 500 chains to sample the parameter space. To perform the simultaneous fit of RV and BIS, we need both values for each observation; we therefore had to exclude the RV measurements for which we could not measure the BIS from the fit (see Sect. 3.2). This analysis yielded the detection of an orbital modulation in the RVs with a semi-amplitude of $K=1.{25}_{-0.71}^{+0.83}$ km s⁻¹ and a period of $P_{\mathrm{orb}}=125.8_{-55.2}^{+31.5}$ d, in addition to a RV modulation by stellar activity consistent with the period of Paper I, and a systemic velocity of RV${}_{\mathrm{syst}}=17.5_{-1.0}^{+0.9}$ km s⁻¹, which is consistent with the median velocity of the Taurus cloud (v = 16.4 ± 1.1 km s⁻¹, Galli et al. 2019). Using the binary mass function and the Kepler’s third law, this means that a companion, a brown dwarf with a minimal mass of 50.1${}_{-32.0}^{+29.6}{M}_{\mathrm{jup}}$, is orbiting HQ Tau with a semi-major axis of about 0.61 AU. The complete set of parameters inferred by this analysis and the prior used are provided in Table 1. In addition, the posterior distribution of the orbital period and semi-amplitude is shown in Fig. 5, and the resulting curves for RV and BIS are shown in Fig. 6. We did not succeed in constraining the orbit’s eccentricity and the argument at periastron (we obtained flat posterior distributions), and therefore we assume a circular orbit and ω = 0. The hyperparameters of the GP are consistent with expectations, meaning that the GP’s period recovers the stellar rotation period, the A₁ and A₂ coefficients are of opposite sign, indicating a phase opposition of the RV and BIS modulations, and λ_e (representing the spot’s lifetime) is larger than the rotation period, indicating that this activity region can be used to recover the rotation period.

Fig. 5. Posterior distribution of the orbital period (left) and of the semi-amplitude (right). The vertical red line shows the mean value and the grey-shaded area illustrates the 68.8% confidence level.

Fig. 6. Fit of the RV (top) and BIS (bottom) curves with pyaneti. The blue points are the measurements with their uncertainties. On the top panel the black curve is the GP fit of the stellar activity (purple) plus the Keplerian modulation by the companion (red). The grey-shaded area is the 1σ uncertainty of the fit. On the bottom panel the black curve is showing the GP fit of the stellar activity only with its 1σ uncertainty in grey.

In order to decipher whether or not the stellar activity alone can explain the observed variability, we performed the same analysis, but this removing the Keplerian modulation. The results are presented in Appendix B. We then compared the Akaike information criterion (AIC; Akaike 1974) in its corrected prescription (AICc; Sugiura 1978) for a small number of data points; the AICc estimates the quality of a given model compared to another one through an estimate of the amount of information lost by the model. The AICc obtained for the model including the Keplerian modulation is 304.77, while that obtained without the Keplerian modulation is 322.07, indicating that the model that includes the Keplerian modulation is preferred. Finally, we performed a leave-one-out cross validation (LOOCV) test to assess the robustness of this solution. This test consists of performing the analysis including all data points except one and computing the mean squared error (MSE) of the model at this point. This is done as many times as the number of data points, excluding a different measure each time. The final MSE is then derived as the average of all the individual MSEs computed. The MSEs obtained for our solution are 6.96 and 3.41 km s⁻¹ for the RV and the BIS, respectively. The corresponding values for the model assuming only the stellar activity are 7.16 and 3.62 km s⁻¹, confirming that the Keplerian needs to be included to explain the observed variability. Together with the MSE, we obtain the Keplerian parameters, in particular T0 and K, which yield HJD 2 459 081.42 ± 4.83 and 1.35 ± 0.20 km s⁻¹, respectively. These values are perfectly consistent with the parameters of Table 1, demonstrating the stability and consistency of this solution.

4. Discussion

HQ Tau was selected for this semester-long RV follow-up because of the surprising intermediate results of Paper I on a two-week high-resolution spectropolarimetric 2017 time series. In order to study the activity of this object, we derived the RV curve, from which we draw the conclusion that a surface spot is modulating the RV on the stellar period. However, the mean RV is not consistent with previous measurements in the literature, and the curve shows a slight downward trend over the time series. Two additional snapshot observations taken in late 2019 and early 2020 reveal a larger velocity modulation around the already published values, suggesting an additional source of RV modulation superimposed on the activity signal, presumably a companion.

To disentangle the RV signal of the activity from the modulation induced by the suspected companion, we performed a multi-variate GPR using the Python tool pyaneti (Barragán et al. 2019), reproducing the activity signal of the BIS and RV modulation with a GP, and adding a Keplerian signal to the latter. The MCMC procedure to recover the Keplerian’s parameters and GP hyperparameters yields a Keplerian period of $P_{\mathrm{orb}}=125.8_{-55.2}^{+31.5}$ d, consistent with the Lomb-Scargle periodogram of the RV curve, and a semi-amplitude of $K=1.{25}_{-0.71}^{+0.83}$ km s⁻¹. Together with the activity-induced modulation reported in Paper I (∼6 km s⁻¹ peak-to-peak), this yields a peak-to-peak amplitude of about 10 km s⁻¹, consistent with the dataset’s peak-to-peak variability (∼13 km s⁻¹). Furthermore, the derived time at periastron and the fitted orbital parameters yield a minimum at about HJD 2 458 103, which is consistent with the downward trend of the RV modulation observed in Paper I (between HJD 2 458 055 and HJD 2 458 067).

Concerning the stellar activity, the period of the GP, which is expected to be the stellar rotation period in our case, is perfectly consistent with our previous results and the Lomb-Scargle periodogram of the present dataset as well, and the λ_e > P_GP confirms that this period is relevant and that the quasi-periodic kernel can be used for this analysis. The harmonic complexity is relatively low (i.e. λ_p ≳ 1), which indicates a moderate to low activity of the star, as expected from the more evolved status of HQ Tau, and the close-to-sinusoidal modulation of its RV (see Paper I).

These orbital parameters, together with the stellar parameters, yield a semi-major axis of 0.61 AU, which is fully consistent with the first estimation of Simon et al. (1987) (4.9 mas separation; i.e. 0.78 AU at the distance of HQ Tau), and is inside the inner disc. The presence of such a companion within the inner disc could possibly disrupt it (en reference therein Long et al. 2018). As discussed in Paper I, this is consistent with the ALMA observation modelled by Long et al. (2019), where the authors highlight dust depletion towards the inner disc region, indicating the presence of a dust cavity, which is nevertheless not well resolved in their data. Assuming that this dust cavity is produced by the companion, this means that the orbit inclination is the disc inclination (i.e. 53.8°), constraining the companion’s mass to about 62 M_jup.

Although the results are consistent with the present dataset, the derived semi-amplitude does not allow us to reach the mean velocity obtained in Paper I. There are several possible explanations for this. Firstly, this six-month follow-up analysis may not be representative, given the three-year gap between the two datasets. The derived λ_e of 7.51 days, although larger than the rotation period as needed for the quasi-period kernel used, remains short regarding the time range of the studied phenomenon. This short spot lifetime indicates that the activity signal might change quickly, especially between 2017 and 2020, where the modulation by the stellar activity might have changed. Long-term RV variations may occur in active stars due to magnetic cycles that change the cool spot configuration at the stellar surface (Grankin et al. 2008), or due to a slowly evolving surface magnetic field that modifies the distribution of surface hot spots (e.g. see the multi-epoch studies of GQ Lup, DN Tau, and V2129 Oph by Donati et al. 2011, 2012, 2013, respectively). However, in all cases, the amplitude of activity-related RV variations is expected to amount to a fraction of the star’s v sin i and, more importantly, the stellar RV averaged over a rotational cycle should vary very little as the spot modulation induces nearly symmetrical RV variations relative to the star’s rest velocity. This is not what is observed for HQ Tau, where the rotationally averaged RV in 2017 strongly departs from that observed in 2020/2021. A second explanation could be the presence of a second companion; with a much longer period, the RV modulation cannot be observed in the measurements analysed here, but could yield significantly lower velocities three years earlier.

Aiming to begin the exploration of these hypotheses, we investigated the multi-variate GPR using pyaneti on the combined 2017 and 2020/2021 data sets. We performed the same analysis as in Sect. 3.3, in two configurations for the two hypotheses, meaning one or two companions modulating the RV curve in addition to the stellar activity, and the results are shown in Appendix C. The results from the one-companion hypothesis appear inconsistent with the velocities reached in the 2017 measurements. The velocity modulation induced by the companion is still about 5 km s⁻¹ above that measured, which is compensated for by a global decrease in the modulation by the stellar activity. Furthermore, the orbital period obtained (∼250d) is inconsistent with the periodogram analysis. However, the systemic velocity derived (RV_syst = 15.5 ± 1.5 km s⁻¹) is consistent with the median velocity of the Taurus cloud (Galli et al. 2019). The two-companion hypothesis yields consistent results with both the analysis on the 2020 data set only and the measurements of Paper I. HQ Tau B would be a 47.7 M_jup brown dwarf, orbiting with a period of 126 d, which is consistent at 1σ with the results of Sect. 3.3. HQ Tau C would exhibit a much larger mass (465 M_jup) and be on a 767 d orbital period (meaning a 2.16 AU semi-major axis), allowing it to reach the velocity determined in 2017, with the corresponding downward trend making its influence almost invisible during the six-month 2020 measurements. The systemic velocity resulting from this fit, of RV${}_{\mathrm{syst}}=12.1_{1.8}^{+1.4}$ km s⁻¹, is consistent at only 2σ with the median velocity of the Taurus cloud derived by Galli et al. (2019).

We performed the same test as in Sect. 3.3 to assess the robustness of the different solutions, that is, we compared the AICc of the regression and the MSE from a LOOCV test. For completeness, we also performed the multi-variate GPR assuming the stellar activity only (see Appendix C) and included it in the comparison. The results are summarised in Table 2, and are all in favour of the two-companion hypothesis.

In addition to a good fit of all our RV measurements, the orbital period of the tertiary is close to the occurrence of the long-lasting UX Orionis events reported by Rodriguez et al. (2017), which take place every 700–1000 days. One such event occurred immediately before the observations studied in Paper I, and we discussed them therein, putting forward possible interpretations, such as a change of the magnetospheric radius relative to the dust sublimation radius, producing dipper-like occultations, or as a sudden change of the disc’s vertical scale height, occulting the central star. This latter change in scale height could result from a disc wrapped by this component. As modelled by Cuello et al. (2019, 2020), stellar inclined retrograde flybys efficiently produce such a warp, in addition to a tilt of the disc, which remains after the warp dissipation, and could explain the discrepancy between the inclination of the stellar rotation axis (∼75°, Paper I) and the disc inclination (53.8°; Long et al. 2019).

To fully confirm the detection of such a tertiary companion, a set of observations similar to our 2020 dataset outside the inflexion point of the tertiary orbital motion is needed. Infrared high-resolution spectroscopy observations with a significantly higher S/N than the present dataset might also allow the direct detection of the tertiary’s spectra. Finally, given the mass(es) of the companion(s) and/or the semi-major axis of its (their) orbit(s), we assume that their influence on the accretion process of the primary studied in Paper I can be neglected. Even if a flyby by the hypothetical tertiary were able to induce mass-accretion-rate enhancements (Cuello et al. 2019), they only occur at low orbital inclination, contrary to the disc-warping phenomenon, and the variation is so drastic that the values derived in Paper I would have been clearly affected if HQ Tau had been observed during this phase.

5. Conclusions

In this paper, we present an investigation of the companion hypothesis to explain the low RVs measured for HQ Tau in our previous work (Pouilly et al. 2020, Paper I). We monitored the RV of HQ Tau over 6 months using high-resolution spectroscopy from four instruments in optical and IR frames (ESPaDOnS and SPIRou at CFHT, Neo-NARVAL at TBL, and SOPHIE at OHP). Finally, we modelled the stellar activity signal from the RV and BIS of the LSD profiles using multi-dimensional GPR to extract the Keplerian RV signal thanks to the pyaneti tool (Barragán et al. 2019, 2022).

This procedure yields the detection of a brown dwarf companion HQ Tau B ($M_{\mathrm{B}}={50}_{-32}^{+30}{M}_{\mathrm{jup}}$) orbiting the primary on a period of 126${}_{-55}^{+32}$ days without affecting its spectrum. Although this companion is consistent with the lunar occultation detection of Simon et al. (1987) and with the inner cavity in the disc suspected by Akeson et al. (2019), it does not allow us to recover the velocities measured in Paper I.

To address this issue, we investigated the possible presence of a third component and performed the same analysis including both datasets. This yielded a secondary component showing parameters consistent we the previous analysis and a tertiary of a mass of 465${}_{-106}^{+186}{M}_{\mathrm{jup}}$ and an orbital period of 767${}_{-75}^{+91}$ days. This orbital period is reminiscent of the quasi-period of the strong dimming events observed by Rodriguez et al. (2017), and this third component may therefore be linked to this behaviour, as well as to a long-term RV modulation, which would explain the low velocities measured in 2017.

To complete this study, we repeated the same analysis on both datasets assuming one companion only. This resulted in a component of ∼188 M_jup, with an orbital period of about 247 days. Although the fit itself is as good as the two-companion hypothesis, the period found is not consistent with the apparent modulation of about 100 d found in the 2020-2021 dataset using a periodogram analysis. Furthermore, the minimum RV modulation expected for the epoch of the 2017 observations based on the derived orbital parameters is still significantly higher than the values measured. The goodness of the fit is only reached thanks to a global decrease in the modulation induced by the stellar activity, which is unexpected. As both the ∼100-day modulation and the minimum at the epoch of the 2017 observations agree with the two-companion hypothesis, we favour this latter, although we cannot completely exclude the possibility that there is only one companion.

We therefore conclude that there is at least one companion to HQ Tau. We also strongly suspect that this companion is a brown dwarf accompanied by a third component, which would explain the low velocities reported in Paper I.

Acknowledgments

We thanks the anonymous referee for the very detailed and pertinent comments which drastically increased the robustness of this work. Based on observations obtained at the Canada–France– Hawaii Telescope (CFHT) which is operated from the summit of Maunakea by the National Research Council of Canada, the institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. The observations at the Canada–France–Hawaii Telescope were performed with care and respect from the summit of Maunakea which is a significant cultural and historic site. Based on observations made at Observatoire de Haute Provence (CNRS), France. We thank the TBL team for providing service observing with Neo-Narval. We thank A. Carmona for his help in reducing the SPIRou observations, especially for the tellurics correction. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 742095; SPIDI: Star-Planets-Inner Disk- Interactions). This work has made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna.

References

Appendix A: Additional tables

In this section we present the observation logs for the SOPHIE, Neo-Narval, SPIRou, and ESPaDOnS data sets in Table A.1, A.2, A.3, and A.4, respectively. The corresponding RV and BIS measurements are summarised in Table A.5.

Appendix B: Multi-variate GPR assuming the stellar activity only

Here we present the results of the multi-variate GPR analysis considering that the variations are induced by the stellar activity only. We are thus fitting only the GP, without the Keplerian modulation.

Fig. B.1. Fit of the RV (top) and BIS (bottom) curves with pyaneti assuming stellar activity only. The blue points are the measurements with their uncertainties, and the black curve is the GP fit of the stellar activity. The grey-shaded area is the 1σ uncertainty of the fit.

Appendix C: Multi-variate GPR on 2017 and 2020 radial velocity measurements

We present here the results of the analysis of the combined 2017 and 2020 datasets of HQ Tau discussed in Sect. 4. We first present the results obtained assuming the stellar activity only (Table C.1 and Fig. C.1). The results of the one-companion hypothesis are shown in Fig. C.2 and summarized in Table C.2. The similar figure and table for the two-companion hypothesis are presented in Fig. C.3 and Table C.3. The two hypothesis are discussed in Sect. 4.

Fig. C.1. Joint fit of the 2017 and 2020 RV curve with pyaneti assuming stellar activity only. The blue points are the measurements with their uncertainties, the black curve is the GP fit of the stellar activity. The grey-shaded area is the 1σ uncertainty of the fit.

Fig. C.2. Joint fit of the 2017 and 2020 RV curve with pyaneti assuming one companion. The blue points are the measurements with their uncertainties, the black curve is the GP fit of the stellar activity (purple) plus the Keplerian modulation by the companion (red). The grey-shaded area is the 1σ uncertainty of the fit. The x-axis has been cut for more visibility.

Fig. C.3. Joint fit of the 2017 and 2020 RV curve with pyaneti assuming two companions. The blue points are the measurements with their uncertainties, the black curve is the GP fit of the stellar activity (purple) plus the two Keplerian modulation by the two companions (red). The grey-shaded area is the 1σ uncertainty of the fit. The x-axis has been cut for more visibility.

Appendix D: Choice of the priors

In this section we present how we choose the priors indicated in Tables 1, B.1, C.1, C.2, and C.3. We are focusing this section on the Keplerian parameters for which we choose a normal distribution.

D.1. 2020 dataset

Orbital period P_orb

The Lomb-Scargle periodogram (Fig. 1) pointed out an orbital period at 111 days, given the low peak and high FAP, we used only 1 significant digit, meaning 100 days. Furthermore, the apparent modulation of the RV curve is consistent with such a period. Finally, from our first work on this object, we suspected an orbital period of the same order. Given these, the probability of having an orbital period of around 100 days seems higher than 50 or 150 days, explaining why we used a normal distribution instead of a linear one. We choose a σ of 50 to reach 0 days at 2-sigma and allow to exceed the 170 days of observation available.

Transit epoch T0

Assuming that the low modulation of the RV curve in Fig. 1 is an orbital modulation, we estimated the time at periastron around HJD 2 459 170, from which we removed the estimated 100-day orbital period to avoid negative cycles. Passing from the time at periastron to the transit epoch in our simple configuration means adding one-fourth of the period, meaning T0 = HJD 2 459 095. Here again, this value is motivated by previous measurements/estimates, explaining why we used a normal distribution. As T0 is modulo the orbital period, we choose the sigma to be lower than half the orbital period to avoid a misinterpretation of the posterior between T0 and T0+P_orb.

D.2. Joint 2017 and 2020 datasets

2-companion model

For this hypothesis, the B component is supposed to be the same as the one modulation the 2020 dataset only, we thus kept the same priors for this component. Concerning the orbital motion of the companion C, we did not have any indication, justifying the use of a uniform distribution for the orbital period. We used the approximate range of occurrence of the dimming events reported by Rodriguez et al. (2017) of about 750 days, which is not reported to be strictly periodic, that we broadened symmetrically to reach the 1000 days reported between two of these events. For the T0 priors, we first ran an MCMC with a uniform distribution, then we set the prior with the normal distribution around one of the peaks obtained to avoid multiple peaks in the posterior.

All Tables

All Figures

	Fig. 1. Radial velocity curve (top) and corresponding Lomb-Scargle periodogram analysis (bottom) of HQ Tau over the 20B semester. The colors of the radial velocity curve correpond to the different instruments used: NeoNarval at TBL (purple), SOPHIE at OHP (pink), ESPaDOnS (orange), and SPIRou (brown) at CFHT. On the bottom panel the Lomb-Scargle periodogram is shown in black and the FAP in grey.
In the text

	Fig. 2. Measurement of the BIS of the Stokes I profiles. The HJD is indicated at the top of each plot. The black lines are the Stokes I profiles and the grey dotted line delimits the top and bottom regions selected for the BIS computation. The blue points are the bisector in these regions and the red dots are the mean bisector of each region used to compute the BIS.
In the text

	Fig. 3. Lomb-Scargle periodogram of the BIS. The peak at f = 0.41 d−1 corresponds to a period of 2.42 d and has a FAP of 0.09.
In the text

	Fig. 4. BIS versus RV. The colour scale corresponds to the HJD, the black line is the linear fit, and the grey-shaded area is the 1σ uncertainty on the slope of the fit.
In the text

	Fig. 5. Posterior distribution of the orbital period (left) and of the semi-amplitude (right). The vertical red line shows the mean value and the grey-shaded area illustrates the 68.8% confidence level.
In the text

Fig. 6. Fit of the RV (top) and BIS (bottom) curves with pyaneti. The blue points are the measurements with their uncertainties. On the top panel the black curve is the GP fit of the stellar activity (purple) plus the Keplerian modulation by the companion (red). The grey-shaded area is the 1σ uncertainty of the fit. On the bottom panel the black curve is showing the GP fit of the stellar activity only with its 1σ uncertainty in grey.

In the text

	Fig. B.1. Fit of the RV (top) and BIS (bottom) curves with pyaneti assuming stellar activity only. The blue points are the measurements with their uncertainties, and the black curve is the GP fit of the stellar activity. The grey-shaded area is the 1σ uncertainty of the fit.
In the text

	Fig. C.1. Joint fit of the 2017 and 2020 RV curve with pyaneti assuming stellar activity only. The blue points are the measurements with their uncertainties, the black curve is the GP fit of the stellar activity. The grey-shaded area is the 1σ uncertainty of the fit.
In the text

	Fig. C.2. Joint fit of the 2017 and 2020 RV curve with pyaneti assuming one companion. The blue points are the measurements with their uncertainties, the black curve is the GP fit of the stellar activity (purple) plus the Keplerian modulation by the companion (red). The grey-shaded area is the 1σ uncertainty of the fit. The x-axis has been cut for more visibility.
In the text

	Fig. C.3. Joint fit of the 2017 and 2020 RV curve with pyaneti assuming two companions. The blue points are the measurements with their uncertainties, the black curve is the GP fit of the stellar activity (purple) plus the two Keplerian modulation by the two companions (red). The grey-shaded area is the 1σ uncertainty of the fit. The x-axis has been cut for more visibility.
In the text