© The Authors 2023

1. Introduction

Accretion and ejection processes are at the origin of most of the peculiar properties of young stellar systems. The structure and dynamics of the accretion flows within the disk and from the inner disk to the star, as well as the properties of the multiple outflows arising from the disk, from the star-disk interface, and from the stellar surface, remain to be fully deciphered, however. Low-mass pre-main-sequence stars, the so-called T Tauri stars (TTS), accrete from their circumstellar disks for a few million years, while contemporaneous planet formation impacts the disk structure and evolution. In the inner regions of the system, the disk is disrupted by the strong stellar magnetosphere that channels the accretion flow toward the star along magnetic field lines (see, e.g., the review by Hartmann et al. 2016). Thus, accretion funnel flows develop that connect the inner disk to the stellar surface, where the material is accreted at nearly free-fall velocity and is eventually halted in a strong accretion shock. Simultaneously, outflows are produced at the star-disk interface close to the magnetospheric truncation radius through the inflation and reconnection of magnetic field lines that are twisted by differential rotation (e.g., Zanni & Ferreira 2013). Ultimately, the release of gravitational energy delivered by the accretion process may trigger accretion-powered stellar winds (Matt & Pudritz 2005). The torque balance between accretion and ejection processes is a central issue for understanding the spin evolution of young stars (e.g., Pantolmos et al. 2020; Ireland et al. 2021).

The star-disk interaction takes place on a distance of a few stellar radii (e.g., Bessolaz et al. 2008), that is, on a scale of about 0.1 au or smaller. MHD models developed by several groups predict the structure and dynamics of the magnetospheric accretion region and associated outflows (see, e.g., the review by Romanova & Owocki 2015). Observationally, two main directions have been explored so far to investigate the properties of this region. On one hand, monitoring the spectroscopic and photometric variability of the system over a few rotational periods, that is, typically over a few weeks, allows identifying the signature of funnel flows, hot spots, and outflows, and relating them to the strength and topology of the surface magnetic field that is measured from spectropolarimetry (e.g., Pouilly et al. 2020, 2021; Bouvier et al. 2020a; Donati et al. 2019, 2020a; Alencar et al. 2018). On the other hand, a direct approach attempts to spatially resolve the star-disk interaction region on a scale of a few milliarcsecond on the sky, using long-baseline near-infrared interferometry (e.g., Eisner et al. 2014; Gravity Collaboration 2020; Bouvier et al. 2020b). Both approaches have been successful in mapping the inner region of accreting systems and have provided strong support to the magnetospheric accretion scenario and its MHD modeling. Following previous studies, we report here the results of a new observing campaign devoted to the young stellar system GM Aur.

GM Aur (RA = 04h55, Dec = +30°21, V = 12.1 mag) is a solar-type pre-main-sequence star located in the Taurus-Auriga molecular cloud at a distance of 157.9 ± 1.2 pc (Gaia Collaboration 2021). This classical T Tauri star (cTTS) has a spectral type K6 (Herczeg & Hillenbrand 2014) and is surrounded by a circumstellar disk from which it actively accretes material at a rate of 0.6–2.0 × 10⁻⁸ M_⊙ yr⁻¹ (Robinson & Espaillat 2019). Based on its spectral energy distribution, which exhibits a small near-infrared excess compared to a significant mid-infrared one, the system has long been suspected to be in a transitional stage, that is, that it is surrounded by a disk whose inner regions are relatively devoid of matter (Strom et al. 1989). High-resolution ALMA images of the circumstellar disk indeed reveal that it is highly structured. The large-scale disk, inclined at ∼53° on the line of sight, features a large inner dust cavity extending over ∼35–40 au and a succession of annular gaps and dusty rings on a wider scale up to 200 au (Macías et al. 2018; Huang et al. 2020). Much closer to the central star, long-baseline VLTI/GRAVITY interferometric observations unveil a compact dusty disk, whose inner edge was recently reported to be located at r_in = 0.013${}_{-0.008}^{+0.015}$ au from the central star (Bohn et al. 2022) and that extends over at least a few 0.1 au (Akeson et al. 2005) and possibly up to 6.6 au (Varga et al. 2018; Woitke et al. 2019). The gaseous component of the inner disk has been detected from CO 4.7 micron emission down to 0.5 ± 0.2 au (Salyk et al. 2009). The inclination and position angle of the major axis of the inner dusty disk (i = 68° ${}_{-28}^{+16}$, PA = 37° ${}_{-22}^{+31}$) are found to be consistent with those of the outer disk, which suggests that the inner and outer disks are aligned (Bohn et al. 2022).

In an attempt to decipher the physical processes at work at the heart of the system, GM Aur has been the subject of several multiwavelength monitoring campaigns. The long-term light curve presented by Grankin et al. (2007) over the period 1986–1995 exhibits relatively low-level variability, with a V-band magnitude ranging from 11.74 to 12.35 mag. Photometric variations are modulated by surface spots at the stellar rotation period of 6.0–6.1 days (Percy et al. 2010; Artemenko et al. 2012). Ingleby et al. (2015) reported variability over the full wavelength range from the far-UV to the near-infrared, which they attributed in part to an accretion rate that varies by about a factor of 2–3 on a timescale of months, and for another part to dust inhomogeneities that are located in the inner disk close to the truncation radius. Variations in the mass accretion rate of similar amplitude have also been reported on a shorter timescale of about a week by Robinson & Espaillat (2019), and a connection between mass loss and mass accretion has been further suggested by Espaillat et al. (2019). McGinnis et al. (2020) presented the results of a high-resolution optical spectroscopic monitoring campaign performed on a timescale of a week that illustrated the variability of the Hα, Hβ, and HeI emission line profiles of the system. From the measured radial velocity variations of the HeI 5876 Å line profile, whose narrow component traces the accretion shock, they deduced that GM Aur accretes material from its circumstellar disk through an inclined magnetosphere, whose axis is tilted by about 13° relative to the stellar rotational axis. GM Aur indeed harbors a strong surface magnetic field, with a mean value of 2.2 kG (Johns-Krull 2007; Symington et al. 2005). Finally, from a recent multiwavelength X-UV-optical campaign, Espaillat et al. (2021) reported evidence for a transverse density stratification within the accretion shock at the base of the magnetic funnel flow.

We report here the results of a new coordinated monitoring campaign on GM Aur that combines high-resolution optical spectroscopy and near-infrared spectropolarimetry, multiwavelength optical and near-infrared photometry, and long-term low-resolution near-infrared spectroscopy. Part of the observations have been obtained simultaneously over a timescale of a few weeks, while the total duration of the campaign amounted to six months. The goal of the campaign was to investigate the physical processes that cause variability in GM Aur on a scale of a few stellar radii, and in particular, to constrain the structure and dynamics of the magnetospheric accretion flow from the inner disk to the star. We devised a long-term campaign in order to be able to probe various timescales, from days to months, and obtain a sufficiently long temporal baseline to investigate the relation between accretion and ejection processes on small spatial scales from the stellar surface to the inner disk regions.

The campaign whose results are reported here took place in the framework of a larger project led by the ODYSSEUS team¹ (see Espaillat et al. 2022), which uses the Hubble UV Legacy Library of Young Stars as Essential Standards program (ULLYSES², Roman-Duval et al. 2020), on HST Director’s Discretionary time, to monitor a sample of T Tauri stars in the UV domain, which includes GM Aur. Additional follow-up observations were acquired for this project at ESO in the framework of the PENELLOPE Large Program³ (Manara et al. 2021).

In Sect. 2 we describe the observational techniques we implemented to perform the campaign. In Sect. 3 we derive the properties of the system and analyze its photometric and spectroscopic variability over timescales from days to months, including veiling measurements and emission line profiles. We infer the global structure of the magnetospheric accretion flow from the observed variability and characterize associated outflows. In Sect. 4 we discuss the dynamics of the accretion and ejection structure and show that short-lived episodic outflows coexist with a stable magnetospheric accretion pattern. In Sect. 5 we conclude on the ability of multiwavelength, multi-technique coordinated observational campaigns to unveil the physical processes at work in young stellar systems at the sub-au scale.

2. Observations

In this section, we describe the acquisition and data-reduction processes of photometric, spectroscopic, and spectropolarimetric datasets obtained during the large-scale campaign we performed on the cTTS GM Aur from September 6, 2021, to March 8, 2022, using CFHT/SPIRou, OHP/SOPHIE, ESO/ExTrA, LCOGT, and ESO/REM. A summary plot of the GM Aur observing campaign reported here is provided in Fig. 1.

Fig. 1. Temporal sampling of the GM Aur campaign from September 6, 2021, to March 8, 2022. Bottom: g′-band light curve from LCOGT (black dots) and from REM (magenta crosses). The REM g′-band magnitudes are offset in this figure by +0.1 mag to match the LCOGT measurements. The mean photometric error on both the LCOGT and REM g′-band measurements is 0.025 mag. Vertical arrows: the vertical arrows show the dates of OHP/SOPHIE (red), CFHT/SPIRou (blue), and ESO/ExTrA (black) observations. The core of the campaign took place during October 2021 with contemporaneous measurements from the five instruments.

2.1. LCOGT: Multiwavelength optical photometry

GM Aur was observed at Las Cumbres Observatory Global Network (LCOGT, Brown et al. 2013) from September 6 to December 30, 2021. We acquired 850 images in the Sloan u′g′r′i′ filters over two runs with a sub-day cadence (LCO2021B-001, PI L. Rebull; CLN2021B-003, PI A. Bayo). The u′ images were obtained with the Sinistro 1m telescopes of the network using an exposure time of 180 seconds and reading the 2K × 2K central window of the detector with a 2 × 2 binning, resulting in a 13 × 13 arcmin field of view on the sky. The g′r′i′ images were obtained with the 0.4m SBIG telescopes, offering a field of view of 29.2 × 19.5 arcmin, with exposure times of 60, 20, and 20 seconds, respectively. We retrieved the BANZAI-reduced images from the LCOGT archive service and the noncalibrated photometric catalogs provided in the image headers for all detected stars in the field.

In order to compute differential photometry, we considered two stars, HD 282625 and HD 282626, both located within 3 arcmin of GM Aur. The first star was used as a reference star to calibrate the differential light curve, and the second was used as a check star to assess that these are nonvariable sources. These field stars have spectral types F2 and F5, respectively, and are only slightly brighter than GM Aur. We confirmed that the two stars are nonvariable from their differential light curves, and we deduced a mean rms photometric error of 0.025 mag in the u′g′r′ filters and 0.033 mag in the i′ filter. We proceeded to compute the differential light curve between GM Aur and HD 282625 in all four filters. We adopted the mean magnitude of HD 282625 listed in the APASS and Pan-STARRS surveys, namely g′ = 11.331 mag, r′ = 10.916 mag, and i′ = 10.756 mag to calibrate the GM Aur light curve in the g′r′i′ filters to within an accuracy of 0.02 mag.

We were not able to find an estimate of the u′-band magnitude for the comparison stars in the literature. Instead, we assumed the intrinsic (u′−g′) colors of an F2 star (Covey et al. 2007; Kraus & Hillenbrand 2007) for HD 282625, to which we applied interstellar reddening. From the observed versus intrinsic color indices of the HD 282625 (g′−r′) and (r′−i′) bands, we derived A_V = 0.8 ± 0.1fmag, using the R = 3.1 reddening law from Fiorucci & Munari (2003). This procedure yielded an estimate of the reddened (u′−g′) color of the comparison star from which we derived its u′-band magnitude. The photometric calibration of the GM Aur u′-band light curve is thus relatively indirect and probably not accurate to better than 0.1 mag⁴.

2.2. REM: Optical and near infrared photometry

Observations were performed with the 60 cm robotic REM telescope located at the ESO La Silla Observatory (Chile), on 15 nights from JD 2 459 497 to JD 2 459 520 (October 9 to November 1, 2021). By means of a dichroic, REM simultaneously feeds two cameras at the two Nasmyth focal stations, one camera for the near-infrared (REMIR), and the other for the optical (ROSS2). The cameras have nearly the same field of view of about 10′×10′ and use wide-band filters (J, H, and K′ for REMIR and Sloan/SDSS g′, r′, i′, and z′ for ROSS2).

Exposure times were 60 s for ROSS2, which simultaneously acquires images in the four Sloan bands, and five ditherings of 3 s each were adopted for each filter of REMIR. For the ROSS2 camera, we generated master flats using the twilight flat-fields taken during the observing run, which are available in the REM archive. The latter were used to correct for pixel-to-pixel sensitivity variations, as well as for the vignetting and illumination of the field of view. After subtracting the dark-frame, each scientific image was divided by the proper master-flat, depending on the filter. The prereduction of the REMIR images is automatically done by the AQuA pipeline (Testa et al. 2004), and the coadded and sky-subtracted frames, resulting from five individual ditherings, are made available to the observer.

The adopted comparison stars are reported in Table A.1 along with their griz (Tonry et al. 2018) and JHK′ (Cutri et al. 2003) magnitudes. Aperture photometry for all the stars listed in Table A.1 was performed with DAOPHOT by using the IDL⁵ routine APER. For each frame and filter, we used the instrumental magnitudes of the stars listed in Table A.1 to generate an artificial comparison, weighting them with the flux corresponding to their standard magnitude in a way similar to the ensemble photometry. This procedure also allowed us to evaluate a standard error based on the differences between the magnitudes calculated with different comparison stars.

The optical photometry gathered at REM is listed in Table A.2. In the common g′r′i′ bands, we found it to agree well with that obtained at LCOGT with a tighter sampling rate, and therefore, we did not use it further in the analysis below. The individual JHK′ measurements are listed in Table 1. The median errors on JHK′ measurements are 0.05, 0.05, and 0.06 mag, respectively. However, some measurements were affected by the nearby bright moon around JD 2 459 512, and had an error larger than 0.1 mag. We chose to discard these measurements. We eventually derived the following values for the median near-infrared magnitudes of the GM Aur system and their rms variations at the time of the observations: J = 9.41 ± 0.10 mag, H = 8.71 ± 0.03 mag, and K′ = 8.40 ± 0.07 mag.

2.3. OHP SOPHIE: High-resolution optical spectroscopy

Observations were carried out from October 12 to 29, 2021, at Observatoire de Haute-Provence using the fiber-fed SOPHIE spectrograph (Perruchot et al. 2008) in high-efficiency mode, which delivers a spectral resolution of R ∼ 40 000 over the wavelength range 387–694 nm. We obtained 15 spectra over 18 nights, with an exposure time of 3600 s, yielding a signal-to-noise ratio ranging from 42 to 67 at 600 nm. The raw spectra were fully reduced at the telescope by the SOPHIE real-time pipeline (Bouchy et al. 2009). The data products include a resampled 1D spectrum with a constant wavelength step of 0.01 Å, corrected for barycentric radial velocity, an order-by-order estimate of the signal-to-noise ratio, and a measurement of the source radial velocity, V_r, and projected rotational velocity, v sin i. The latter two quantities are derived from a cross-correlation analysis of nearly 7000 spectral lines between the observed spectrum and a K5 spectral mask template (e.g., Boisse et al. 2010). We list the values of these parameters in Table 2. The mean formal error provided by the SOPHIE pipe-line on the V_r measurement is 0.013 km s⁻¹. This accuracy is well suited to investigating the significantly larger amplitude of photospheric line profile variability induced by surface spots and/or accretion flows in young stars (e.g., Petrov et al. 2001). No error is provided by the pipeline for the v sin i measurements, and we assumed that an upper limit is given by the rms deviation of the individual measurements, excluding JD 2 459 512 (see below), namely 0.38 km s⁻¹.

The SOPHIE spectrograph includes a second fiber that simultaneously records the spectrum of the nearby sky. Inspection of the cross-correlation function (CCF) of the sky fiber with a synthetic mask of spectral type G2 revealed that the signature of the moon becomes apparent at the expected barycentric Earth radial velocity from JD 2 459 505 onward because the growing moon approaches the target. The lunar contamination culminates on JD 2 459 512, as the bright moon is located about 10 degrees away from GM Aur, which explains the discrepant values measured for V_r and v sin i on this date. Table 2 suggests that except for JD 2 459 512, the contamination of the CCF by the moon only marginally impacts the V_r and v sin i measurements. However, to be conservative, we only considered the V_r and v sin i measurements obtained from the first six spectra of the observing run for the subsequent analysis, from JD 2 459 499 to JD 2 459 504, where no lunar contamination is present.

The journal of observations is given in Table 2. It lists the Julian date, the signal-to-noise ratio of individual spectra at 600 nm, the radial and rotational velocities derived from each spectrum, and the bisector span computed from the cross-correlation function (Queloz et al. 2001).

2.4. CFHT SPIRou: Near-infrared spectropolarimetry

Near-infrared spectropolarimetric observations of GM Aur were performed at CFHT using the SPIRou near-infrared spectropolarimeter. It has a spectral range covering from 0.95 to 2.50 μm in a single exposure at a spectral resolution of 70 000 (Donati et al. 2020b). The observations were completed in the framework of the CFHT SPIRou Legacy Survey over four observing runs extending from September 15 to December 18, 2021. Each monthly run was scheduled around the full moon, and we aimed at obtaining one spectrum per night during the run. An additional single spectrum was obtained on January 6, 2022. We thus gathered 34 spectra, whose temporal sampling is shown in Fig. 1.

Each spectrum consists of four polarimetric sub-exposures⁶ for a total integration time of 2,200 s. Individual exposures were combined to yield a single spectrum with a signal-to-noise ratio ranging from ∼60 to 100. In one instance, on JD 2 459 503.08, the polarimetric sequence was aborted, and the spectrum consists of a single sub-exposure. The raw data were reduced within the SPIRou consortium, using version V6.132 of the APERO pipeline (Cook et al. 2022). Spectra were cross-correlated with a K2 spectral mask template over about 6,700 spectral lines, and the radial velocity of the object was derived with sub-km s⁻¹ accuracy (0.08 km s⁻¹ mean rms uncertainty) by fitting a Gaussian to the resulting CCF. The median V_r amounts to 14.65 ± 0.27 km s⁻¹. Using TWA 9A, a WTTS of spectral type K5, as a template, the veiling was derived in the JHK bands following the procedure described in Sousa et al. (2023) The median rms errors on veiling measurements in the JHK bands are 0.011, 0.025, and 0.027, respectively.

The journal of observations is presented in Table 3. It lists the Julian date, the signal-to-noise ratio at 2.16 μm, the photospheric radial velocity and its uncertainty, and the veiling in the JHK bands.

2.5. ESO ExTrA: Low-resolution near-infrared spectroscopy

The ExTrA facility (Bonfils et al. 2015), located at La Silla Observatory in Chile, consists of three 60 cm telescopes and a single near-infrared (0.88–1.55 μm) fibre-fed spectrograph. We observed GM Aur on 89 nights between October 13, 2021, and March 8, 2022, using either one telescope (21 nights) or two telescopes simultaneously (68 nights). Five fiber units are located at the focal plane of each telescope, each consisting of two 8″ aperture fibers. One fiber is used to observe a star and the other is used to observe the nearby sky background. We observed GM Aur with one fiber unit and used another fiber unit to simultaneously observe 2MASS J04535474+3021441 (J = 8.450 ± 0.027 mag) as a comparison star to compute differential photometry. We used the higher-resolution mode of the spectrograph (R ∼ 200) and 300-second exposures. We obtained between 1 and 30 exposures per night for a total of 1898 spectra with a median signal-to-noise ratio of 105 at 1.05 μm for GM Aur. The ExTrA data were corrected for dark current, extracted using the flat-field, corrected for sky background emission, and were wavelength calibrated using custom data reduction software. Median spectra of GM Aur were computed for each night and telescope, yielding a total of 157 spectra with a median signal-to-noise ratio of 179 and a standard deviation of 62.

We computed differential photometry of GM Aur relative to the comparison star by integrating the individual ExTrA spectra over the J-filter passband. The UKIRT-WFCAM filter transmission curves were retrieved from the SVO Filter Profile Service⁷ (Rodrigo et al. 2012; Rodrigo & Solano 2020). We multiplied each corrected individual spectrum of GM Aur and the comparison star by the filter transmission curve, integrated the flux, and obtained a magnitude difference from the flux ratio. We computed a differential magnitude measurement for each night as the mean and standard deviation of the individual measurements taken on that night. We then derived the J-band magnitude of GM Aur from the 2MASS magnitude of the comparison star. We obtained a median value and 68.3% confidence interval of J = 9.417 ± 0.061 mag for the ExTrA observations. The results are listed in Table C.1.

3. Results

3.1. Multicolor photometry

The LCOGT u′g′r′i′ light curves of GM Aur span nearly four months and are shown in Fig. 2. A TESS light curve extending over 50 days, from JD 2 459 474 to JD 2 459 524, is also shown. We derive a mean magnitude and rms variability of u′ = 14.32 ± 0.35 mag, g′ = 12.77 ± 0.17 mag, r′ = 11.64 ± 0.11 mag, and i′ = 11.21 ± 0.09 mag, respectively. The variability amplitude is much larger than the photometric errors in all bands; it amounts to 0.025 mag in the g′r′i′ bands and 0.033 mag in the u′ band, and it decreases toward longer wavelengths, from u′ to i′. We ran two period-search algorithms on the light curves: a CLEAN periodogram analysis (Roberts et al. 1987), which is conceptually similar to a Fourier transform, and the String-Length method (Dworetsky 1983), which finds the period that minimizes the average distance between consecutive points in the phased light curve. Both yielded the same period at all wavelengths, namely P = 6.04 ± 0.15 days, with the uncertainty being estimated from the standard deviation of a Gaussian fitted to the periodogram peak. The results of the period search are shown in Fig. 3. The light curves folded in phase at the 6.04 d period are shown in Fig. 2. They clearly show the modulation of the brightness level at this period, particularly in the u′ band, where the amplitude is largest. We ascribe this low-level modulation to surface spots and the P = 6.04 d period to stellar rotation. This estimate agrees within the error bars with the previously reported periods for GM Aur, namely P = 6.1 d by Percy et al. (2010) and P = 6.02 d by Artemenko et al. (2012). We therefore adopt the following ephemeris:

$$\begin{array}{c}\hfill \mathrm{HJD}(\mathrm{d})=2459460.80+6.04\times E,\end{array}$$(1)

Fig. 2. GM Aur light curves. Top: GM AUr’s u′g′r′i′ light curves from LCOGT observations that extended over nearly four months. The mean photometric error is 0.025 mag in the g′r′i′ bands and 0.033 mag in the u′ band. A scaled TESS light curve (black) obtained contemporaneously is overplotted onto the LCOGT r′-band light curve. Middle: LCOGT light curves folded in phase with a period of 6.04 days with the ephemeris of Eq. (1). Bottom: same as above, with a color scale for data points that reflects the Julian date. The amplitude of the light variations appears to increase slightly toward the end of the observing run. In all panels, the error bars on the measurements are smaller than the symbol size.

Fig. 3. Period search results. Top: CLEAN periodogram analysis of the u′g′r′i′ light curves. A peak occurs at the frequency of 0.166 day−1, which corresponds to a period of 6.04 ± 0.15 days. Bottom: string-length analysis of the u′g′r′i′ light curves. A clear minimum appears for a period of 6.03 ± 0.15 days.

which defines the rotational phase Φ_rot = (HJD-2 459 460.80)/P_rot (modulo P_rot), where phase zero (Φ_rot = 0) is chosen as the epoch of maximum optical brightness in the rotational cycle.

Superimposed onto spot modulation, additional signs of intrinsic variability are visible. This is the case, for instance, of an apparent brightness event that is centered on JD 2 459 509 and lasted for a few days, as well as more pronounced wide dips toward the end of the observing period, from JD 2 459 544 on. Interestingly, the first half of the light curve exhibits several day-long brightening events, while the last third is dominated by wide dimming events. We note that most of the brightening events, with g′≤12.5 mag, tend to occur at rotational phases shorter than 0.15 or longer than 0.85, that is, close to the maximum brightness of the spot modulation. If the photometric variations of the system are modulated by the visibility of a bright accretion spot at the stellar surface, these brightening events seen around the time of maximum accretion shock visibility most likely reflect varying accretion on the stellar surface. The wide dips in the last part of the light curve exhibit the same periodicity and phase as the spot modulation. They reach a minimum brightness close to phase 0.5, with an amplitude that steadily decreases from 0.6 mag to 0.3 mag in the g′ band over a timescale of a few weeks, from JD 2 459 550 to JD 2 459 570.

Figure 4 shows the color behavior of the system. As the system fades, it becomes redder, with a color slope larger than expected from ISM-like extinction at least in the (u′−g′) and (g′−r′) color indices. According to Venuti et al. (2015), a large color slope at short wavelengths is characteristics of accretion-driven photometric variability. However, both the large scatter seen in the (u′−g′) color index at low brightness levels and the changing shape of the light curve during the semester suggest that several sources of variability might be present, such as a combination of accretion and obscuration events.

Fig. 4. Color-magnitude relation. The (u′−g′), (g′−r′), and (r′−i′) colors are plotted as a function of the system brightness in the g′-band. The dashed lines indicate the expected ISM reddening slope computed from Fiorucci & Munari (2003) with R = 3.1.

3.2. High-resolution optical spectroscopy

We took advantage of the wide wavelength range covered by the SOPHIE spectrograph at high spectral resolution to derive the stellar parameters of GM Aur, to measure optical veiling, and to investigate the emission line profiles and their variability across the optical range. These results are described in the next subsections.

3.2.1. Stellar parameters

To derive the stellar parameters (T_eff, v sin i, V_r, and V_mic), we fit synthetic ZEEMAN spectra (Landstreet 1988; Wade et al. 2001; Folsom et al. 2012) to the average of the first six SOPHIE spectra, which are not contaminated by the moon. Synthetic spectra were computed from MARCS atmosphere models (Gustafsson et al. 2008) and the VALD line list database (Ryabchikova et al. 2015). We applied a χ² minimization procedure based on a Levenberg-Marquart algorithm over seven wavelength windows ranging from 455 to 649 nm, excluding the regions affected by tellurics, emission, or molecular lines. We set the macroturbulent velocity to 2.0 km s⁻¹ and the surface gravity log g to 4.0, and we assumed solar metallicity. These are typical parameters for low-mass TTSs. During the fitting procedure, we applied the mean veiling value derived below for the GM Aur mean optical spectrum (r_0.55 = 0.3, see Sect. 3.2.2) to the synthetic spectra. Finally, we averaged the results obtained from the various wavelength windows, except for one window that yielded values higher than 2σ from the mean. We thus derived T_eff = 4287 ± 35 K, V_r = 14.94 ± 0.14 km s⁻¹, V_mic = 3.3 ± 0.4 km s⁻¹, and v sin i = 14.9 ± 0.3 km s⁻¹. We also used the ROTFIT package (Frasca et al. 2003, 2006) applied to the mean GM Aur’s SOPHIE spectrum, from which we derive T_eff = 4505 ± 53 K, V_r = 15.37 ± 0.26 km s⁻¹, and v sin i = 13.0 ± 0.7 km s⁻¹.

While all these values are consistent within 3σ, the large T_eff difference derived from ZEEMAN and ROTFIT illustrates model-dependent uncertainties that are likely related to the use of different model templates (MARCs for ZEEMAN vs. BTSetll for ROTFIT) and possibly to wavelength-dependent systematics induced by starspots (Gangi et al. 2022; Flores et al. 2022). For consistency with similar observing campaigns that we previously performed on young stars (e.g., Pouilly et al. 2020, 2021; Bouvier et al. 2020a; Alencar et al. 2018), we adopted the results of the ZEEMAN fitting. This estimate also agrees better with the K6 spectral type derived by Herczeg & Hillenbrand (2014), from which they deduced T_eff = 4115 K. When the T_eff-SpT conversion tables from Herczeg & Hillenbrand (2014) and Pecaut & Mamajek (2013) are used, the T_eff value derived above corresponds to a K4-K5 spectral type, which agrees fairly well with previous estimates obtained from optical and near-infrared spectroscopy (K5-K6; e.g., Espaillat et al. 2010; Herczeg & Hillenbrand 2014).

The V_r and v sin i values can be compared to those derived from the uncontaminated CCF of the first six SOPHIE spectra, namely ⟨V_r⟩ = 14.79 ± 0.40 km s⁻¹ and ⟨v sin i⟩ = 12.25 ± 0.26 km s⁻¹. The quoted uncertainties are the rms of the six measurements and therefore include intrinsic variability of the CCF profiles due to spot modulation, for example. The two estimates of V_r agree within the errors as well as with the median V_r value derived from the SPIRou spectra (see Sect. 2.4), while the v sin i value derived from the CCF is significantly lower than the value deduced from spectral fitting. We suspect that the discrepancy may arise from the color-dependent FWHM-v sin i relation used for SOPHIE CCFs, which is calibrated on main-sequence stars (Boisse et al. 2010) and may not be fully adequate for pre-main-sequence objects.

From the average g′r′i′ colors reported for GM Aur in Sect. 2.1, the intrinsic colors of a K4-K5 dwarf listed by Kraus & Hillenbrand (2007) and Covey et al. (2007), and using ISM extinction coefficients from Fiorucci & Munari (2003), we derive a visual extinction on the line of sight A_V = 0.3 ± 0.3 mag. This agrees with previous determinations (e.g., Herczeg & Hillenbrand 2014). The REM photometry reported above yields a median J-band magnitude of 9.42 ± 0.11 mag, which is close to that of 2MASS (J = 9.34 mag). We dereddened the median J-band magnitude with A_j = 0.28 × A_V = 0.084 ± 0.084 mag and used the J-band bolometric correction listed in Pecaut & Mamajek (2013) for a K4-K5 dwarf, BC_J = 1.55 ± 0.03 mag, to derive the stellar luminosity, L_⋆ = 0.9 ± 0.2 L_⊙, assuming the Gaia distance of 157.9 ± 1.2 pc (Gaia Collaboration 2021). We thus derive a stellar radius R_⋆ = 1.7 ± 0.2 R_⊙ from Stefan’s law, and a stellar mass M_⋆ = 1.05 ± 0.05 M_⊙ from the Siess et al. (2000) pre-main-sequence evolution models, while CESAM models (Marques et al. 2013) yield M_⋆ = 0.88 ± 0.12 M_⊙ (E. Alécian, priv. comm.). We therefore adopt M_⋆ = 0.95 ± 0.13 M_⊙, in agreement with Baraffe et al. (2015) models. This estimate is also consistent within the errors with the dynamical mass estimate, M_dyn = 1.00 ± 0.02 M_⊙, reported by Guilloteau et al. (2014), which was later revised to M_dyn = 1.14 ± 0.02 M_⊙ by Simon et al. (2019). Table 4 summarizes the derived stellar parameters.

Finally, we combined the v sin i and rotational period measurements with the stellar radius estimate to derive the stellar inclination sin i = P_rot × v sin i/(2 πR_⋆) = 1.05 with the values listed in Table 4. Accounting for 1σ uncertainties on the stellar parameters, we derive a lower limit of i_⋆≥ 63° for the stellar rotational axis onto the line of sight. This value is significantly higher than the inclination inferred from high-resolution ALMA images of the outer disk of GM Aur observed at millimeter wavelength, which yield i_disk = 53.2 deg (Huang et al. 2020). It agrees better, however, with the inclination value derived for the disk seen in scattered light with adaptive optics on a scale of a few dozen au, for which Oh et al. (2016) obtained i_disk = 64 ± 2 deg. On the much smaller scale of 0.013 au, Bohn et al. (2022) measured an inner-disk inclination i_disk = 68${}_{-28}^{+16}$ deg from long-baseline K-band interferometry using VLTI/GRAVITY. They did not find evidence for an inner and outer disk misalignment in this system. As outlined by Appenzeller & Bertout (2013), the uncertainty on the determination of stellar inclinations from rotation measurements rapidly increases at large angles and is prone to systematic errors. For GM Aur, inferring the stellar inclination from the disk inclination might therefore be more reliable than estimating it from rotation measurements, owing in particular to the significant uncertainty on the stellar radius. Nevertheless, all independent measurements indicate a moderate to high inclination for the system.

3.2.2. Veiling

At optical wavelengths, accreting T Tauri stars exhibit an additional source of continuum flux, which presumably arises from the accretion shock at the stellar surface. This optical excess continuum partly veils the stellar photospheric spectrum and is therefore referred to as “veiling” (Hartigan et al. 1995). We measured the optical veiling from the high-resolution OHP/SOPHIE spectra by comparing the photospheric spectrum of GM Aur to that of the unveiled nonaccreting template V819 Tau, a WTTS with T_eff = 4250 ± 50 K, V_r = 16.6 km s⁻¹, and v sin i = 9.5 km s⁻¹ (Donati et al. 2015). We retrieved an archival CFHT/ESPaDOnS spectrum of V819 Tau, which we resampled at the spectral resolution of OHP/SOPHIE spectra, translated into the radial velocity of GM Aur, and rotationally broadened using the rotational function from Gray (1973) to match the rotational velocity of GM Aur. The template spectrum was then fit to the GM Aur spectra over the same wavelength windows as discussed in Sect. 3.2.1 by adjusting the veiling using the following formula:

$$\begin{array}{c}\hfill I=\frac{I_0+r}{1+r},\end{array}$$(2)

where I is the veiled spectrum, I₀ is the spectrum without veiling, and r is the continuum veiling. Rei et al. (2018) showed that an additional line-dependent veiling component may be present in the strongest photospheric lines. We therefore retained only weak to moderate lines with an EW between 0.01 and 0.1 Å to perform the fit. The veiling measured for each spectrum in each spectral window is shown in Fig. 5. Systematic offsets are clearly seen between the veiling values measured in the different spectral windows. These offsets may partly reflect a wavelength-dependent veiling, as veiling seems to decrease toward longer wavelengths for all but the bluest spectral window, but it may also result from systematic errors depending on the specific sample of lines included in each window. We therefore computed an average veiling value, r_0.55, over all spectral windows for individual GM Aur spectra. This value is listed in Table 5 with its associated uncertainty. The optical veiling is moderate, ranging from 0.17 to 0.55 at 5500 Å. Similar mean values were derived from ROTFIT, namely r = 0.5, 0.4, and 0.3 at 4500, 6000, and 6500 Å, respectively, using a library of 400 spectral templates with spectral types FGKM from the OHP/ELODIE database, which yielded a best χ² fit for spectral types ranging from K3.5 to K5. Figure 6 shows the mean optical veiling plotted along the rotational phase. A hint of higher veiling values towards phases 0 and 1, and minimum values around phase 0.5–0.6 appears. A periodogram analysis of the mean veiling variation did not yield significant results, however.

Fig. 5. Optical veiling measurements from the OHP/SOPHIE spectra. Top: part of a spectral window showing the template spectrum (gray), the observed spectrum (black), and the velocity broadened and veiled template spectrum (red) that fits the observed spectrum. Bottom: veiling measured in several spectral windows (see text) as a function of Julian date. The central wavelength of the spectral windows is indicated in the top left corner of the panel.

Fig. 6. Mean optical veiling plotted as a function of rotational phase. The color code indicates the Julian date.

The lack of significant veiling variability is confirmed by the examination of the LiI 6707 Å photospheric line profile shown in Fig. 7. The shape and depth of the line appear to remain quite stable over the duration of the observing run, and a periodogram analysis across the line profile (Giampapa et al. 1993) revealed no signs of periodic modulation. This suggests that the source of optical veiling remains at least partly in view throughout the rotational cycle. This might arise from a high-latitude accretion shock.

Fig. 7. LiI 6707 Å line profile. Left: the 34 line profile measurements from SOPHIE spectra are shown superimposed. The color code corresponds to successive rotational cycles. Right: 2D periodogram across the line profile. The dotted horizontal line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve displays the mean line profile. The color code reflects the periodogram power.

3.2.3. Emission lines

The main emission lines seen in the optical spectrum of GM Aur, namely Hα, Hβ, and HeI 5876 Å, are displayed in Fig. 8⁸. The Balmer lines show a broad emission peak and a slightly blueshifted absorption component, whose depth varies from being hardly discernable to reaching below the continuum in the Hβ profile. The wide, slightly blueshifted absorption component peaks at −30 to −20 km s⁻¹ and covers a velocity range from about −90 to +40 km s⁻¹ in the Hβ line profile. Additional absorption components appear at higher blueshifted velocities, peaking from −110 to −90 km s⁻¹, and extending down to about −160 km s⁻¹. These blueshifted absorption components cause most of the variability in the Balmer line profile. The Hβ profile also exhibits significant variability over the red wing, up to velocities of +200 km s⁻¹. However, none of the profiles exhibits redshifted absorption components reaching below the continuum level.

Fig. 8. Optical line profile variability. Top: series of optical line profiles Hα, Hβ, and HeI plotted as a function of Julian date (left subpanels) and rotational phase (right subpanels). The color code corresponds to successive rotational cycles. Middle: same profiles, superimposed to illustrate their variability. Bottom: 2D periodograms across the line profiles. The color code reflects the periodogram power, from zero (blue) to 0.6 (red). The dotted horizontal red line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve is the mean line profile.

Owing to the complex line shapes, equivalent widths (EW) of the Hα, Hβ, and HeI 5876 Å lines were computed by directly integrating below the line profile. The results are listed in Table 5. The measurement accuracy is 10% or better for EW(Hα) and EW(Hβ), and about 20% for EW(HeI) due to the more uncertain continuum location. We note that on JD 2 459 512, the EWs measurements are systematically lower than during the rest of the observations, which might be due to lunar contamination, as discussed in Sect. 2.3. The average and rms values we obtain are EW(Hα) = 83 ± 8 Å, EW(Hβ) = 9.3 ± 2.5 Å, and EW(HeI) = 0.40 ± 0.13 Å.

The HeI 5876 Å line profile is roughly symmetric and consists of a narrow component (FWHM ∼ 30–40 km s⁻¹) superimposed on a broad pedestal, as previously reported by McGinnis et al. (2020). The peak intensity of the narrow component varies significantly, while the broad component appears relatively stable (see Fig. 8). We fit the HeI line profile with a two-component Gaussian model to extract the properties of the narrow (NC) and broad (BC) components. The EWs were derived from the Gaussian fit of the components. The radial velocity, FWHM, and EW of the NC and BC, as well as their uncertainty derived from the covariance matrix of the Levenberg-Marquart fitting algorithm, are listed in Table 6.

We derived the radial velocity of the narrow component and found it to be variable and redshifted by ∼5–10 km s⁻¹ relative to the stellar velocity. Figure 9 shows the HeI NC radial velocity curve plotted as a function of Julian date and rotational phase. As previously reported by McGinnis et al. (2020), HeI NC V_r appears to be rotationally modulated with a full amplitude of ∼6 km s⁻¹, as expected if the NC component of the HeI 5876 Å line profile were produced in a high-latitude accretion shock at the stellar surface. We fit the observed NC V_r curve with the geometrical accretion shock model described in Pouilly et al. (2021). The model computes the variation of the HeI NC radial velocity as the combination of the rotational modulation of the accretion shock and the intrinsic inflow velocity. The free parameters of the model are the inflow velocity, the latitude of the accretion shock, the phase at which it faces the observer, and the stellar inclination. The HeI NC V_r curve is best reproduced with an accretion shock located at a latitude of 83 ± 1.5° that faces the observer at phase 0.2 ± 0.08 and has a radial post-shock velocity of 18.3 ± 1.0 km s⁻¹ in the stellar rest frame. Because the model now includes the inflow velocity, the HeI NC radial velocity curve does not reach the mean velocity at the time when the spot faces the observer, as it would for the case of static stellar spots. The stellar inclination we derive from the model is i = 64 ± 2.2°, which is consistent with the inner disk inclination derived from K-band VLTI/GRAVITY data (see Sect. 3.2.1). According to the model, the accretion shock faces the observer close to the origin of the phase, which is consistent with the photometric behavior described above (see Sect. 3.1). The HeI line profile is also strongest over rotational phases ranging from 0.75 to 1.09 (excluding the probable accretion burst occuring at JD 2 459 510, see Sect. 3.1), and this is also when the highest veiling values are measured in the JHK bands (see Fig. D.1). These two results further support a maximum visibility of the accretion shock close to Φ_rot = 0.

Fig. 9. Radial velocity of the narrow component of the HeI 5876 Å line profile in the stellar rest frame plotted as a function of Julian date (top) and rotational phase computed from the ephemeris of Eq. (1), (bottom). The color code corresponds to successive rotational cycles. The fit by a geometrical accretion shock model (see text) is shown (dash-dotted curve) together with its 1σ uncertainty (gray area).

A bidimensional periodogram analysis (Giampapa et al. 1993) of the Balmer and HeI 5876 Å line profiles reveals a periodic modulation of part of the profiles (see Fig. 8). The intensity of the narrow component of the HeI line profile is modulated at a frequency of ∼0.15 d⁻¹, corresponding to a period of 6.7 d, with a large uncertainty due to the limited temporal sampling of the spectral series, however, that translates into a poor frequency resolution. As the HeI line profile NC component arises in the accretion shock (Beristain et al. 2001), it is expected to be modulated at the stellar rotational period or close to it in case of latitudinal differential rotation at the stellar surface. In the Balmer line profiles, a modulation close to the stellar period also appears over three distinct locations: in the highly redshifted wing over the velocity range ∼200–400 km s⁻¹, at slightly redshifted velocities between 0 and ∼50 km s⁻¹, and at highly blueshifted velocities from −400 to −200 km s⁻¹. While the maximum power of the periodogram in the blue wing appears at the stellar rotation period, it seems to drift to longer periods, similar to that seen in the HeI NC component, toward the red wing. If the red wing modulation of Balmer line profiles is caused by the absorption of shock emission by a funnel flow crossing the line of sight, this might be an indication of differential rotation along the funnel flow. Noticeably, no sign of periodic variability is seen in the Balmer line profiles in the velocity channels in which the variable blueshifted absorption components arise. This suggests that these components either result from sporadic ejection processes or that they vary on periods longer than ten days.

Correlation matrices (Johns & Basri 1995) between line profiles were computed and are presented in Appendix B. They display the degree to which temporal flux variations in a pair of spectral lines are correlated. Matrices can be computed for a single profile (autocorrelation), for instance, Hα ⋆ Hα, to investigate how the different parts of the profile vary with respect to each other, or between two profiles (cross-correlation), for example, Hα ⋆ Hβ, to compare the intensity variations of different lines. The Hα ⋆ Hα and Hβ ⋆ Hβ matrices shown in Fig. B.1 are quite similar. The blue and red wings of the profiles vary in a correlated way. The high-velocity red wings are anticorrelated with the emission peak region. This may be a sign that high-velocity redshifted absorption components appear when the peak intensity is higher, although the absorptions do not reach below the continuum. Many absorption features are superimposed on the emission line component. These features do not present a periodicity and are not correlated with each other, nor with the emission part of the profile. The Hβ ⋆ Hα matrix mimics the matrices of Hα ⋆ Hα and the Hβ ⋆ Hβ, showing that both lines are formed in the same region, as expected from magnetospheric accretion models (e.g., Muzerolle et al. 2001).

3.3. High-resolution near-infrared spectroscopy

The main emission lines seen in the high-resolution near-infrared spectrum of GM Aur, namely HeI 10830 Å, Paβ, and Brγ, are depicted in Figs. 10 and 11⁹. The near-infrared hydrogen lines exhibit broad emission profiles, with FWHM ∼ 200 km s⁻¹, whose peaks are slightly blueshifted compared to the stellar velocity, and are located at −30 and −10 km s⁻¹ for Paβ and Brγ, respectively. The profiles are roughly symmetric, but have pronounced high-velocity redshifted absorptions that extend up to +400 km s⁻¹ and reach below the continuum level. Paβ and Brγ exhibit inverse P Cygni (IPC) profiles of varying depth in most SPIRou spectra, which is in stark contrast to the optical hydrogen line profiles, Hα and Hβ, which do not exhibit such features (see Sect. 3.2.3).

Fig. 10. Near-infrared line profiles: HeI (left), Paβ (center), and Brγ (right), plotted with arbitrary offsets as a function of Julian date. Each color represents a SPIRou run, namely September (blue), October (orange), November (green), December 2021 (red), and January 2022 (purple). The October SPIRou run is contemporaneous to the OHP/SOPHIE observations.

Fig. 11. Near-infrared line profile variability. Top: series of residual near-infrared line profiles HeI (left), Paβ (center), and Brγ (right), plotted superimposed to illustrate their variability. Each color represents a SPIRou run as in Fig. 10. Bottom: 2D periodograms across the line profiles. The color code reflects the periodogram power from zero (blue) to 0.5 (red). The dotted red horizontal line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve displays the mean line profile.

Correlation matrices for the near-infrared hydrogen line profiles are shown in Fig. B.2. The Paβ ⋆ Paβ, Brγ ⋆ Brγ, and Paβ ⋆ Brγ matrices are quite similar. The emission part of the profiles is overall correlated, and anti-correlated with the high-velocity redshifted absorption component. This indicates that the redshifted absorption is deeper when the emission line is more intense. Cross-correlation matrices between optical and near-infrared line profiles observed in the same night during the October runs are shown in Fig. B.3. The Brγ and the Paβ emission components correlate well with the main Hα and Hβ emission components (−100 km s⁻¹ < v < 200 km s⁻¹), which suggests that at least part of the line emission forms in the same extended region. The Paβ redshifted absorption, and to a lesser extent, the Brγ redshifted absorption, correlates well with the Hα and Hβ blue and red high-velocity wings. This indicates that they all form in the same region close to the star. The near-infrared profiles are strongly anticorrelated with the highly variable blueshifted absorption components that appear at velocities of about −100 km s⁻¹ in the Balmer line profiles and presumably trace outflows. Except for this feature, the overall correlated variations of optical and near-infrared line profiles suggest that they are all good diagnostics of the accretion flow, as previously suggested by Alcalá et al. (2014), for example. Detailed modeling of the line profile shapes and variability is needed, however, to ascertain the exact location from which they arise (Tessore et al. 2023).

The HeI 10830 Å line profile is more complex. It is usually dominated by an emission peak at low velocities that becomes quite weak at specific times, however. The profile also often displays highly redshifted absorption features, similar to the IPC components seen in the Paβ and Brγ lines. Unlike the near-infrared hydrogen profiles however, the HeI line additionally exhibits various absorption components in the blue wing of the profile. At least two systems of blueshifted absorptions can be identified: a low-velocity system extending between −20 to −80 km s⁻¹, and a high-velocity system ranging from −100 to −250 km s⁻¹. Figures 10 and 11 show that the two systems are quite variable. These two absorption systems are reminiscent of those observed in the Hα and Hβ lines, although they occur at slightly bluer velocities in the HeI profile, with an offset of about −40 km s⁻¹ compared to the optical profiles.

Correlation matrices involving the HeI 10830 Å line are displayed in Fig. B.4. The HeI autocorrelation matrix shows several correlated components. Over the velocity channels extending from −100 km s⁻¹ up to 200 km s⁻¹, the line flux varies in a correlated fashion, but this region of the line does not correlate with the rest of the profile. Similarly, the redshifted absorption region, which extends from 200 km s⁻¹ to 400 km s⁻¹, varies as a whole and does not correlate with the rest of the profile. The blue wing of the profile presents many short-duration absorption components superimposed on the emission component, and, probably due to this, each velocity bin varies independently of the other. Comparison of the HeI 10830 Å line profile to optical and near-infrared hydrogen profiles indicates that the emission core and peak intensity are correlated, as are the high-velocity redshifted absorptions seen in the HeI, Paβ, and Brγ line profiles. However, the HeI blue wing shows a strong anticorrelation with the core of the Hβ profile.

The periodogram analysis of the line profiles is shown in Fig. 11. In the three lines, the IPC components appear to be periodically modulated at the stellar rotation period. As the SPIRou dataset extends over four months, this suggests that a quite stable structure, presumably the accretion funnel flow, gives rise to this spectral feature¹⁰. The redshifted absorption component is deeper when the emission line is more intense, which agrees with the assumption that the two components trace the densest part of the accretion column close to the star. The low-velocity red wing of the HeI line profile also shows a periodic modulation at the same period, which might trace the visibility of the accretion shock, as seen in the optical HeI line profile. The periodogram power in the blue wing of the lines is weaker, except perhaps over the high-velocity channels of the Paβ line. In particular, the variable blueshifted absorption systems seen in the HeI line profile are not periodically modulated on this long timescale. We performed a similar analysis of each SPIRou run individually, and the results are shown in Appendix E. Figures E.1–E.4 reveal that the most conspicuous high-velocity redshifted absorptions in the HeI, Paβ, and Brγ line profiles occur preferentially around Φ_rot = 0. During the SPIRou October run, which was simultaneous with the OHP/SOPHIE observations, we realized that the highly blueshifted velocity channels of the near-infrared lines appear to be modulated at the stellar rotation period (see Fig. E.2). Hence, while the modulation of high-velocity blueshifted absorptions seen in the HeI profile disappear on a timescale of several months, the modulation may survive over a few rotational periods. Finally, similar to what is seen in the Balmer line profiles over a much shorter time span (see Sect. 3.2.3), there is a marginal indication from the periodograms of the near-infrared lines that the modulation period drifts from the blue to the red wing of the profiles. This might be a sign of differential rotation in the source of the variability.

The EW of the HeI 10830 Å, Paβ, and Brγ lines was computed on the residual profiles using two methods. The first method consists of adjusting a Gaussian fit to the line profile, and the second method is integrating below the line profile. The two methods were applied to the Paβ and Brγ lines that most often exhibit a strong Gaussian-like emission and, at times, a pronounced redshifted absorption. The difference between the EW derived from the Gaussian fit and from the profile integration then provides an estimate of the strength of the redshifted absorption component. The HeI line exhibits a complex profile, with pronounced absorptions appearing in both the blue and red wings, and it cannot be fit by a Gaussian. We report here HeI EWs measured through profile integration only. EW measurements are usually accurate to within 5% because of the well-defined adjacent stellar continuum level.

The results are listed in Table 7, and the evolution of EW during the observing campaign is illustrated on Fig. 12. While the EW variations of the three lines are usually correlated, there are notable exceptions. For instance, during the first SPIRou run around JD 2 459 478, the HeI line goes into absorption, driven by a broad redshifted absorption component that reaches half the continuum value, while the Paβ and Brγ lines reach a local intensity maximum, even though they also exhibit a pronounced redshifted absorption (see Fig. E.1). In contrast, during the second SPIRou run, the variations of the three lines are well correlated. Line EWs and near-infrared veiling measurements (listed in Table 3) are both shown as a function of Julian date and rotational phase in Fig. D.1. Both quantities appear to be rotationally modulated, with maximum values occurring close to Φ_rot = 0.

Fig. 12. EW of near-infrared lines. Top: EW of the HeI, Paβ, and Brγ lines measured on SPIRou spectra plotted as a function of Julian date. The measurement uncertainties are about the symbol size. Bottom: difference between the EWs measured by direct line profile integration (i.e., including the redshifted absorption component) and the EW derived from the Gaussian fitting of the emission component only plotted as a function of Julian date. This differential quantity measures the strength of the redshifted absorption component in the Paβ and Brγ line profiles. The more negative the differential quantity, the deeper the redshifted absorption.

Finally, the photospheric radial velocity curve derived from SPIRou spectra is shown in Fig. 13. V_r is found to vary between 13.96 and 15.04 km s⁻¹, with a median value of 14.65 km s⁻¹. A CLEAN periodogram analysis reveals a period P = 5.94 ± 0.11 d, and the String-Length method yields P = 5.98 ± 0.12 d, both consistent with the 6.04 ± 0.15 d photometric period derived in Sect. 3.1. The photospheric radial velocity curve folded in phase at the stellar rotation period is shown in Fig. 13. V_r exhibits a roughly sinusoidal variations in rotational phase, which suggests that it is modulated by a surface spot. It reaches the median velocity going blueward around Φ_rot = 0, as expected from a stellar spot facing the observer at this phase. Because this is also the phase of maximum brightness of the system (see Sect. 3.1), this suggests that a hot spot modulates the photospheric V_r curve. The V_r curve derived from the near-infrared photospheric lines (see Fig. 13) is inverted compared to the V_r curve derived for the HeI 5876 Å NC line profile (see Fig. 9). Similar antiphase radial velocity variations between absorption and emission lines have been reported for other accreting T Tauri stars and were interpreted as caused by the modulation of the line profiles by a hot spot at the stellar surface (Petrov et al. 2001, 2011; Gahm et al. 2013). Alternatively, the photospheric radial velocity variations may also be produced by a cool spot that coexists with the accretion shock at the same location on the stellar surface, as Doppler images of accreting T Tauri stars suggest (e.g., Donati et al. 2010, 2011, 2019, 2020a).

Fig. 13. Radial velocity variations measured in the near-infrared photospheric lines of the SPIRou spectra. Top: radial velocity as a function of Julian date. Bottom: radial velocity curve folded in phase at the stellar rotational period P = 6.04 days. The color code reflects the Julian date.

3.4. Low-resolution near-infrared spectroscopy

The median nightly spectra obtained from the ExTrA-T2 and ExTrA-T3 telescopes are shown in Fig. 14. We computed the EW of the HeI 10830 Å, Paβ, Paγ, and Paδ lines from the ExTrA spectra using specutils¹¹. We analyzed each telescope independently because the point spread function (and therefore the resolution) of the spectrograph depends on the position on the detector. First, we fit a Gaussian to each line on the median spectra to measure the line center and full width, the latter we defined as amounting to six times the standard deviation of the Gaussian fit in order to isolate the line profile from the nearby continuum. The local continuum around each line was modeled with a third-degree polynomial, adjusted on three line widths centered on the line, but excluding the line. Then, for each line in each of the 1898 individual spectra, we computed the EW by integrating the flux over the full line width using the parameters and the normalization region defined from the median spectra. Finally, we computed the median of the individual measurements for each night, regardless of the telescope. The median is less affected by outliers than the mean, and the differences between the mean and the median values are within the error bars. Table C.1 lists the results.

Fig. 14. Median spectrum for each night from the ExTrA-T2 telescope (lower part, 88 spectra) and ExTrA-T3 telescope (upper part, 69 spectra). The main emission lines are indicated. The color is proportional to the EW of the Paβ line shown in the top panel as a function of Julian date.

In order to estimate the reliability of the procedure, we compared EWs derived from ExTrA and from SPIRou spectra for the Paβ line using 20 measurements obtained on each instrument less than one day apart. The results show an excellent correlation, with a slight tendency for the EW measured from ExTrA to exceed those measured from SPIRou, with a mean difference of 0.7 ± 1.1 Å. A similar result is obtained for the HeI line, with a mean difference of 0.31 Å and an rms of 1.1 Å between ExTrA and SPIRou estimates. The comparison with high-resolution measurements thus validates the EWs obtained from low-resolution spectra.

The line variability is illustrated on Fig. 14, where the median nightly spectra are superimposed. The median and extreme values of EWs measured for the HeI, Paβ, Paγ, and Paδ lines are listed in Table 8, and the night-by-night measurements are listed in Table C.1. We note that the HeI line at times appears to be in absorption at this low spectral resolution.

Figure 15 shows the EW measurements plotted as a function of time and its generalized Lomb-Scargle periodogram (GLS; Zechmeister & Kürster 2009). The EW of the four lines is found to be modulated with a period of 6.028 ± 0.087 days, consistent with the stellar rotation period, where the error estimate is the standard deviation of a Gaussian fit to the periodogram peak. Figure 16 shows the HeI and Paβ line EWs folded in phase at the stellar rotation period. The modulation of the line strength clearly appears, with maximum flux around phase zero. As shown in the same figure, these variations follow the modulated brightness level of the system in the J band.

Fig. 15. EW variability of the near-infrared line profiles. Left: EW measurements plotted as a function of Julian date for the HeI, Paβ, Paγ, and Paδ lines from the ExTrA spectra. Right: GLS periodogram of the EW measurements. The period on the x-axis is displayed on a log scale.

Fig. 16. EWs of near-infrared line profiles and J-band magnitude variability. Top: EW of the near-infrared HeI and Paβ lines folded in phase at the stellar rotational period. EW(HeI) is offset by −10 Å for clarity. Bottom:J-band light curve, deduced from ExTrA spectra, folded in phase at the stellar period. The brightness modulation is similar to that observed at optical wavelengths (see Fig. 2).

Longer-term EW variations of higher amplitudes are also clearly seen in Fig. 15. These variations seem to be correlated with the multicolor photometry presented in Sect. 3.1, in particular, during the brightness event centered on JD 2 459 509 and the wide dip around JD 2 459 549. In Fig. 17, a clear correlation appears between the photometric variations in the u′ band and the Paβ line variations. If the brightening of the system in the u′ band is linked to accretion, this correlation suggests that most of the emission line flux is connected to the same process.

Fig. 17. Mid-term variation of EW(Paβ), (red) compared to the system u′-band light curve (blue) for measurements taken less than one day apart. To facilitate the comparison, the EW measurements are plotted on a magnitude scale and are offset, namely −2.5 log EW(Paβ) + 5.

Fig. 18. Magnetospheric ejection model. Development of a magnetospheric ejection computed from star-disk interaction MHD simulations. The snapshots shown here are extracted from model 3 of Pantolmos et al. (2020), where the magnetospheric truncation radius amounts to 54% of the corotation radius. The three snapshots are shown at 0.2, 0.5, 0.8, and 1.1Prot. The white curves indicate expanding magnetic field lines that give rise to the ejection of plasmoids. The color scale indicates density. The green arrows show the velocity field of the stellar and disk winds and of MEs at their interface.

3.5. Mass-accretion rate

The combination of optical and near-infrared spectroscopy offers a number of emission lines from which we can estimate the mass accretion rate onto the star, using the empirical relations between line luminosity and accretion luminosity proposed by Alcalá et al. (2017). Combining the range of Hα and Hβ EWs reported above with the nearby continuum fluxes computed from the r′ and g′ magnitudes corrected for extinction, we obtain the line fluxes and luminosities as follows: $F_{\text{line}}=F_{\lambda }^o\times \text{EW}(\text{line})\times {10}^{-0.4(m_{\lambda }-A_{\lambda })}$ and L_line = 4πd²F_line, where $F_{\lambda }^o$ is the flux of a zero-magnitude star¹², m_λ and A_λ are the magnitude and extinction in the photometric band of interest, and d is the distance to the star. The accretion luminosity was derived from the line luminosity using the relation reported by Alcalá et al. (2017), and Ṁ_acc was deduced from L_acc assuming a magnetospheric radius of 5R_⋆ (see Alcalá et al. 2017, Eq.(1)). Taking the uncertainties on all involved quantities into account, we thus derive Ṁ_acc = 0.7 ± 0.3 and 0.5 ± 0.4 × 10⁻⁸ M_⊙ yr⁻¹ from the Hα and Hβ line fluxes, respectively.

We performed a similar analysis using the extensive measurements of Paβ EWs obtained from the ExTrA spectra. From the extreme values, namely EW(Paβ) = 2.8–18.6 Å, and the mean REM J-band magnitude of the system during the observing period, J = 9.41  ±  0.10, we derive Ṁ_acc = 0.3–2.0 × 10⁻⁸ M_⊙ yr⁻¹, with a median value of Ṁ_acc = 1.0 × 10⁻⁸ M_⊙ yr⁻¹. Finally, we derived additional estimates of Ṁ_acc from the Brγ line flux, using the median and extreme values of EW(Brγ) measured on SPIRou spectra, namely 5.1, 1.7, and 9.7 Å, which we calibrated with the mean REM K′-band magnitude of the system, K′ = 8.40. Using the relations reported by Alcalá et al. (2017) between line and accretion luminosity, we thus derived Ṁ_acc = 0.2–1.8 × 10⁻⁸ M_⊙ yr⁻¹, with a median value of Ṁ_acc = 0.8 × 10⁻⁸ M_⊙ yr⁻¹.

The dispersion in the Ṁ_acc estimates partly results from intrinsic Ṁ_acc variability. For instance, the SOPHIE spectra were obtained during a brightening of the system that occurred around JD 2 459 512, at a time at which all optical and infrared lines were relatively strong. We thus derive a relatively high value of Ṁ_acc from these spectra compared, for example, to the smaller Ṁ_acc estimate derived from the Brγ line in the SPIRou spectra that were obtained during more quiescent phases of the system. We cannot exclude, however, that some of the dispersion may also arise from systematic uncertainties in the empirical line-to-accretion luminosity relations over the optical and near-infrared wavelength ranges. In any case, the various estimates agree globally, with Ṁ_acc typically varying between 0.2 and 2.0 × 10⁻⁸ M_⊙ yr⁻¹.

4. Discussion

The monitoring campaign we performed on GM Aur reveals significant but relatively low-level temporal variability over a timescale of six months. The photometric variations are mild, as might be expected for this moderately accreting young system (Ṁ_acc ∼ 0.8 × 10⁻⁸ M_⊙ yr⁻¹), with amplitudes ranging from 1.5 mag in the u′ band to 0.3 mag in the i′ band, and 0.1 mag in the J band. The brightness of the system varies smoothly because it is modulated by the visibility of surface spots at the stellar rotational period of 6.04 d. Except for the HeI 10830 Å line profile, the spectral appearance of the system does not change drastically on a timescale of months, and the veiling is low and stable, amounting to about 0.3 in the optical range. The strength of the main emission lines (Hα, Hβ, HeI 5876 Å, Paβ, and Brγ) varies by a factor of 2–3 over the course of the semester. In contrast, the HeI 10830 Å line profile exhibits extreme variability and is sometimes barely noticeable in emission. The line profile shapes are strongly variable on a timescale of days, and the development of both blueshifted and redshifted absorption components is superimposed onto a broad emission component. Blueshifted absorption components indicate outflows (e.g., Edwards et al. 2003; Kwan et al. 2007), and redshifted components probe funnel flows (e.g., Edwards et al. 2006; Fischer et al. 2008). Remarkably, the redshifted absorption features that reach below the continuum level, that is, inverse P Cygni profiles (IPC; see Calvet & Hartmann 1992), which are seen in the high-resolution near infrared line profiles (HeI 10830 Å, Paβ, Brγ), are steadily modulated by stellar rotation over an extended observing time span of 3 months. Rotational modulation is also clearly detected in the strength of the near-infrared lines and is continuously observed at low resolution for more than five months. Assuming that most of the line flux arises from the magnetospheric accretion region, as suggested by the periodic appearance of the IPCs and the correlation between line flux and u′-band excess, this indicates a stable large-scale accretion structure on this timescale.

In contrast, high-velocity blueshifted absorption components are neither periodic nor stable on this timescale. While they are ubiquitous in the optical line profiles, most notably Hα and Hβ, and are also quite conspicuous in the HeI 10830 Å line profile, their signature evolves on a timescale of a few days, sometimes drifting in velocity before disappearing altogether. As an example, Fig. E.3 shows the evolution of the blueshifted absorption components in the HeI line profile over the course of the SPIRou November run. On JD 9537, a high-velocity blueshifted feature appears in the profile and drifts toward lower velocities over the next several days, from −200 km s⁻¹ to −120 km s⁻¹. Soon after, on JD 9540, a new high-velocity component appears and follows the same trend. A similar behavior is seen in the high-velocity blueshifted component appearing in the HeI 10830 Å line profile during the SPIRou September run, which drifted from −240 km s⁻¹ to −150 km s⁻¹ on a timescale of five days, from JD 9475 to JD 9480 (see Fig. E.1). This suggests episodic outflows lasting for a few days only. We see no sign of a steady, constant velocity wind in the HeI line profiles over the observing period, nor do we find evidence for a rotational modulation of the outflow signatures. The only stable outflow signature seen in the optical and near-infrared line profiles consists of a narrow, low-velocity blueshifted absorption in the Hα profile that peaks at −20 km s⁻¹ and remains visible over more than two weeks.

Observations thus suggest that we witness a globally stable accretion structure and a succession of short-lived episodic outflows. The contrasting behavior of accretion and outflow diagnostics observed on a timescale of months thus raises the question whether the two processes are physically connected on the scale of a few stellar radii that we probe here. To examine this issue, we investigated the aftermath of the brightening event GM Aur underwent around JD 2 459 509. On this date, the system exhibited a significant brightening at optical and near-infrared wavelengths, most notably in the u′ band (∼1 mag), as well as some of the strongest line fluxes and highest optical and near-infrared veiling values measured during the campaign. A simultaneous TESS light curve of the system recorded the brightening event (see Fig. 2), which started on JD 2 459 508 and ended on JD 2 459 511, and exhibited a flat peak lasting for two days with a 20% flux increase. We therefore interpret this episode as an accretion burst that occurred around the rotational phase Φ_rot = 0, corresponding to the maximum visibility of the accretion shock, and lasted for several days. From the measured continuum level and Balmer line EW during the burst, we derive an increase of a factor of 2 in the accretion rate, reaching ∼2 × 10⁻⁸ M_⊙ yr⁻¹. Inspection of the optical and near-infrared line profiles during and after this event reveals high-velocity blueshifted absorption components that appear in the Balmer and HeI 10830 Å line profiles. On JD 2 459 509 a new blueshifted absorption component appears in the HeI infrared line at a velocity of −280 km s⁻¹ and extends down to −350 km s⁻¹ before it reaches the continuum level again, and it drifts to lower velocities (∼−150 km s⁻¹) over the next few days. The closest optical spectrum was recorded only toward the end of the accretion burst, on JD 2 459 510.6. In this spectrum, the Hα and Hβ profiles feature a new blueshifted absorption component at a velocity of −110 km s⁻¹ that lasts for a few days. The contemporaneous occurrence of an accretion burst rapidly followed by outflow signatures in the line profiles therefore suggests that the accretion and ejection processes are physically connected on small scales.

We propose that the most likely scenario that accounts for these episodic events is magnetospheric ejections, possibly triggered by magnetic reconnections in the accreting magnetosphere. Magnetospheric ejections (Zanni & Ferreira 2013; Sauty et al. 2022), or nonstationary conical winds (Romanova et al. 2009), are caused by the expansion and reconnection of the field lines that connect the star with the disk. The inflation process is the result of the star-disk differential rotation and the consequent build-up of toroidal magnetic field pressure (e.g., Goodson et al. 1997). Quasi-periodic ejections of plasmoids are predicted to occur throughout the magnetospheric inflation cycle on a timescale of about the rotational period (Hayashi et al. 1996). The speed and variability of the outflows likely depend on various parameters, such as the magnetic field strength and topology, the thermal disk pressure, nonideal MHD effects, and the interaction of the magnetospheric-ejection region with the surrounding outflows, that is, stellar and disk winds (Miller & Stone 1997; Romanova et al. 2009; Zanni & Ferreira 2013). Previous monitoring campaigns on young stellar objects have reported evidence for magnetospheric inflation cycles (Bouvier et al. 2003; Alencar et al. 2018). We show a 2.5D MHD simulation of the interaction between an inner accretion disk and a dipolar magnetosphere from Pantolmos et al. (2020) in Fig. 18. The figure illustrates the ejection of plasmoids along expanding field lines at the distance of a few stellar radii from the stellar surface. The timescale for successive ejections is about that of the stellar rotation period. The outflow speed, shown in Fig. 19, reaches more than 200 km s⁻¹ in the early phases of the ejection, then decelerates to about 150 km s⁻¹ on a timescale of days (0.3 × P_rot), and finally vanishes altogether. A direct comparison of the model to observations is not straightforward. The terminal speed we derive from observations is higher than predicted by the simulation, and its temporal evolution on a timescale of a few days might be dominated by projection effects as the system rotates and does not trace the evolution of the plasmoid velocity field. Moreover, the sporadic ejections we observe do not appear to have the quasi-periodic character of the model. Both the wind speed and ejection timescale may, however, depend on numerical effects. In any case, the behavior of the magnetospheric ejection model qualitatively matches the dynamics of the high-velocity blueshifted absorptions seen in the line profiles of GM Aur on a timescale of days to weeks. Moreover, in the scenario of a magnetospheric inflation cycle, magnetic reconnection leads to an accretion burst and simultaneously triggers an ejection episode (e.g., Goodson & Winglee 1999). This is quite reminiscent indeed of what is suggested by the variability of GM Aur.

Fig. 19. Magnetospheric ejection model. Velocity of the gas in the ejected plasmoid at a distance of 10 (magenta), 20 (blue), and 30 (cyan) stellar radii as a function of time in units of the stellar rotational period. Successive ejections of plasmoids are featured. The vertical lines correspond to the snapshots shown in Fig. 18.

We also noted a change in the GM Aur light curve whose first part is dominated by successive low-level brightening events, up to the major accretion burst described above, while it exhibits luminosity dips toward the end of the observations. It is unclear whether the contrasting behavior of the system luminosity observed before and after the main accretion burst is a consequence of the burst itself, perhaps inducing a structural change in the star-disk interaction process. It is conceivable that the rearrangement of the magnetic topology after the inflation or reconnection event that led to the burst has impacted the large-scale geometry of the star-disk interaction region. Either a modest increase in the inner disk scale-height (Nagel 2017) or a reduction of the extent of the truncation radius is prone to trigger a periodic obscuration of the central star by circumstellar dust, that is, a dipper phenomenon (Cody et al. 2014), especially in young systems seen at high inclination (McGinnis et al. 2015; Bodman et al. 2017). The magnetic topology and the mass accretion rate may thus have slightly evolved over the six-month time span of the campaign, as suggested by the long-term variations of the emission line EWs shown in Fig. 15. However, the results we report here clearly indicate that the large-scale geometry of the star-disk interaction was not drastically modified over this timescale. This is evidenced by the phase stability of the modulated light curve, the smooth variability of the emission line profiles around their mean shape, and the strictly periodic appearance of IPC profiles over at least three months. All these accretion diagnostics support a globally stable magnetospheric accretion structure during the campaign.

Finally, it is interesting to compare the line profile shape and variability reported here to those reported by McGinnis et al. (2020), which were obtained in 2011, that is, ten years prior to our observing campaign. The shapes of the Hα and Hβ profiles are quite different in the two studies. In McGinnis et al. (2020), these are pure emission profiles with a triangular shape, without any significant absorption components, and they exhibit little variability over the timescale of a week (see their Figs. S5 and S6). Here, the same profiles appear to be much more structured, with highly variable absorption components on the same timescale. The HeI 5876 Å line profile variability also differs between the two studies, but in the opposite direction. In both studies, the line profile consists of a broad and a narrow component. However, in McGinnis et al. (2020), the intensity of the broad component clearly varies, especially on the blue wing, while we found it to be quite stable here. It seems that the behavior of the system was different between the two epochs. The Hα and Hβ line EWs are indeed systematically higher in McGinnis et al. (2020) and were comparable to the highest values we measured here during the JD 2 459 509 accretion burst. It is therefore likely that GM Aur was in a state of more active accretion during the 2011 observations. This is consistent with the more triangular shape and lack of structure of the Balmer emission line profiles, which are predicted to become more optical thick as the funnel flow density increases (Muzerolle et al. 2001). It is also consistent with the higher level of variability seen in the blue wing of the broad component of the HeI 5876 Å line profile that may betray the existence of a hot accretion-driven wind at times of enhanced accretion (Beristain et al. 2001). Long-term changes in the system behavior driven by a varying mass accretion rate and/or a change in the magnetic topology are therefore likely to occur. Over the six-month span of our observing campaign, the significant variation observed in the EW of the Paβ line profile, which by the end of the campaign reaches similar levels to those measured during the JD 2 459 509 accretion burst (see Fig. 15), suggests that such changes may occur on a timescale of a few weeks to months.

5. Conclusion

By combining optical and near-infrared high-resolution spectroscopic time series, seconded by a long-term monitoring of the photometric variability of the system and low-resolution near-infrared spectrophotometry, we were able to characterize the accretion and ejection process occurring in the young system GM Aur on a timescale ranging from days to months. We report a stable accretion pattern according to which the large-scale magnetic field of the star controls the accretion of gas from the inner disk onto the central star along funnel flows. The appearance of inverse P Cygni profiles that signal the crossing of funnel flows on the line of sight is remarkably periodic at the stellar rotation period of 6.04 days. Similarly, the photometric and line flux variations, both driven by the visibility of the accretion shock located at the foot of the main accretion column, are modulated at the same period. While the amplitude varies, the phase of variability of all these accretion diagnostics remains stable over the 30 rotational periods covered by the campaign. This suggests that the underlying magnetic topology that controls the non-axisymmetric accretion flow, presumably an inclined dipole on the large scale, did not undergo major changes over a timescale of six months. In stark contrast, high-velocity blueshifted absorption components that indicate outflows appear at random times in the emission line profiles. They are not rotationally modulated, and their signatures last for a few days only. We argue that these transient outflows associated with a stable accretion pattern are best accounted for by magnetospheric ejection models, as predicted by MHD simulations. Thus, by probing the dynamics of the star-disk interaction region, these results show that the physical connection between accretion and ejection processes that has long been established on large scales also appears to be valid on the much smaller sub-au scales.

Acknowledgments

We thank the referee for a prompt and detailed report. This study is based on observations obtained at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. The observations at the CFHT 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; based on data collected under the ExTrA project at the ESO La Silla Paranal Observatory. ExTrA is a project of Institut de Planétologie et d’Astrophysique de Grenoble (IPAG/CNRS/UGA), funded by the European Research Council under the ERC Grant Agreement n. 337591-ExTrA. We thank Ágnes Kóspál for providing a reduced TESS light curve of GM Aur. Funding for the TESS mission is provided by NASA’s Science Mission directorate. 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, http://www.spidi-eu.org). We acknowledge funding from the French National Research Agency (ANR) under contract number ANR-18-CE31-0019 (SPlaSH). SHPA acknowledges financial support from CNPq, CAPES and Fapemig. JFD acknowledges funding from the European Research Council (ERC) under the H2020 research & innovation programme (grant agreement no. 740651 NewWorlds). AF acknowledges support by the PRIN-INAF 2019 STRADE (Spectroscopically TRAcing the Disk dispersal Evolution) and by the Large Grant INAF YODA (YSOs Outflow, Disks and Accretion). JFG was supported by fundação para a Ciência e Tecnologia (FCT) through the research grants UIDB/04434/2020 and UIDP/04434/2020. This work benefited from discussions with the ODYSSEUS (HST AR-16129) and PENELLOPE teams. Some of the plots presented in this paper were built using TOPCAT (Taylor 2005).

References

Appendix A: REM comparison stars and optical photometry

Table A.1 lists the comparison stars we used to calibrate the REM optical photometry (see Section 2.2). Table A.2 lists the derived g’r’i’z’ magnitudes for GM Aur.

Appendix B: Line profile correlation matrices

This section provides correlation matrices between line profiles (Johns & Basri 1995). Correlation coefficients are computed across the line profiles for every pair of velocity channels common to the two profiles. Fig. B.1 to B.4 present the correlation matrices of the optical and near-infrared line profiles studied here.

Fig. B.1. Correlation matrices for hydrogen optical line profiles computed from the 15 OHP/SOPHIE spectra obtained during the campaign: Hα ⋆ Hα (left), Hβ ⋆ Hβ (center), and Hα ⋆ Hβ (right).

Fig. B.2. Correlation matrices for the hydrogen near-infrared line profiles computed from the 34 CFHT/SPIRou spectra obtained during the campaign : Paβ ⋆ Paβ (left), Brγ ⋆ Brγ (center), and Paβ ⋆ Brγ (right).

Fig. B.3. Correlation matrices between optical and near-infrared hydrogen line profiles computed from ten OHP/SOPHIE and ten CFHT/SPIRou spectra obtained over the same nights during the October runs: Hα ⋆ Paβ (left), Hβ ⋆ Paβ (center left), Hα ⋆ Brγ (center right), and Hβ ⋆ Brγ (right).

Fig. B.4. Correlation matrices for the HeI 10830 Å and hydrogen lines: HeI ⋆ HeI (left), Paβ ⋆ HeI (center), and Hβ ⋆ HeI (right). The Brγ ⋆ HeI and Hα ⋆ HeI matrices are not shown as they are similar to those of Paβ ⋆ HeI and Hβ ⋆ HeI, respectively.

Appendix C: EWs and J-band measurements from the ExTrA spectra

We provide the HeI 10830 Å, Paβ, Paγ, and Paδ line EWs and the J-band photometry measured from the ExTrA spectra in Table C.1. For each night, the table lists the mean observation time, the median EW measurement and its standard deviation for each spectral line, and the J-band magnitude and its error.

Appendix D: nIR veiling and line EWs

We present in Figure D.1 the evolution of near-infrared veiling in the JHK bands and of the HeI, Paβ, and Brγ line EWs over the course of the campaign.

Fig. D.1. Veiling measured on SPIRou spectra in the JHK bands (triangles) and the EW of the HeI, Paβ, and Brγ lines (dots) plotted as a function of Julian date (top) and rotational phase (bottom).

Appendix E: Line profile variability on successive SPIRou runs

Figures E.1 to E.4 show the line profile variability of the HeI 10830 Å, Paβ, and Brγ line profiles for the successive SPIRou runs in September, October, November, and December 2021. The figures also include 2D periodograms for each line. We note, however, that only the October SPIRou run is long enough, extending over 14 days, to yield significant results when searching for a modulation of the line profiles around the rotational period of 6.04 days. The September, November, and December runs lasted for only 9, 6, and 9 days, respectively, which is too short to reliably investigate periods longer than 5 days.

Fig. E.1. Near-infrared HeI (left), Paβ (center), and Brγ (right) line profiles obtained over nine days during the September 2021 SPIRou run. Top: Line profiles are plotted as a function of Julian date (left subpanels) and rotational phase (right subpanels). The colors represent successive rotational cycles. Bottom: 2D periodograms across the line profiles. The dotted horizontal red line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve displays the mean line profile. The color code reflects the periodogram power from zero (blue) to 1 (red).

Fig. E.2. Near-infrared HeI (left), Paβ (center), and Brγ (right) line profiles obtained over 14 days during the October 2021 SPIRou run. Top: Line profiles are plotted as a function of Julian date (left subpanels) and rotational phase (right subpanels). The colors represent successive rotational cycles. Bottom: 2D periodograms across the line profiles. The dotted horizontal red line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve displays the mean line profile. The color code reflects the periodogram power from zero (blue) to 0.9 (red).

Fig. E.3. Near-infrared HeI (left), Paβ (center), and Brγ (right) line profiles obtained over six days during the November 2021 SPIRou run. Top: Line profiles are plotted as a function of Julian date (left subpanels) and rotational phase (right subpanels). The colors represent successive rotational cycles. Bottom: 2D periodograms across the line profiles. The dotted horizontal red line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve displays the mean line profile. The color code reflects the periodogram power from zero (blue) to 1 (red).

Fig. E.4. Near-infrared HeI (left), Paβ (center), and Brγ (right) line profiles obtained over nine days during the December 2021 SPIRou run. Top: Line profiles are plotted as a function of Julian date (left subpanels) and rotational phase (right subpanels). The colors represent successive rotational cycles. Bottom: 2D periodograms across the line profiles. The dotted horizontal red line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve displays the mean line profile. The color code reflects the periodogram power from zero (blue) to 1 (red).

All Tables

All Figures

Fig. 1. Temporal sampling of the GM Aur campaign from September 6, 2021, to March 8, 2022. Bottom: g′-band light curve from LCOGT (black dots) and from REM (magenta crosses). The REM g′-band magnitudes are offset in this figure by +0.1 mag to match the LCOGT measurements. The mean photometric error on both the LCOGT and REM g′-band measurements is 0.025 mag. Vertical arrows: the vertical arrows show the dates of OHP/SOPHIE (red), CFHT/SPIRou (blue), and ESO/ExTrA (black) observations. The core of the campaign took place during October 2021 with contemporaneous measurements from the five instruments.

In the text

Fig. 2. GM Aur light curves. Top: GM AUr’s u′g′r′i′ light curves from LCOGT observations that extended over nearly four months. The mean photometric error is 0.025 mag in the g′r′i′ bands and 0.033 mag in the u′ band. A scaled TESS light curve (black) obtained contemporaneously is overplotted onto the LCOGT r′-band light curve. Middle: LCOGT light curves folded in phase with a period of 6.04 days with the ephemeris of Eq. (1). Bottom: same as above, with a color scale for data points that reflects the Julian date. The amplitude of the light variations appears to increase slightly toward the end of the observing run. In all panels, the error bars on the measurements are smaller than the symbol size.

In the text

	Fig. 3. Period search results. Top: CLEAN periodogram analysis of the u′g′r′i′ light curves. A peak occurs at the frequency of 0.166 day−1, which corresponds to a period of 6.04 ± 0.15 days. Bottom: string-length analysis of the u′g′r′i′ light curves. A clear minimum appears for a period of 6.03 ± 0.15 days.
In the text

	Fig. 4. Color-magnitude relation. The (u′−g′), (g′−r′), and (r′−i′) colors are plotted as a function of the system brightness in the g′-band. The dashed lines indicate the expected ISM reddening slope computed from Fiorucci & Munari (2003) with R = 3.1.
In the text

Fig. 5. Optical veiling measurements from the OHP/SOPHIE spectra. Top: part of a spectral window showing the template spectrum (gray), the observed spectrum (black), and the velocity broadened and veiled template spectrum (red) that fits the observed spectrum. Bottom: veiling measured in several spectral windows (see text) as a function of Julian date. The central wavelength of the spectral windows is indicated in the top left corner of the panel.

In the text

	Fig. 6. Mean optical veiling plotted as a function of rotational phase. The color code indicates the Julian date.
In the text

Fig. 7. LiI 6707 Å line profile. Left: the 34 line profile measurements from SOPHIE spectra are shown superimposed. The color code corresponds to successive rotational cycles. Right: 2D periodogram across the line profile. The dotted horizontal line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve displays the mean line profile. The color code reflects the periodogram power.

In the text

Fig. 8. Optical line profile variability. Top: series of optical line profiles Hα, Hβ, and HeI plotted as a function of Julian date (left subpanels) and rotational phase (right subpanels). The color code corresponds to successive rotational cycles. Middle: same profiles, superimposed to illustrate their variability. Bottom: 2D periodograms across the line profiles. The color code reflects the periodogram power, from zero (blue) to 0.6 (red). The dotted horizontal red line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve is the mean line profile.

In the text

Fig. 9. Radial velocity of the narrow component of the HeI 5876 Å line profile in the stellar rest frame plotted as a function of Julian date (top) and rotational phase computed from the ephemeris of Eq. (1), (bottom). The color code corresponds to successive rotational cycles. The fit by a geometrical accretion shock model (see text) is shown (dash-dotted curve) together with its 1σ uncertainty (gray area).

In the text

	Fig. 10. Near-infrared line profiles: HeI (left), Paβ (center), and Brγ (right), plotted with arbitrary offsets as a function of Julian date. Each color represents a SPIRou run, namely September (blue), October (orange), November (green), December 2021 (red), and January 2022 (purple). The October SPIRou run is contemporaneous to the OHP/SOPHIE observations.
In the text

Fig. 11. Near-infrared line profile variability. Top: series of residual near-infrared line profiles HeI (left), Paβ (center), and Brγ (right), plotted superimposed to illustrate their variability. Each color represents a SPIRou run as in Fig. 10. Bottom: 2D periodograms across the line profiles. The color code reflects the periodogram power from zero (blue) to 0.5 (red). The dotted red horizontal line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve displays the mean line profile.

In the text

Fig. 12. EW of near-infrared lines. Top: EW of the HeI, Paβ, and Brγ lines measured on SPIRou spectra plotted as a function of Julian date. The measurement uncertainties are about the symbol size. Bottom: difference between the EWs measured by direct line profile integration (i.e., including the redshifted absorption component) and the EW derived from the Gaussian fitting of the emission component only plotted as a function of Julian date. This differential quantity measures the strength of the redshifted absorption component in the Paβ and Brγ line profiles. The more negative the differential quantity, the deeper the redshifted absorption.

In the text

	Fig. 13. Radial velocity variations measured in the near-infrared photospheric lines of the SPIRou spectra. Top: radial velocity as a function of Julian date. Bottom: radial velocity curve folded in phase at the stellar rotational period P = 6.04 days. The color code reflects the Julian date.
In the text

	Fig. 14. Median spectrum for each night from the ExTrA-T2 telescope (lower part, 88 spectra) and ExTrA-T3 telescope (upper part, 69 spectra). The main emission lines are indicated. The color is proportional to the EW of the Paβ line shown in the top panel as a function of Julian date.
In the text

	Fig. 15. EW variability of the near-infrared line profiles. Left: EW measurements plotted as a function of Julian date for the HeI, Paβ, Paγ, and Paδ lines from the ExTrA spectra. Right: GLS periodogram of the EW measurements. The period on the x-axis is displayed on a log scale.
In the text

	Fig. 16. EWs of near-infrared line profiles and J-band magnitude variability. Top: EW of the near-infrared HeI and Paβ lines folded in phase at the stellar rotational period. EW(HeI) is offset by −10 Å for clarity. Bottom:J-band light curve, deduced from ExTrA spectra, folded in phase at the stellar period. The brightness modulation is similar to that observed at optical wavelengths (see Fig. 2).
In the text

	Fig. 17. Mid-term variation of EW(Paβ), (red) compared to the system u′-band light curve (blue) for measurements taken less than one day apart. To facilitate the comparison, the EW measurements are plotted on a magnitude scale and are offset, namely −2.5 log EW(Paβ) + 5.
In the text

Fig. 18. Magnetospheric ejection model. Development of a magnetospheric ejection computed from star-disk interaction MHD simulations. The snapshots shown here are extracted from model 3 of Pantolmos et al. (2020), where the magnetospheric truncation radius amounts to 54% of the corotation radius. The three snapshots are shown at 0.2, 0.5, 0.8, and 1.1Prot. The white curves indicate expanding magnetic field lines that give rise to the ejection of plasmoids. The color scale indicates density. The green arrows show the velocity field of the stellar and disk winds and of MEs at their interface.

In the text

	Fig. 19. Magnetospheric ejection model. Velocity of the gas in the ejected plasmoid at a distance of 10 (magenta), 20 (blue), and 30 (cyan) stellar radii as a function of time in units of the stellar rotational period. Successive ejections of plasmoids are featured. The vertical lines correspond to the snapshots shown in Fig. 18.
In the text

	Fig. B.1. Correlation matrices for hydrogen optical line profiles computed from the 15 OHP/SOPHIE spectra obtained during the campaign: Hα ⋆ Hα (left), Hβ ⋆ Hβ (center), and Hα ⋆ Hβ (right).
In the text

	Fig. B.2. Correlation matrices for the hydrogen near-infrared line profiles computed from the 34 CFHT/SPIRou spectra obtained during the campaign : Paβ ⋆ Paβ (left), Brγ ⋆ Brγ (center), and Paβ ⋆ Brγ (right).
In the text

	Fig. B.3. Correlation matrices between optical and near-infrared hydrogen line profiles computed from ten OHP/SOPHIE and ten CFHT/SPIRou spectra obtained over the same nights during the October runs: Hα ⋆ Paβ (left), Hβ ⋆ Paβ (center left), Hα ⋆ Brγ (center right), and Hβ ⋆ Brγ (right).
In the text

	Fig. B.4. Correlation matrices for the HeI 10830 Å and hydrogen lines: HeI ⋆ HeI (left), Paβ ⋆ HeI (center), and Hβ ⋆ HeI (right). The Brγ ⋆ HeI and Hα ⋆ HeI matrices are not shown as they are similar to those of Paβ ⋆ HeI and Hβ ⋆ HeI, respectively.
In the text

	Fig. D.1. Veiling measured on SPIRou spectra in the JHK bands (triangles) and the EW of the HeI, Paβ, and Brγ lines (dots) plotted as a function of Julian date (top) and rotational phase (bottom).
In the text

Fig. E.1. Near-infrared HeI (left), Paβ (center), and Brγ (right) line profiles obtained over nine days during the September 2021 SPIRou run. Top: Line profiles are plotted as a function of Julian date (left subpanels) and rotational phase (right subpanels). The colors represent successive rotational cycles. Bottom: 2D periodograms across the line profiles. The dotted horizontal red line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve displays the mean line profile. The color code reflects the periodogram power from zero (blue) to 1 (red).

In the text

Fig. E.2. Near-infrared HeI (left), Paβ (center), and Brγ (right) line profiles obtained over 14 days during the October 2021 SPIRou run. Top: Line profiles are plotted as a function of Julian date (left subpanels) and rotational phase (right subpanels). The colors represent successive rotational cycles. Bottom: 2D periodograms across the line profiles. The dotted horizontal red line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve displays the mean line profile. The color code reflects the periodogram power from zero (blue) to 0.9 (red).

In the text

Fig. E.3. Near-infrared HeI (left), Paβ (center), and Brγ (right) line profiles obtained over six days during the November 2021 SPIRou run. Top: Line profiles are plotted as a function of Julian date (left subpanels) and rotational phase (right subpanels). The colors represent successive rotational cycles. Bottom: 2D periodograms across the line profiles. The dotted horizontal red line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve displays the mean line profile. The color code reflects the periodogram power from zero (blue) to 1 (red).

In the text

Fig. E.4. Near-infrared HeI (left), Paβ (center), and Brγ (right) line profiles obtained over nine days during the December 2021 SPIRou run. Top: Line profiles are plotted as a function of Julian date (left subpanels) and rotational phase (right subpanels). The colors represent successive rotational cycles. Bottom: 2D periodograms across the line profiles. The dotted horizontal red line drawn at a frequency of 0.166 day−1 indicates the stellar rotational period. The white curve displays the mean line profile. The color code reflects the periodogram power from zero (blue) to 1 (red).

In the text