© K. Pouilly et al. 2020

1. Introduction

Classical T Tauri stars (cTTs) are young stellar objects still surrounded by an accretion disk. They possess a strong magnetic field that truncates the inner disk at a distance of a few stellar radii above the stellar surface and drives accretion through funnel flows, a process referred to as magnetospheric accretion (see reviews in, e.g., Bouvier et al. 2007a; Hartmann et al. 2016). The kinetic energy of the infalling material is dissipated in a shock at the stellar surface, creating a localized hot spot. The star-disk interaction takes place within 0.1 au or less, a scale hardly resolved by current interferometers. An alternative approach to study this compact region is to monitor the variability of the system through contemporaneous photometric, spectroscopic, and spectropolarimetric observing campaigns focused on specific targets. Over the past decade, our group has reported several such studies (e.g., Bouvier et al. 2007b; Donati et al. 2007, 2019; Alencar et al. 2012, 2018).

This new study focuses on HQ Tau (RA = 04h35m, Dec. = +22°50′), a 1.9 M_⊙ T Tauri star located in the Taurus star forming region at a distance of 159 pc (Gaia Collaboration 2018). The system is moderately bright (V  ≃  12.5; Norton et al. 2007), has a K0-K2 spectral type (Nguyen et al. 2012; Herczeg & Hillenbrand 2014), and is so far considered as single (Pascucci et al. 2015, but see Sect. 4 below). From the SuperWASP campaign, Norton et al. (2007) measured a light curve modulation with a peak-to-peak amplitude of about 2 mag in the V-band and a period of 2.4546 days. Rodriguez et al. (2017) analyzed a 9 yr long KELT light curve of the system and reported several long-duration, non-periodic dimming events with an amplitude of about 1.5 mag and lasting for weeks to months, which led them to classify this object as UXor-like (Grinin 1992). Nguyen et al. (2012) derived a radial velocity of 16.65 ± 0.11 km s⁻¹, vsini of 48 ± 2 km s⁻¹, and a 10% width of the Hα line profile of 442 ± 93 km s⁻¹, indicative of ongoing accretion onto the star (see also Duchêne et al. 2017). Simon et al. (2016) measured EW(Hα) = 2.22 Å, from which they derived a mass accretion rate of Ṁ_acc = 2 × 10⁻⁹ M_⊙yr⁻¹. More recently, from high-resolution ALMA observations, Long et al. (2019) derived an inclination of 53.8 ± 3.2° for the circumstellar disk on a scale of ∼25 au. HQ Tau is the cTTs with the faintest disk of their Taurus sample in the mm range, and they report a depletion of dust toward the inner regions, perhaps indicative of an unresolved inner disk cavity. Akeson et al. (2019) similarly derived a low disk mass M_d = 0.4 × 10⁻³ M_⊙ and mm flux compared to single young stars of similar mass.

We resolved to launch a monitoring campaign on this object as a representative of intermediate-mass T Tauri stars (IMTTs), having a mass at the upper range of cTTs. Moreover, its K2 light curve exhibits a clear and smooth periodicity, which makes this target amenable to the variability approach developed here. Our main goal is to investigate whether the magnetospheric accretion process, which seems to be quite common among low-mass T Tauri stars applies as well to the intermediate-mass range. In Sect. 2, we describe the photometric and spectropolarimetric datasets used in this study. In Sect. 3, we present the results of the analysis, which includes deriving stellar parameters and investigating the photometric, spectral, and polarimetric variability. We discuss the results in Sect. 4, where we provide a framework for their interpretation. We also report long-term radial velocity variations, which indicate that the system is a spectroscopic binary. We conclude on the presence of the magnetospheric accretion process in this young system of the intermediate-mass group in Sect. 5.

2. Observations

We briefly describe here the origin of the photometric and spectropolarimetric observations that were obtained during the HQ Tau campaign. We also describe how they were processed.

2.1. Photometry

HQ Tau (EPIC 247583818) was observed by Kepler-K2 during Campaign 13, which took place over 80 days from March 8, 2017 to May 23, 2017. The observations were performed in a broad band filter (420–900 nm) with measurements taken at a cadence of 30 minutes. The K2 light curve was reduced by A.-M. Cody (Cody & Hillenbrand 2018) and we used the PDC version in this work.

Additional photometric observations of HQ Tau were secured in the Johnson’s VR_JI_J filters at the Crimean Astrophysical Observatory (CrAO) from December 22, 2016 to February 15, 2018, on the AZT-11 1.25m telescope equipped with the CCD camera ProLine PL23042. CCD images were bias subtracted and flat-field corrected following a standard procedure. We performed differential photometry between HQ Tau and a non-variable, nearby comparison star, 2MASS J04361769+2247125, whose brightness and colors in the Johnson system are 10.16 (V), 1.61 (V − R)_J, and 2.98 (V − I)_J. A nearby control star of similar brightness, 2MASS J04362344+2252171, was used to verify that the comparison star was not variable. It also provided an estimate of the photometric rms error in each filter, which amounts to 0.010, 0.012, and 0.011 in the VR_JI_J bands, respectively.

2.2. Spectroscopy and spectropolarimetry

High-resolution optical spectropolarimetry was obtained for HQ Tau using the Echelle SpectroPolarimetric Device for the Observation of Stars (ESPaDOnS; Donati 2003) at the Canada-France-Hawaii Telescope (CFHT) between October 28, 2017 and November 9, 2017. We obtained 14 spectra covering the 370–1000 nm range at a resolution of 68 000, reaching an S/N from 120 to 200 at 731 nm. Each observation consists in the combination of 4 individual spectra taken in different polarization states. All the data were reduced using the Libre-ESpRIT package (Donati et al. 1997), which provides an intensity Stokes I spectrum and a polarized Stokes V spectrum for each observation. The reduced spectra were normalized using a polynomial fit over selected continuum points (see Folsom et al. 2016), which produces a flat continuum, and simultaneously accounts for the spectral order overlap. The journal of observations is provided in Table 1.

3. Results

We describe in the following subsections the analysis of photometric, spectral, and spectropolarimetric variations observed in the HQ Tau system. We also derive the system’s fundamental properties.

3.1. Photometric variability

The 80 day-long detrended Kepler-K2 light curve of HQ Tau is shown in Fig. 1. It exhibits clear periodic modulation. A Lomb-Scargle periodogram analysis (Scargle 1982; Press & Rybicki 1989) yields a period P = 2.424 d, where the period estimate is derived from a Gaussian fit to the periodogram’s main peak, with an associated false alarm probability (FAP) of less than 10⁻⁴, computed following the Baluev (2008) approximation. Assuming that the formal error on the period is measured as the σ of the Gaussian fit (VanderPlas 2018), it would amount to 0.028 d. However, when the K2 light curve is folded in phase, it is clear that the phase coherency between the first and last minima is lost, meaning that their location differs by at least 0.05 in phase, as soon as the period is varied by more than 0.002 d. We therefore use this empirical estimate as the true error on the period. The light curve folded in phase with this period (see Fig. 1) displays a stable, sinusoidal-like pattern with a nearly constant amplitude, which we ascribe to the modulation of the stellar flux by a cool surface spot (Herbst et al. 1994). We therefore assume the photometric period P = 2.424 ± 0.002 d is the rotational period of the star. This period is consistent with the P1 = 2.42320 d period found by Rebull et al. (2020) whose work also highlighted a secondary period P2 = 5.02495 d. We do see a secondary peak in the periodogram of the detrended light curve at a period of 5.04 ± 0.13 d, with a FAP of less than 10⁻⁴ (see Fig. 1). We argue below that HQ Tau is probably a binary system, and we cannot exclude that the secondary period is to be associated with the rotation period of the companion.

Fig. 1. Top: HQ Tau’s detrended Kepler K2 light curve. Low frequency variations have been removed from the original light curve shown in Fig. 2. The spike appearing at J.D. 2 457 890.0 is presumably an instrumental defect. Middle: HQ Tau’s light curve folded in phase with P = 2.424 d. The origin of phase is taken as JD 2 457 823.81 to correspond to the epoch of maximum brightness. Bottom: lomb-Scargle periodogram obtained from HQ Tau’s Kepler K2 light curve.

Figure 2 shows the full HQ Tau light curve over 2 epochs, from December 2016 to March 2018, including the K2, before detrending, and CrAO datasets, as well as publicly available photometry from the ASAS-SN and AAVSO surveys. The K2 light curve was rescaled to the 7 CrAO measurements taken contemporaneously. The normalized K2 fluxes were thus converted to magnitudes applying a 12.3 mag zero-point offset. Over the time frame where they overlap, the amplitude of the K2 light curve appears somewhat smaller than that of the V-band CrAO light curve, presumably reflecting the longer effective wavelength of the K2 wide-band filter. Interestingly, the low-frequency part of the K2 light curve, most notably the brightness “bump” seen over its first part, is recovered over the slightly overlapping parts of the CrAO and ASAS-SN light curves, which indicates intrinsic longer term variations superimposed onto the spot modulation.

Fig. 2. HQ Tau’s V-band light curve. CrAO measurements are shown as filled grey circles, the K2 light curve as a continuous line, AAVSO measurements as filled triangles, and the ASAS-SN rescaled g- and V-band datasets as crosses and open circles, respectively. The upper panel shows the light curve over 2 observing seasons while the lower one displays the photometric variability of the source around the CFHT/ESPaDOnS observations.

Photometric data points were also available from the ASAS-SN survey (Jayasinghe et al. 2019). HQ Tau was observed in both the V-band and g-band filters during this period. Taking advantage of the overlap between the V- and g-band light curves, the g-band measurements were arbitrarily shifted by −0.8 mag in order to match the V-band ones. The resulting ASAS-SN light curve contains a few measurements during the first half of the ESPaDOnS observations and many additional measurements before and after the run was completed. It thus provides an estimate of the amount of variability the source exhibited at the time of the spectroscopic observations, performed about 6 months after the K2 monitoring. Overall, it suggests relatively mild variability at this epoch, amounting to about 0.2 mag in amplitude in the V-band. A significant dimming event, which lasted for about a month just prior to the start of the spectroscopic observations had apparently ceased by the time of the ESPaDOnS measurements. We notice, however, that the mean flux level was changing over the spectroscopic run, with the source being about 0.3 mag brighter toward the end of the run. Both the amplitude of variability and the low-frequency variations are not unlike those seen in the spot-driven K2 light curve obtained several months earlier. This suggests that at the time of the ESPaDOnS observations, the source was in a state of relatively mild periodic variability, and probably not in a drastic dipper state as reported at some other epochs by Rodriguez et al. (2017), in spite of the occurence of a recent dimming event.

Indeed, we verified that the ASAS-SN g-band and V-band datasets restricted to the 7 days following the ESPaDOnS run (JD 2 458 067.8–2 458 073.7) exhibit a smoothly varying light curve consistent with a period of 2.424 days when folded in phase. Using the same epoch for the origin of phase as for the K2 light curve above (i.e., JD 2 457 823.81), we find that the photometric minimum of the ASAS-SN dataset occurs around phase ∼0.56, which we estimated by interpolating the two lowest photometric measurements occuring at phase 0.47 and 0.66, respectively. Hence, the ASAS-SN photometric minimum is not far from phase 0.50 of the K2 photometric minimum, the slight difference being easily accounted for by the 0.002 d uncertainty on the K2 period. The near conservation of phase over the 155 day-long temporal gap stretching between the end of the K2 observations and the beginning of the ESPaDOnS run suggests the modulation results from a long lived, relatively stable spot structure. For the rest of the paper, we thus use the following ephemeris:

$$\begin{array}{c}\hfill \mathrm{HJD}(d)=2457823.81+2.424E,\end{array}$$(1)

where E is the rotational phase of the system.

Color variations associated to the brightness changes are shown in Fig. 3. In both the (V − R)j and (V − I)j colors, the source becomes redder when fainter. The color slope of the small scale variations likely results from spot modulation (Herbst et al. 1994; Venuti et al. 2015). A single deep faintening event was recorded on JD 2 458 025.6 in the VR_JI_J filters of the CrAO dataset, and is confirmed by the single-filter ASAS-SN light curve (see Fig. 2). The color plot shows that as the brightness decreased by about 1.2 mag in the V-band, the system became much redder, with a color slope close to that expected for extinction by ISM-like grains. This suggests that the dimming event was caused by circumstellar dust crossing the line of sight.

Fig. 3. Color-magnitude plots from CrAO VRJIJ photometry. Left: (V − R)J vs. V. Right: (V − I)J vs. V. The source gets redder when fainter. The dashed line is the expected reddening line for ISM extinction. The photometric measurements are color-coded by Julian date to distinguish the 2 observing seasons: lighter points correspond to the first epoch, darker points to the second one (see Fig. 2). We note that the deep faintening event of about 1.2 mag in the V-band is accompanied by a strong reddening of the system.

3.2. Stellar parameters

We used the ESPaDOnS high resolution spectra to derive HQ Tau’s stellar parameters, namely effective temperature (T_eff), and the radial (V_r), rotational (vsini), and microturbulent (v_mic) velocities. We averaged the 14 ESPaDOnS spectra gathered during the campaign and fit synthetic spectra calculated with the ZEEMAN code (Landstreet 1988; Wade et al. 2001; Folsom et al. 2012) based on MARCS stellar atmosphere grids (Gustafsson et al. 2008), VALD line lists (Ryabchikova et al. 2015), and including the same oscillator strength corrections as those used in Folsom et al. (2016). We explored a range of T_eff, vsini, V_r, and v_mic, and obtained the best fit to HQ Tau’s mean spectrum through a χ² minimization procedure using a Levenberg-Marquardt algorithm (see Folsom et al. 2013). The fit was performed on 11 independent spectral windows spanning the range from 422 to 754 nm, each with a width ranging from 4 to 10 nm. Each window contains well resolved, relatively unblended photospheric lines, and is devoid of emission features. The spectral windows are shown in Fig. 4.

Fig. 4. Upper six rows: spectral windows selected from HQ Tau’s mean ESPaDOnS spectrum to derive stellar parameters. Bottom row: as an example, HQ Tau’s mean spectrum (blue) fitted with a ZEEMAN synthetic spectrum (orange) in the 547–557 nm window.

In order to derive the stellar properties from spectral fitting, we first fixed the macroturbulent velocity to 2 km s⁻¹, log g to 4.0, and assumed solar metallicity, in other words, values that are typical of low-mass PMS stars in the solar neighborhood (Padgett 1996; James et al. 2006; Santos et al. 2008; Taguchi et al. 2009; D’Orazi et al. 2011). An example of the resulting fit is shown in Fig. 4. Once stellar parameters were obtained for each spectral window, we averaged them and computed the rms dispersion, and then removed the windows that yielded results beyond 1σ from the mean value. Usually, between 3 or 4 windows were thus rejected. We thus derive T_eff = 4997 ± 160 K, vsini = 53.9 ± 0.9 km s⁻¹, V_r = 6.64 ± 0.71 km s⁻¹, and v_mic = 1.4 ± 0.2 km s⁻¹. The T_eff of 4997 K we derive is consistent with a spectral type K0-K1 according to Pecaut & Mamajek (2013) conversion scale.

We next obtained an estimate of HQ Tau’s luminosity, based on its 2MASS J-band magnitude (J = 8.665 ± 0.024, Skrutskie et al. 2006), using the Gaia parallax (π = 6.304 ± 0.213 mas, Gaia Collaboration 2018). Extinction was computed from the (V − J) color excess, using intrinsic colors from Pecaut & Mamajek (2013), which yields A_V = 2.60 ± 0.1 mag and A_J = 0.73 ± 0.03 mag, consistent with the value found by Herczeg & Hillenbrand (2014). We applied bolometric corrections from Pecaut & Mamajek (2013), and obtained L_⋆ = 3.91 ± 0.65 L_⊙, and R_⋆ = 2.64  ± 0.31 R_⊙.

We then plotted HQ Tau in an Hertzsprung-Russell (HR) diagram and used a grid of CESTAM evolutionary models (Marques et al. 2013; Villebrun et al. 2019) to derive its mass and internal structure (see Fig. 5). From bilinear interpolation in model grids, we obtain $M_{\star }=1.{87}_{-0.55}^{+0.21}{M}_{\odot }$, and an age of ∼2 Myr. According to evolutionary models, the star, is partly radiative with $M_{\mathrm{rad}}=0.{55}_{-0.55}^{+0.35}{M}_{\star }$ and $R_{\mathrm{rad}}=0.{51}_{-0.51}^{+0.14}{R}_{\star }$, where M_rad and R_rad are the mass and radius of the radiative core, respectively.

Fig. 5. HQ Tau’s position in the HR diagram. The red cross illustrates HQ Tau’s position with corresponding uncertainties. The black curves are the evolutionary tracks from CESTAM PMS models, with the corresponding mass indicated on the left. The blue curves show the position where the radiative core of the star reaches 75%, 50% and 25% of the stellar mass.

We find that HQ Tau is a fast rotator with vsini = 53.9 km s⁻¹. Combining the stellar radius with the star’s vsini and the K2 rotational period, we derive the inclination angle of the rotational axis onto the line of sight, $i={75}_{-17}^{+15}$ deg, indicative of a highly inclined system.

3.3. Spectroscopic properties and variability

We investigate the spectral variability of the system by analyzing the 14 ESPaDOnS spectra taken over 12 nights. Table 2 summarizes radial velocity and line equivalent width measurements.

3.3.1. Radial velocity variations

Radial velocity measurements were obtained by cross-correlating each HQ Tau’s spectrum with the spectrum of a spectral template. We used the ESPaDOnS observations of the Hyades cluster member Melotte 25-151 as a template, a slowly rotating (vsini = 4.8 km s⁻¹) K2-type star with V_r = 37.98 km s⁻¹ and T_eff = 4920 K (Folsom et al. 2018), reduced and normalized in the same way as HQ Tau. We computed the cross correlation function (CCF) over 5 spectral windows (542–547, 558–563, 585–587, 608–613, and 639–649 nm), and fit a Gaussian profile to derive the radial velocity difference between HQ Tau and the template. We averaged the results over the 5 spectral windows to get a mean value of V_r and its rms uncertainty. The results are listed in Table 2 and the radial velocity curve is shown in Fig. 6.

Fig. 6. From top to bottom: radial velocity, Hα line flux, and average longitudinal magnetic field as a function of Julian date (left) and rotational phase (right). The sinusoidal fit, forced at the K2 photometric period, to radial velocity variations is shown as a black dash-dotted curve in the top panels. The Bl = 0 G line is drawn for reference in the bottom panels.

The radial velocity appears to be modulated and a sinusoidal fit yields a period of 2.48 ± 0.16 d using the 1σ confidence level on χ² minimization, consistent with the stellar rotation period within uncertainties. We therefore ascribe this modulation to surface spots. The phased radial velocity curve is shown in Fig. 6. Indeed, the sinusoidal fit indicates that the mean radial velocity (⟨V_r⟩=7.22 ± 0.27 km s⁻¹) occurs around phase 0.6, which is expected when the spot faces the observer (Vogt & Penrod 1983) and is consistent with the photometric minimum of the contemporaneous ASAS-SN light curve. The amplitude of the V_r modulation amounts to hardly a tenth of the star’s vsini, which suggests it is mostly driven by a large area and high latitude cold spot. We also notice a regular downward drift of V_r with an amplitude of 2–3 km s⁻¹ over a timescale of 10 days. We will come back to this feature in Sect. 4.

3.3.2. Emission line profiles

We computed residual emission line profiles by subtracting the rotationally broadened photospheric profiles of the template Melotte 25-151 from HQ Tau’s profiles. Hα, Hβ, and the CaII infrared triplet (IRT) are the only lines in HQ Tau’s spectrum exhibiting significant emission flux. The 3 lines of the CaII IRT are similar in shape but the line at 866.2 nm may be affected by a hydrogen line of the Paschen series and the line at 849.8 nm exhibited a variation far lower than the other two. We therefore focused on the 854.2 nm component. Figure 7 shows the Hα, Hβ, and the selected CaII IRT residual line profiles. The raw line profiles are shown in Fig. A.1.

Fig. 7. Top: residual Hα (left), Hβ (middle), and CaII 854.2 nm (right) emission lines of HQ Tau spectra. Bottom: corresponding 2D periodograms with the mean profile superimposed as a white curve. The maximum power of the periodogram is in yellow and the minimum in dark blue. The red and blue wings of Hα and the red wing of Hβ and CaII, as well as the center of the 3 line profiles, exhibit a periodic modulation at f ≈ 0.4 d−1, corresponding to the rotational period of the star. The power seen at f ≈ 0.6 d−1 is the 1-day alias induced by the observational window.

The Hα line profile exhibits high velocity wings to the blue and red sides, up to about 300 km s⁻¹. The red wing displays large variability, with at times signs of high velocity redshifted absorptions reaching below the continuum, which are Inverse P Cygni (IPC) profiles. The blue wing also displays significant variability. The Hβ line is dominated by deep IPC signatures extending up to about 300 km s⁻¹ seen in nearly half of the observations. The CaII line also exhibits strong variability from the line center to its red wing, with however only marginal IPC components below the continuum level.

We computed periodograms in each velocity bin across the line profiles, the result of that is a 2D periodogram for each line, which is shown in Fig. 7. In all 3 lines, a clear periodicity is seen at a frequency of 0.4 d⁻¹ (P = 2.5 d) with a FAP reaching ∼10⁻⁴ and a typical value of ∼10⁻², corresponding to the stellar rotational period, extending from the line’s central peak all the way to the red wing. This is a clear indication that the IPC components are modulated by stellar rotation, as expected for funnel flow accretion. We do notice significant power as well at about the same frequency at highly blueshifted velocities (∼ − 200 km s⁻¹) in the Hα profile, which might also be present in the CaII line profile. The peak at a frequency of 0.6 d⁻¹ is the 1-day alias of the 0.4 d⁻¹ frequency. It appears clearly here, though at a higher FAP of 0.1, due to the night-to-night sampling of ESPaDOnS spectra.

Figure 8 shows the residual Hα and Hβ profiles ordered by day and by phase. The phase ordering illustrates well the periodicity of IPC components: they appear over the 5 rotational cycles covered by the spectral series at specific phases, from 0.69 to 0.02 in both profiles, with a maximum depth around phase 0.90. The depth of the IPC components appears to slightly vary from one rotational cycle to the next, being stronger at phase 4.90 than at phase 2.85 for instance.

Fig. 8. Residual Hα (left panels) and Hβ (right panels) profiles sorted by rotational phase and by day of observation. The successive rotational cycles are shown by different colors, as in Fig. 6.

We investigate the relationship between the various components seen in the line profiles by computing correlation matrices (Johns & Basri 1995a; Oliveira et al. 2000; Alencar & Batalha 2002; Kurosawa et al. 2005). Correlation matrices consist of the computation of Pearson’s linear correlation coefficient on line intensity between 2 velocity channels of the same or different line profiles. The coefficient approaches 1 for a strong correlation, 0 for no correlation or −1 for anticorrelated intensity variations. Correlation matrices are shown in Fig. 9.

Fig. 9. Correlation matrices of the line profiles. (a) Hα vs. Hα, (b) Hβ vs. Hβ, (c) Hβ vs. Hα, (d) CaII vs. CaII, (e) Hα vs. CaII, (f) NaD vs. NaD.

A strong correlation is seen between velocity channels within the red wing of all profiles, and within the blue wing for Hα and CaII. These correlations appear as bright squares in the correlation matrices, located on each side of the line center and extending to high velocities. As discussed above, the strong correlation within the red side of the line profiles is probably related to the periodic appearance of IPC components, which extends from close to the line center up to +300 km s⁻¹. Strikingly, little correlation is seen between the red and blue wings of the line profiles, suggesting that their variability is driven by unrelated physical processes, presumably accretion for redshifted velocity channels and outflows for blueshifted ones. Yet, an interesting anti-correlation appears between a restricted range of blueshifted velocity channels, from about −130 to −220 km s⁻¹, and the redshifted part of the line profile. This is clearly seen in the Hα correlation matrix and also appears in the CaII matrix. This indicates that, as the IPC components appear on the red side of the line profile, the high velocity part of the blue wing becomes more intense, perhaps the signature of an accretion-driven high velocity outflow (e.g., Johns & Basri 1995b).

3.3.3. The NaI D profile

The NaI doublet is sensitive to wind signatures (Mundt 1984), and can thus help to better understand the accretion-ejection connection. The NaI D lines are seen in absorption in HQ Tau’s spectrum and exhibit significant variability (see Fig. 10). As the two lines of the doublet have similar profiles and behavior, we focus here on the 589.0 nm line (NaI D2), leaving aside the 589.6 nm twin line, which lies at the edge of a spectral order.

Fig. 10. Top: NaI D2 line appearing in the 14 spectra of HQ Tau. The black dash-dotted line is a synthetic NaI D2 line produced using the ZEEMAN code with HQ Tau’s parameters. Bottom: 2D periodogram of the NaI D2 line.

Two prominent, relatively narrow, and apparently stable absorption components are seen around the line center. The stable component located at +10 km s⁻¹ (V_helio ≃ 17 km s⁻¹) was previously reported by Pascucci et al. (2015) and ascribed to local interstellar absorption. This narrow and stable absorption component is also present at the same redshifted velocity in all our spectra in the KI 770 nm line profile. We notice however significant variability on both the blue and red sides of this component in the NaI D2 line, over a range of about ±10 km s⁻¹, which indicates an additional source of variable absorption, linked to the stellar system itself. These transient components are also seen in the KI profile. The second narrow absorption component is blueshifted and located at −20 km s⁻¹ (V_helio ≃ −13 km s⁻¹). This component is also apparent in the NaI D line profile reported in Pascucci et al. (2015), but is not discussed there. CI Tau, which is located only 0.5 degrees away from HQ Tau, does not exhibit this second component, while it does display the redshifted one. This suggests the blueshifted narrow absorption component is not due to interstellar cloud absorption. Yet, its narrow width compared to photospheric lines and its stability suggest it is not related to the inner variable system. This component could conceivably be the signature of the distant, cold part of a stellar wind, as discussed in Mundt (1984).

The variability seen across the NaD line profile is relatively complex, with different components. The 2D periodogram shown in Fig. 10 reveals that the red side of the profile, from about +40 to +150 km s⁻¹, varies periodically, at the stellar rotation period. The comparison of successive NaI D2 profiles ordered in phase, shown in Fig. 11, suggests this is due to the appearance of IPC components seen between phase 0.69 and 0.02, as reported above for Balmer and CaII emission lines. In contrast, the variability seen in the blue wing of the line profile is not periodic. Sporadic blueshifted absorption components appear in nearly half of the observations, over velocity channels ranging from about −120 to −30 km s⁻¹. Figure 11 shows the occurrence of deep, high velocity blueshifted absorption components from JD 2 458 060 to 2 458 066 covering velocity channels from about −120 to −50 km s⁻¹, with evidence for a gradual shift in velocity over this timescale. Lower velocity absorption components are also seen for a couple of days from JD 2 458 056 to 2 458 057 around −60 to −30 km s⁻¹. None of these blueshifted components appear to last over the whole observing period. They appear to be transient phenomena occurring only during specific cycles, are of limited duration, and are not modulated by stellar rotation. These are presumably the signatures of episodic outflows.

Fig. 11. NaI D2 line profiles sorted by rotational phase (left) and by day of observation (right). The successive rotational cycles are shown by different colors, as in Fig. 6.

Finally, Fig. 9f presents the correlation matrix of the Na D2 line. It resembles the matrices of the other emission lines in the red part of the profile, being similarly modulated by the IPC components. It does not show any correlated variability in the blue wing, as none is expected from the sporadic, independent outflows components analyzed above.

3.3.4. Mass accretion rate

We computed the mass accretion rate using the relationship between the accretion luminosity and the luminosity in the Hα residual line (Alcalá et al. 2017). Duchêne et al. (2017) notice that the Hα emission line of HQ Tau might not be a good accretion indicator, due to its weakness, though its width and double peaked shape favor an accretion origin. In order to confirm the mass accretion rate values found with Hα, we also computed the mass accretion rate using the Hβ, and CaII IRT residual lines. The line flux was computed from the line equivalent width as F_line = F₀ ⋅ EW × 10^−0.4m_λ, where F_line is the line flux, F₀ is the reference flux in the selected filter, EW the line equivalent width, and m_λ is the extinction-corrected magnitude of HQ Tau in the selected filter, namely R_C for Hα, R_J for the Ca II IRT, and B for Hβ. The line luminosity is then derived from L_line = 4πd²F_line, where d is the distance to HQ Tau. The accretion luminosity follows from the line luminosity by log(L_acc) = alog(L_line)+b, where a and b are coefficients taken from Alcalá et al. (2017). The mass accretion rate is obtained from:

$$\begin{array}{c}\hfill L_{\mathrm{acc}}=\frac{G{M}_{\star }{\dot{M}}_{\mathrm{acc}}}{R_{\star }}[1-\frac{R_{\star }}{R_t}],\end{array}$$(2)

where R_t is the magnetospheric truncation radius, typically 5R_⋆ for cTTs (Bouvier et al. 2007a).

As the equivalent width of Balmer lines is strongly affected by the IPC components, we computed the mean value of L_acc and Ṁ_acc by averaging the results excluding the observations showing strong IPC profiles. The mass accretion rates we deduce from the various lines are consistent within 2σ. This yields L_acc = 0.019 ± 0.005 L_⊙ and Ṁ_acc = 1.26 ± 0.35 × 10⁻⁹ M_⊙yr⁻¹ for the system. As EW(Hα) is relatively weak in HQ Tau, the emission line flux might include a non negligible chromospheric contribution. In order to quantify it, we considered a sample of wTTSs with similar spectral type and vsini and (or) rotation period as HQ Tau, which were chosen from the list of James et al. (2006). We thus derived ⟨EW(Hα)⟩ = 0.7 ± 0.1 Å for these wTTSs, assumed this value to be the chromospheric contribution, and subtracted it from HQ Tau’s EW(Hα). We then recomputed Ṁ_acc from the different line estimates to obtain Ṁ_acc = 1.16 ± 0.35 × 10⁻⁹ M_⊙yr⁻¹. The temporal variation of the Hα line flux, a proxy for Ṁ_acc, is shown in Fig. 6. It is clearly modulated and a sinusoidal fit yields a period of 2.4 d, consistent with the stellar rotational period. The peak-to-peak amplitude of the modulation is about a factor 3, typical of cTTs Ṁ_acc variability on these timescales (Venuti et al. 2014).

Another signature of accretion is the so-called veiling. Veiling is an additional continuum component emitted by the accretion shock, which fills the photospheric lines. In order to compute the amount of veiling in HQ Tau’s spectrum, we used the non-accreting star Melotte 25-151 as a template. Varying the amount of continuum line filling over the spectral window 639–649 nm in the template spectrum, we matched HQ Tau’s spectrum by a χ² minimizing method:

$$\begin{array}{c}\hfill I(\lambda )=\frac{I_t(\lambda )(1+r)}{1+I_t(\lambda )r},\end{array}$$(3)

where I_t is the intensity spectrum of Melotte 25-151, the template, and r the fractional veiling that is the excess continuum flux divided by the stellar continuum flux at a given wavelength. We derived a weak veiling at 640 nm ranging from 0.15 to 0.20 in the 14 HQ Tau spectra, as expected for the modest mass accretion rate and bright photosphere of the system, with uncertainties of order of 0.20, thus preventing us from detecting any significant temporal variations.

3.4. Spectropolarimetry

We analyzed the ESPaDOnS spectropolarimetric data to study the magnetic field properties of HQ Tau. We used the Least Square Deconvolution method (LSD) from Donati et al. (1997) to compute the mean Stokes I and V photospheric line profiles. This method increases the S/N by extracting the Zeeman signature in many photospheric lines and averaging them. The parameters used for the LSD computation are the mean wavelength, intrinsic line depth, and Landé factor, set at 640 nm, 0.2, and 1.2, respectively, as in Donati et al. (2010). We extracted a list of spectral lines from VALD atomic database for a star with HQ Tau’s fundamental parameters, and selected absorption lines located between 450 and 850 nm. We then removed lines contaminated with emission, blended with strong broad lines, or affected by telluric lines. The LSD Stokes I and V profiles are shown in Fig. 12. The typical rms is 3.5 × 10⁻⁴ and 1.4 × 10⁻⁴ for the Stokes I and V profiles, respectively, computed from more than 9000 photospheric lines.

Fig. 12. HQ Tau’s Stokes I(left), V(middle), and Null (right) profiles (black) ordered by increasing fractional phase, and the fit produced by the ZDI analysis (red). The fractional number next to each profile indicates the rotational phase.

The shape of the Stokes I profiles displays a clear modulation over the rotational cycle. A periodogram analysis shows the same peak at a frequency 0.4 d⁻¹ as for emission lines, ranging from −25 km s⁻¹ to 60 km s⁻¹, with a FAP reaching 10⁻⁴. This indicates that the modulation is periodic at the stellar rotational period, as expected from a large rotating dark spot at the stellar surface. The profile is nearly symmetric at phase 3.69, confirming the spot is facing the observer around phase 0.65, as defined above from the V_r curve (see Sect. 3.3.1). The Stokes V signal is much weaker than the Stokes I profile. Nevertheless, a clear magnetic signature is seen at several phases, with a complex shape.

The average surface longitudinal magnetic field, B_l, was computed from the Stokes I and V profiles, as:

$$\begin{array}{c}\hfill B_l=-2.14\times {10}^{11}\times \frac{\int vV(v)\mathrm{d}v}{\lambda gc\int (1-I(v))\mathrm{d}v},\end{array}$$(4)

where B_l is given in Gauss, v is velocity relative to line center, λ in nm is the mean wavelength, and g the Landé factor chosen for the LSD computation (Donati et al. 1997; Wade et al. 2000). We computed the error by propagating the uncertainties in the trapezoidal integration over the range ±60 km s⁻¹.

Figure 6 shows the variation of the longitudinal magnetic field as a function of time and rotational phase. B_l varies from −70 G to +68 G. We notice that the minimum of B_l is reached at ϕ ∼ 0.7, consistent with the rotational phase at which the dark spot that modulates the radial velocity faces the observer, and reaches a maximum around phase 0.9 when the IPC components are best seen in the emission line profiles. Apart from these extrema, there is little evidence for B_l to be modulated, possibly the result of a complex magnetic topology at the stellar surface and (or) the strongest field regions being hidden in dark spots. We ran a similar analysis in the CaII IRT 854.2 nm line profile. However, we did not detect any significant Stokes V signature within the line profile and our data can only provide a 3σ upper limit of 450 G to the longitudinal magnetic field component in the line.

The full Zeeman Doppler Imaging (ZDI, Donati et al. 2011, 2012) analysis was performed based on the LSD Stokes I and V profiles, using the code described in Folsom et al. (2018). The resulting best fit of the Stokes I and V profiles is shown in Fig. 12.

This first step provides a Doppler image shown in Fig. 13. The fitting process starts from a uniform brightness distribution across the stellar disk and iteratively adds darker and brighter features onto the stellar surface to reproduce the shape of Stokes I profile. For each cell on the stellar surface, the local profile is assumed to be a Voigt profile and the brightness of the pixel is adjusted by minimizing χ² and maximizing entropy. The solutions for brightness and magnetic field distributions are not unique based on χ². Therefore, an additional constraint is added using entropy. ZDI uses the algorithm of Skilling & Bryan (1984) for maximizing entropy while minimizing χ². We adopted 0.66, 2.8 km s⁻¹, and 1.9 km s⁻¹ for the limb darkening coefficient, the Gaussian and Lorentzian widths of the Voigt profile, respectively. The latter values were chosen based on theoretical values found in Gray (2005) for a K0 spectral type star around 600 nm and adapted to fit the observed profiles. The Doppler image reveals a large cold spot extending over phases 0.4 to 0.8 in longitude, ranging from 90° to about 50° in latitude, and covering about 20% of the stellar visible hemisphere.

Fig. 13. Upper panel: surface brightness map of HQ Tau. The color scale goes from dark brown to light yellow, for dark to bright features. Lower panels, from top to bottom: radial, azimuthal, and meridional surface magnetic field distribution. Ticks on top of the panels indicate the phase of observations.

ZDI then takes the brightness map as input and fits the Stokes V profiles by adjusting the spherical harmonic coefficients that describe the magnetic field (see Donati et al. 2006). It is also possible to fit both intensity and polarisation profiles simultaneously but the resulting maps are essentially the same. Taking advantage of the large vsini of the star, which provides spatial resolution on the stellar surface, the spherical harmonic expansion was carried out to the 15th order in l (Folsom et al. 2016). Letting the stellar parameters freely vary in the fitting process, we obtain an inclination of 50 ± 5 °, a period of 2.453${}_{-0.008}^{+0.006}$ d, and a vsini of 51.2 ± 0.5 km s⁻¹. All uncertainties were obtained using the 68.27 % confidence level on χ². All these values are within 3σ of those we derived in the previous sections, although the inclination we derive from ZDI is on the lower side of the estimate we obtained in Sect. 3.2 by combining the rotational velocity and period with the stellar radius. This suggests that the stellar radius may be underestimated by about 20%. The solution that minimizes χ² did not allowed us to constrain the differential rotation dΩ. We thus fixed dΩ = 0.0 rad d⁻¹ to derive the other parameters, as the derived period was consistent within 3σ with the photometric one.

The large-scale magnetic reconstruction we recover from the analysis of LSD profiles reveals a mostly poloidal field, with the toroidal component contributing 25% to the total magnetic energy. The main structure of the radial magnetic field extends from phase 0.3 to 0.8, from 90° to 60° in latitude, and reaches a strength of −562 G locally. The dipolar field amounts to only 14% of the poloidal component, is tilted by 31.5° from the rotation axis, and reaches a maximum strength of 63 G at the stellar photosphere. The quadrupole and octupole components contribute respectively 9.7% and 8.3% to the total magnetic flux. The complex topology and moderate strength of HQ Tau’s magnetic field is reminiscent of those reported by Hussain et al. (2009) and Villebrun et al. (2019) for IMTTs.

We caution that the values reported here from the LSD analysis are lower limits to the strength of the magnetic field that interacts with the disk. Previous similar studies have shown that the ZDI analysis of emission lines arising from the accretion shock located at the foot of funnel flows, most notably HeI 588 nm, yields a much stronger magnetic field strength than that deduced for the large-scale field at the photospheric level from LSD analysis (e.g., Donati et al. 2019). Unfortunately, HQ Tau’s spectrum does not show the post-shock HeI line in emission, owing to the modest accretion rate producing a relatively weak shock seen against a bright photosphere.

4. Discussion

We selected HQ Tau among many Taurus sources observed by K2 (Rebull et al. 2020) for a spectropolarimetric follow-up campaign in order to investigate the magnetospheric accretion process in a representative member of the class of relatively massive T Tauri stars, the so-called IMTTs (Calvet et al. 2004). With a mass of M_⋆ = 1.9 M_⊙, a short rotational period of 2.424 d and vsini = 53.9 km s⁻¹, HQ Tau’s properties are intermediate between the cool low-mass T Tauri stars and the hotter intermediate-mass Herbig Ae stars. Indeed, PMS models suggest HQ Tau’s interior is already partly radiative, with R_rad/R_⋆ ∼ 0.51, as the star is transiting from the Hayashi to the Henyey track in the HR diagram. The mass accretion rate we derive, Ṁ_acc ∼ 1.16 × 10⁻⁹ M_⊙ yr⁻¹, is however relatively low and more typical of low-mass TTS (e.g., Mendigutía et al. 2011). Only few such IMTTs had been monitored for their magnetic properties so far (e.g., Hussain et al. 2009) and the goal here was to explore the possible extension of the magnetospheric accretion process that is ubiquitous among low-mass T Tauri stars to the higher mass range.

HQ Tau is undoubtedly a member of the Taurus star forming region (Luhman 2018). Yet, the mean radial velocity we measured, ⟨V_r⟩=7.22 ± 0.27 km s⁻¹, is significantly different from the radial velocity distribution of Taurus members, ⟨V_r⟩=16.3 ± 6.43 km s⁻¹ (Bertout & Genova 2006). This led us to suspect that HQ Tau might not be a single source, as usually assumed, and the small but regular drift in V_r seen during the ESPaDOnS run may be additional evidence for that (see Fig. 6). Historically, Simon et al. (1987) reported the system to be a tight binary from lunar occultation, with a separation of 4.9 ± 0.4 mas, later revised to 9.0 ± 2 mas by Chen et al. (1990). However, the former authors state that “the binary nature of HQ Tau is not obvious” and indeed this result was not confirmed by following studies (Richichi et al. 1994; Simon et al. 1995, 1996; Mason 1996).

We searched for additional radial velocity measurements in the literature. They are summarized in Table 3, together with new measurements, kindly obtained for us recently by L. Hillenbrand and H. Bouy. The results are illustrated in Fig. 14. HQ Tau exhibits clear velocity variations. Measurements obtained from 2006 to 2017 show V_r oscillations between about 7 and 22 km s⁻¹ on a timescale of years. However, the two most recent measurements taken 50 days apart in late 2019, show the same amplitude. This suggests that HQ Tau is a short period type 1 spectral binary (SB1), whose orbit determination awaits additional measurements. We note that the companion must be relatively faint compared to the primary in the optical, as we see no evidence for a double-lined system in the Stokes I LSD profiles. Nevertheless, the detection of a secondary period at 5.03 d in the K2 light curve (see Sect. 3.1) suggests that the companion may account for a low-level contribution to the optical flux of the system. We caution that we have implicitly assumed in the above analysis that the companion’s contribution to the emission line flux and to the ZDI reconstruction process could be neglected. Although this assumption may not be fully verified, we have currently no way to estimate the flux contribution of the companion to the system. Finally, we suggest that the faint mm flux of HQ Tau’s disk and its suspected inner cavity (Long et al. 2019; Akeson et al. 2019) could conceivably be related to a low-mass companion orbiting within the inner disk and, at least partly, clearing it.

Fig. 14. Radial velocity of HQ Tau as a function of date of observation. The error bars are indicated, and can be smaller than the symbol size.

On a much shorter timescale of days, the radial velocity of the system is modulated by surface spots with a period consistent with the photometric period of 2.424 d derived from the K2 light curve. From the additional ASAS-SN photometry obtained at the time of the spectroscopic campaign, we find that the photometric minimum occurs around ϕ = 0.54, which is close to the rotational phase when the spot which modulates radial velocity variations faces the observer (at ϕ = 0.6 in Fig. 6). We deduce that a large, cool spot is mainly responsible for both the photometric and radial velocity variability on a timescale of a few days. From the flat-bottom shape of Stokes I LSD profiles at phase 0.69, the relatively low amplitude of variability (∼5% in both photometric flux and V_r), and the sinusoidal shape of the K2 light curve, we further infer that the cool spot is presumably situated at a high latitude on the star. Indeed, the Doppler image reconstruction indicates the presence of a cool spot straddling the stellar pole and covering about 20% of the stellar surface.

In order to check if a cold spot can produce the observed modulation of radial velocity and photometric flux, we used SOAP2 (Dumusque et al. 2014) to generate synthetic variability curves. We used the stellar parameters of HQ Tau reported above and ran a set of simulations for a range of spot latitudes, spot sizes, and temperature differences between the spotted area and the stellar photosphere. The strong degeneracy between the spot location and its size does not allow us to strictly constrain these parameters. Nevertheless, we found a solution reproducing the amplitude of modulation of both radial velocity and flux with a spot located at a latitude of 60° and whose projected surface covers 30% of the visible hemisphere, qualitatively consistent with the Doppler map. These spot properties are not uncommon for T Tauri stars (e.g., Bouvier & Bertout 1989).

Most of the circumstellar diagnostics also seem to be modulated at the stellar rotation period, consistent with the expectations for magnetospheric accretion. The line profiles exhibit periodically modulated IPC components, whose maximum occurs around phase 0.9, and last for nearly half of the rotational cycle from phase 0.69 to 0.02. This is consistent with magnetospheric accretion funnel flows covering a significant azimuthal extension around the central star, as expected for a tilted large-scale magnetosphere. The shape of the profiles themselves are reminiscent of those computed from magnetospheric accretion models. For instance, the mean Hβ profile of HQ Tau is quite similar to that computed by Muzerolle et al. (2001) for Ṁ_acc = 10⁻⁹ M_⊙yr⁻¹ at an inclination of 60°, which shows a pronounced IPC component for a funnel flow temperature of 8000 K. The shape of the Hα profile is comparable to that computed by Lima et al. (2010) for this mass accretion rate.

The apparent mass accretion rate itself is modulated at the stellar rotational period. As it scales with the equivalent width of the line profiles, this is mostly the result of the periodic appearance of IPC components, which reduce the line flux. Hence, Ṁ_acc variations are primarily related to the geometric projection of the corotating accretion funnel flow. We do not have evidence for significant intrinsic variations of Ṁ_acc over a timescale of a couple of weeks, although we note that the depth of the IPC components slightly varies over this timescale. The deepest IPC components occur around ϕ = 0.9, which is significantly later than the epoch at which the cold spot that modulates the radial velocity faces the observer. If the cold spot marks the magnetic pole where the funnel flow is anchored, it might indicate that the upper part of the funnel flow that produces the redshifted absorptions is trailing. This may occur if the magnetospheric truncation radius is located beyond the disk’s corotation radius.

Such a phase delay has been previously reported for T Tauri stars, such as AA Tau (Bouvier et al. 2003) and V2129 Oph (Alencar et al. 2012). Unfortunately, the veiling is too weak in HQ Tau to trace the location of the accretion shock from the modulation of optical excess flux. Indeed, a weak veiling is not unexpected against a relatively bright photosphere (T_eff = 4997 K) at such a low Ṁ_acc. Similarly, we do not see a clear modulation of the longitudinal component of the magnetic field, presumably due to its complex topology, which prevents us from assigning a rotational phase to the magnetic pole. It may therefore be that the accretion shock at the base of the funnel flow is located at a slightly different longitude than the cold dark spot around the stellar pole, thus accounting for the IPC components appearing around ϕ = 0.9 without having to resort to a twisted magnetospheric accretion column.

Together with evidence for magnetospheric funnel flows, some line profiles reveal signatures of outflows. The Hα line profile exhibits a high velocity blueshifted absorption component, from about −200 to −150 km s⁻¹, which is modulated at the stellar rotational period. Interestingly, this component appears to be anti-correlated with the appearance of IPC features in the line profile. This is reminiscent of the “egg beater” model introduced by Johns & Basri (1995b) for SU Aur, where funnels flow and mass outflows occur on opposite azimuths at the disk inner edge: while accretion funnel flow is favored where the magnetic pole is inclined toward the disk, the situation is reversed at the opposite azimuth and favors interface or inner disk winds. This interpretation would be consistent with the overall variability reported here for the Hα profile of HQ Tau.

Transient mass loss episodes are seen on a timescale of a few days in the NaI D2 line profile. Several deep absorption components are seen over blueshifted velocities ranging from about −40 down to −120 km s⁻¹. These absorptions are not rotationally modulated. Instead, they are seen over a few consecutive days and then disappear altogether (see Fig. 11). These episodic mass outflows could conceivably be related to magnetic reconnections at the star-disk interface, following the magnetospheric inflation scenario put forward by Bouvier et al. (2003) for AA Tau: as the magnetospheric field lines inflate under the shear of differential rotation, they eventually open up and reconnect (Goodson et al. 1997). Such an inflationary cycle first induces a phase of reduced mass accretion onto the star during inflation, followed by a transient outflow during reconnection, and the restoration of funnel flow accretion.

A detailed analysis of the line profile variability may support the magnetospheric inflation scenario. The transient blueshifted absorptions are most conspicuous in the NaI D2 line profile over 2 rotational cycles, from J.D. 8060 to 8065 (see Fig. 11), and the central velocity of these components appear to slowly drift toward the line center over this timeframe. The IPC components seen in Balmer lines are the deepest on J.D. 8066 (ϕ = 4.88, 4.90), meaning right after the transient outflow episode. This would be consistent with the restoration of a magnetospheric accretion funnel after reconnection of the inflated magnetosphere. In contrast, at the start of the transient outflow episode, on J.D. 8061 (ϕ = 2.85), the IPC profiles are the weakest, consistent with an inflated state for the magnetosphere. As discussed in Bouvier et al. (2003) and Alencar et al. (2018), the timescale for the magnetospheric inflation cycle is expected to be several times the rotational period of the star (Zanni & Ferreira 2013), which is presumably the reason why it is rarely observed in spectroscopic time series extending over only a couple of weeks.

Although we detect only relatively weak Stokes V signatures in the LSD spectra of HQ Tau, they are sufficient to derive the longitudinal component of the magnetic field and reconstruct surface magnetic maps. We do not detect rotational modulation of the longitudinal magnetic field, unlike what is usually seen in T Tauri stars (e.g., Donati et al. 2019, 2020). The rotational modulation of B_l in photospheric lines is often complex, except in the dipolar magnetic field case. It may also be that the weakness of the signatures and the complexity of the field topology at the stellar surface combine to hide the modulation. As HQ Tau has already developed a significant radiative core, we do not expect a surface magnetic field dominated by a dipolar component (Gregory et al. 2012).

The full ZDI analysis suggests a weak dipolar component, inclined by 31.5° from the rotational axis, and amounting to about 63 G. However, cancellation of opposite polarities on small spatial scales at the stellar surface may affect the total magnetic energy of the star. Also, strong field regions may be hidden in dark stellar spots. Therefore, the total magnetic energy we measure with ZDI on photospheric LSD profiles is to be taken as a lower limit (Lavail et al. 2019; Sokal et al. 2020). Indeed, the magnetic field strength derived from Stokes V signatures in the CaII and (or) HeI emission line profiles of T Tauri stars usually exceeds that deduced from LSD profiles by at least a factor of 5 to 6 because they probe small and localised highly magnetised regions on the surface of the star (Donati et al. 2019, 2020). In HQ Tau, the only emission line available to investigate additional magnetic tracers is CaII IRT. As this line is formed at least partly in the accretion funnel flow, it is a potential probe of the specific field connection between the star and the disk. From the lack of detectable Stokes V signal, we could merely derived a 3σ upper limit of ∼450 G for the longitudinal field component in the emission line core.

The two intermediate mass T Tauri stars previously imaged in this way, CR Cha (P_rot = 2.3 d) and CV Cha (P_rot = 4.4 d) by Hussain et al. (2009), have similar spectral types and masses as HQ Tau but have larger radiative cores (M_rad/M_⋆ ≈ 0.65 and 0.92 respectively) and are therefore probably further along in their PMS evolution. In fact, HQ Tau and CR Cha also have very similar rotation periods and mass accretion rates, with Ṁ_acc = 2 × 10⁻⁹ M_⊙yr⁻¹ for CR Cha (Nisini et al. 2018). It is interesting to note that there are strong similarities between all three stars in their surface activity, as shown by comparing their brightness maps and their large scale magnetic field distributions. All three stars show a large cool spot near their poles, with CR Cha showing an almost-identical off-centre-polar cap like HQ Tau. The large scale magnetic fields of all three systems are all complex and non-axisymmetric.

Are our results consistent with magnetospheric accretion being at work in HQ Tau? Unfortunately, due to the lack of Zeeman signatures in emission lines, we cannot derive from our data the strength of the dipolar magnetic field in the acretion shock, which would allow us to compute the magnetospheric truncation radius, r_mag, according to Bessolaz et al. (2008)’s prescription:

$$\begin{array}{c}\hfill \frac{r_{\mathrm{mag}}}{R_{\star }}=2m_s^{2/7}{B}_{\star }^{4/7}{\dot{M}}_{\mathrm{acc}}^{-2/7}{M}_{\star }^{-1/7}{R}_{\star }^{5/7},\end{array}$$(5)

where m_s ≈ 1, B_⋆ is the equatorial magnetic field strength, Ṁ_acc is the mass accretion rate, M_⋆ the stellar mass, and R_⋆ the stellar radius, respectively in units of 140 G, 10⁻⁸ M_⊙ yr⁻¹, 0.8 M_⊙, and 2 R_⊙. Instead, we may assume that the truncation radius is located close to the corotation radius, r_cor = 3.55 ± 0.35 R_⋆. We then derive from Eq. (5) a magnetic field strength of ∼120 G at the equator, which translates to 240 G at the pole for the dipolar component (and to ∼3 G at the truncation radius, following Gregory 2011). Although we cannot measure this component directly, we note that the value required at the stellar surface is not inconsistent with the large-scale ZDI magnetic maps derived above nor with the upper limit we set on the CaII line field strength.

An independent support to the magnetospheric accretion scenario comes from the maximum velocity of IPC components measured in the emission line profiles, ${\mathit{v}}_{\mathrm{IPC}}^{\max }\sim 330\pm 19$ km s⁻¹. We measured this velocity by fitting the red part of Hβ’s IPC components by a straight line. The location where the line joins the continuum yields to ${\mathit{v}}_{\mathrm{IPC}}^{\max }$. This provides an estimate of the free-fall velocity of the accreted material projected onto the line of sight. Accounting for projection effects that combines the star’s inclination (i  ∼  60°, as the average between the rotational and ZDI estimates above) and the magnetic obliquity (β  =  31.5°), we have ${\mathit{v}}_{\mathrm{IPC}}^{\max }\simeq {\mathit{v}}_{\mathrm{ff}}cos(i-\beta )$, which yields v_ff  =  376 km s⁻¹. Assuming free fall accretion from r_cor to R_⋆, the material would hit the stellar surface at a velocity v_ff = (2GM_⋆/R_⋆)^1/2(1 − R_⋆/R_cor)^1/2 = 424 km s⁻¹, which is consistent with the estimate obtained from the maximum redshifted velocity of the IPC components, and thus supports a magnetospheric truncation radius being located close to the disk’s corotation radius.

Finally, while we observed HQ Tau during a relatively quiescent phase of variability, we note that its behavior on longer timescales may be more complex. In particular, the deep and long lasting UXOr events reported for the system by Rodriguez et al. (2017) must have another origin than the variability described here. The high inclination we derive for the system, $i={75}_{-17}^{+15}$ deg., could favor transient circumstellar extinction events, such as UXOr and (or) AA Tau type. One of these deep faintening events occurred just before our spectroscopic observations, reaching nearly 1.5 magnitudes and lasting for a month, with significant intra-variability. The system became much redder during the event, consistent with obscuration by circumstellar dust. We may envision two related scenarios to explain this additional component to the system’s variability. One is that it might result from a sudden change in the vertical scale-height of the inner disk. As the inner disk edge lies at a distance of only ∼3 R_⋆ from the stellar surface, and the system is seen at high inclination, a puffed-up inner disk could conceivably obscure the star for the duration of the instability (Turner et al. 2010). Another possibility is a change in the location of the magnetospheric truncation radius relative to the dust sublimation radius. From the expression of Monnier & Millan-Gabet (2002), we derive a sublimation radius r_sub = 5.55 ± 1.00 R_⋆, using Q_R = 1 and T_s = 1500 K. As r_sub is larger then r_mag, there is no dust at the magnetospheric truncation radius, which accounts for the lack of dipper-like occultations at the time of our (and K2) observations. However, should the truncation radius increase, following either an increase in the magnetic field strength (by a factor of 3) or a decrease in the mass accretion rate (by a factor of 8), to eventually reach the sublimation radius, a dusty disk warp would result and could induce occultation events in this highly inclined system. We therefore suggest that the episodic UXOr events regularly observed in this system are due to the combination of close-in circumstellar dusty material around the central star and the specific geometry under which the system is seen.

5. Conclusions

Following the K2 Taurus campaign, we monitored the young, intermediate mass HQ Tau system with spectropolarimetry in order to investigate the accretion and ejection processes on a timescale of days. The results of this campaign provide clear diagnostics of magnetically mediated accretion occurring in the system. We observed redshifted absorptions in emission line profiles that are periodically modulated at the stellar rotational period, indicative of magnetospheric funnel flows passing through the line of sight. This is expected from a global dipolar topology of the stellar magnetic field at few stellar radii, despite the fact that, at the stellar surface, the Stokes V signatures indicate a more complex topology. While the star hosts a modest dipolar magnetic field component, the reduced mass accretion rate and the star’s rapid rotation result in the magnetospheric truncation radius being close to the disk’s corotation radius, as observed in most T Tauri systems investigated so far. Concurrent with accretion diagnostics, spectral signatures of outflows are clearly seen in line profiles, some being modulated by stellar rotation, others being transient phenomena, possibly resulting from instabilities at the star-disk magnetospheric boundary. We thus conclude that intermediate-mass pre-main sequence systems may undergo similar accretion and ejection processes as lower mass T Tauri stars, including stable magnetospheric accretion funnel flows, in spite of moderate magnetic field strength. Systems seen at high inclination also experience longer term faintening events, which result from the occultation of the inner system by circumstellar dusty material.

HQ Tau is a fascinating pre-main sequence system, reported here to possibly be a short period spectroscopic binary whose orbit remains to be determined. It exhibits a mix of variability behaviors, including accretion signatures modulated on a timescale of a few days as seen in lower mass T Tauri stars, as well as longer term UXOr events more typical of Herbig Ae stars. The richness of variability patterns this intermediate-mass young star displays warrants additional multi-wavelength studies of the system and its circumstellar environment on the long term.

Acknowledgments

We thank an anonymous referee whose comments improved the content of this paper. We thank Lynne Hillenbrand and Hervé Bouy for taking additional high resolution spectra of HQ Tau that helped us assessing its binary nature. We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. 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. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 742095; SPIDI: Star-Planets-Inner Disk- Interactions; http://www.spidi-eu.org). This project was funded in part by INSU/CNRS Programme National de Physique Stellaire and Observatoire de Grenoble Labex OSUG2020. S.H.P. Alencar acknowledges financial support from CNPq, CAPES and Fapemig.

References

Appendix A: Raw emission line profiles

The average Hα, Hβ, and CaII IRT line profiles of HQ Tau are shown in Fig. A.1. These figures show the lines before subtracting the photospheric template profiles.

Fig. A.1. From left to right: mean Hα, Hβ, and CaII IRT HQ Tau’s line profiles (blue) overplotted on those of the photospheric template Melotte 25-151 (dash dotted black). The template has been corrected for the radial velocity shift and rotationally broadened to the vsini of HQ Tau.

All Tables

All Figures

Fig. 1. Top: HQ Tau’s detrended Kepler K2 light curve. Low frequency variations have been removed from the original light curve shown in Fig. 2. The spike appearing at J.D. 2 457 890.0 is presumably an instrumental defect. Middle: HQ Tau’s light curve folded in phase with P = 2.424 d. The origin of phase is taken as JD 2 457 823.81 to correspond to the epoch of maximum brightness. Bottom: lomb-Scargle periodogram obtained from HQ Tau’s Kepler K2 light curve.

In the text

Fig. 2. HQ Tau’s V-band light curve. CrAO measurements are shown as filled grey circles, the K2 light curve as a continuous line, AAVSO measurements as filled triangles, and the ASAS-SN rescaled g- and V-band datasets as crosses and open circles, respectively. The upper panel shows the light curve over 2 observing seasons while the lower one displays the photometric variability of the source around the CFHT/ESPaDOnS observations.

In the text

Fig. 3. Color-magnitude plots from CrAO VRJIJ photometry. Left: (V − R)J vs. V. Right: (V − I)J vs. V. The source gets redder when fainter. The dashed line is the expected reddening line for ISM extinction. The photometric measurements are color-coded by Julian date to distinguish the 2 observing seasons: lighter points correspond to the first epoch, darker points to the second one (see Fig. 2). We note that the deep faintening event of about 1.2 mag in the V-band is accompanied by a strong reddening of the system.

In the text

	Fig. 4. Upper six rows: spectral windows selected from HQ Tau’s mean ESPaDOnS spectrum to derive stellar parameters. Bottom row: as an example, HQ Tau’s mean spectrum (blue) fitted with a ZEEMAN synthetic spectrum (orange) in the 547–557 nm window.
In the text

	Fig. 5. HQ Tau’s position in the HR diagram. The red cross illustrates HQ Tau’s position with corresponding uncertainties. The black curves are the evolutionary tracks from CESTAM PMS models, with the corresponding mass indicated on the left. The blue curves show the position where the radiative core of the star reaches 75%, 50% and 25% of the stellar mass.
In the text

	Fig. 6. From top to bottom: radial velocity, Hα line flux, and average longitudinal magnetic field as a function of Julian date (left) and rotational phase (right). The sinusoidal fit, forced at the K2 photometric period, to radial velocity variations is shown as a black dash-dotted curve in the top panels. The Bl = 0 G line is drawn for reference in the bottom panels.
In the text

Fig. 7. Top: residual Hα (left), Hβ (middle), and CaII 854.2 nm (right) emission lines of HQ Tau spectra. Bottom: corresponding 2D periodograms with the mean profile superimposed as a white curve. The maximum power of the periodogram is in yellow and the minimum in dark blue. The red and blue wings of Hα and the red wing of Hβ and CaII, as well as the center of the 3 line profiles, exhibit a periodic modulation at f ≈ 0.4 d−1, corresponding to the rotational period of the star. The power seen at f ≈ 0.6 d−1 is the 1-day alias induced by the observational window.

In the text

	Fig. 8. Residual Hα (left panels) and Hβ (right panels) profiles sorted by rotational phase and by day of observation. The successive rotational cycles are shown by different colors, as in Fig. 6.
In the text

	Fig. 9. Correlation matrices of the line profiles. (a) Hα vs. Hα, (b) Hβ vs. Hβ, (c) Hβ vs. Hα, (d) CaII vs. CaII, (e) Hα vs. CaII, (f) NaD vs. NaD.
In the text

	Fig. 10. Top: NaI D2 line appearing in the 14 spectra of HQ Tau. The black dash-dotted line is a synthetic NaI D2 line produced using the ZEEMAN code with HQ Tau’s parameters. Bottom: 2D periodogram of the NaI D2 line.
In the text

	Fig. 11. NaI D2 line profiles sorted by rotational phase (left) and by day of observation (right). The successive rotational cycles are shown by different colors, as in Fig. 6.
In the text

	Fig. 12. HQ Tau’s Stokes I(left), V(middle), and Null (right) profiles (black) ordered by increasing fractional phase, and the fit produced by the ZDI analysis (red). The fractional number next to each profile indicates the rotational phase.
In the text

	Fig. 13. Upper panel: surface brightness map of HQ Tau. The color scale goes from dark brown to light yellow, for dark to bright features. Lower panels, from top to bottom: radial, azimuthal, and meridional surface magnetic field distribution. Ticks on top of the panels indicate the phase of observations.
In the text

	Fig. 14. Radial velocity of HQ Tau as a function of date of observation. The error bars are indicated, and can be smaller than the symbol size.
In the text

	Fig. A.1. From left to right: mean Hα, Hβ, and CaII IRT HQ Tau’s line profiles (blue) overplotted on those of the photospheric template Melotte 25-151 (dash dotted black). The template has been corrected for the radial velocity shift and rotationally broadened to the vsini of HQ Tau.
In the text