© ESO, 2012

1. Introduction

Classical T Tauri stars (CTTSs) are young, magnetically active low-mass (M < 2 M_⊙) stars that display signs of accretion from a circumstellar disk. They show an emission excess with respect to the photosphere at wavelengths from X-rays to the radio and are both spectroscopically and photometrically variable (Bouvier et al. 2007). The CTTSs have strong permitted emission lines that vary on short timescales (days), a redshifted absorption component associated with the accreting material, and a blueshifted absorption component due to winds. They also show forbidden emission lines, which originate from a wind/jet that is powered by the accretion process. Magnetospheric accretion models have been successful in explaining most of the observed characteristics of CTTSs (Shu et al. 1994; Hartmann et al. 1994; Muzerolle et al. 2001; Kurosawa et al. 2006; Lima et al. 2010) and are currently the most accepted ways to describe these systems. In the past decade, magnetohydrodynamic (MHD) simulations have followed the rotational evolution of CTTSs over as many as hundreds of rotational periods, presenting a dynamical view of the star-disk interaction (Goodson et al. 1999; Romanova et al. 2002; Long et al. 2008; Zanni & Ferreira 2009) that can be tested with synoptic observations.

Accretion has a long-lasting role in early stellar evolution, by providing mass and helping to regulate the angular momentum transfer from the star to the disk. The accretion process also affects the disk evolution, since the high energy radiation produced in UV and X-rays influences the circumstellar disk lifetime (Gorti & Hollenbach 2009; Owen et al. 2010). Furthermore, it has been suggested that accretion may affect the early evolution of the star itself, by modifying the stellar radius at young ages (Baraffe et al. 2009). Therefore, it is important to understand the effects of the accretion process, and CTTSs are excellent laboratories for such a study.

V2129 Oph is a K5 CTTS (T_eff ≈ 4500 K, Donati et al. 2007) that has a maximum visual brightness of m_V = 11.2 (Grankin et al. 2008). It is a young system (2–3 Myr) located in the ρ Oph star-forming region at 120 ± 5 pc (Loinard et al. 2008) and seen at moderate inclination with respect to our line of sight (i ~ 60°). The measured values of rotational periods vary in the literature from 6.35 to 6.6 days (Grankin et al. 2008), which is indicative of differential rotation. V2129 Oph has been the focus of a large observing campaign, including simultaneous or quasi-simultaneous optical and near-IR photometry, high-resolution spectroscopy and spectropolarimetry, and Chandra X-ray observations. The spectropolarimetric data that we use in the present paper, which were obtained with ESPaDOnS, were previoulsy analyzed by Donati et al. (2011). These authors obtained the magnetic field structure at the surface of the star for two observing epochs (in 2005 and 2009), and found that it is composed of a dipole and an octupole tilted by about 20° with respect to the rotation axis. The magnetic field components had intensities that varied between the two epochs, which implies that the magnetic field has a non-fossil origin. Romanova et al. (2011) computed a numerical 3D MHD model of the 2005 magnetic field configuration of V2129 Oph and simulated accretion onto the star. They showed that the disk is truncated by the dipole component and that accretion proceeds towards the star in two main accretion streams. Close to the stellar surface, the flow is redirected by the octupolar component. Argiroffi et al. (2011) analyzed Chandra X-ray observations obtained in 27–29 June 2009, simultaneously with the HARPS data analyzed in the present paper. They showed that the soft X-ray emission observed in V2129 Oph, corresponding to dense cool plasma of a few MK, comes from material heated in the accretion shock. The variability of this soft X-ray emission in the observations is attributed to changes in the viewing angle of the accretion shock, as the system rotates.

We present here an analysis of HARPS, ESPaDOnS, and SMARTS echelle spectra, obtained over a period of 73 days, together with simultaneous photometric measurements. Most of our spectra were taken from June 10 to July 13 2009, providing a coverage of almost four rotational cycles. Although the ESPaDOnS data were previously analyzed to determine the magnetic field configuration of the system by Donati et al. (2011), we use them here to complete our dataset in the time domain and to analyze the structure of the magnetospheric accretion flow. We describe the accretion diagnostics and propose a viable scenario for the magnetospheric configuration and the star-disk interaction. We compare the observed Balmer emission lines with radiative transfer models calculated for the funnel flow structure obtained from 3D MHD simulations, using the 2009 magnetic field configuration. We also calculate Hα line profiles with a hybrid magnetospheric and diskwind model to investigate the importance of the wind to the observed line profiles.

2. Observations

2.1. Spectroscopy

We present in Table 1 a journal of spectroscopic observations. Some of the spectroscopic observations of V2129 Oph were carried out from June 10 to June 30, 2009, at ESO, La Silla with the HARPS dual fiber echelle spectrograph (Mayor et al. 2003), covering the 3800 Å to 6900 Å spectral domain at a spectral resolution of λ/Δλ ~ 80 000. We acquired 28 high-resolution spectra, over 20 non-consecutive nights at the 3.6 m telescope. The data were automatically reduced by the HARPS Data Reduction Software. The reduction procedure includes optimal extraction of the orders and flat-fielding, wavelength calibration, and the removal of cosmic rays. Spectropolarimetric observations were collected from July 01 to July 14, 2009, with ESPaDOnS on CFHT. The reduction of these ESPaDOnS data is described in detail in Donati et al. (2011). Spectra were also obtained with the SMARTS 1.5 m telescope at CTIO (Chile). Two spectrographs were used. The low-to-moderate long slit RC spectrograph was used with the 300′′ long slit, oriented east-west. The detector was a Loral 1K CCD, with a 260 × 1199 readout region. A total of 13 spectra were obtained from May 12 (JD 2 454 964.81) to July 7, 2009 (JD 2 455 020.73) with different settings, covering either the Hβ or Hα spectral regions. We obtained three images at each epoch in order to filter cosmic rays. Each set of images was accompanied by a wavelength calibration exposure of a Ne-Ar or Th-Ar arc lamp. The image were bias-subtracted, trimmed, and flattened using dome flats obtained each night. We co-added the three images using a median filter. We also obtained echelle spectra with the 4.0 m Cassegrain echelle, now bench-mounted and fiber-fed from the SMARTS 1.5 m telescope. A total of 19 spectra were obtained from May 1 (JD 2 454 953.93) to July 10, 2009 (JD 2 455 023.71). A 100 micron slit was used, resulting in a resolution of about 30 000. The 31.6 l/mm echelle was used with a 226 l/mm cross-disperser and the detector was a 2k SITe CCD with an ARCON controller. Full wavelength coverage was achieved from 4020 Å to 7300 Å. We obtained three 20 min exposures at most epochs, and median-filtered the data to reject cosmic rays. We located the orders and extracted the spectra by fitting a Gaussian perpendicular to the dispersion direction at each pixel. We determined the precise wavelength solution by cross-correlating the Th-Ar images obtained along with the spectra, with the template solution. We corrected the velocities to the barycenter of the solar system.

2.2. Photometry

The journal of photometric observations is presented in Table 2. Our BVRIJHK photometric observations were carried out with the ANDICAM dual-channel imager at the SMARTS 1.3 m telescope at Cerro Tololo. The optical data is composed of 10 sets of observations on 10 nights, between 23 June (JD 2 455 007.57) to 8 July 2009 (JD 2 455 021.64). BKLT J162749-242540 (ROC 31) was used as a comparison star. The near-IR data consists of 13 sets of observations on 13 nights, between 21 June (JD 2 455 005.65) and 8 July 2009 (JD 2 455 021.65). Integration times were 4 s (the minimum allowed). We obtained 3 images through each filter at different dither positions that were shifted and co-added. Measurements are differential with respect to the star BKLT J162746-242323 (Elias 2–35; VSSG 13), which was a highly reddened background K5III star viewed through the dark cloud. Published 2MASS magnitudes of the comparison star were used to convert the instrumental magnitudes to apparent magnitudes.

Seven sets of V and R observations on 7 nights were obtained at Palermo in the period between June 22 (JD 2 455 006.40) and June 29 (JD 2 455 013.40). SR21 was used as a comparison star. We also obtained 15 UBVRI observations at Mount Maidanak in 7 nights from 19 June (JD 2 455 002.29) to 30 June 2009 (JD 2 455 013.24).

3. Results

In our analysis, we use the ephemeris published by Donati et al. (2007) for V2129 Oph in the phase calculations $\mathrm{HJD}\mathrm{=}\mathrm{2}\hspace{0.17em}\mathrm{453}\hspace{0.17em}\mathrm{540.0}\mathrm{+}\mathrm{6.53}\mathit{E,}$(1)where 6.53 days is the adopted period of V2129 Oph.

3.1. Photometric and spectral variations

The photometric data show rotational modulation in all filters (BVRI), with an amplitude decreasing from the B to the I band. The V data are presented in Fig. 1 (panel a), after correction for the interstellar extinction with A_V = 0.6 (Donati et al. 2007).

Photospheric radial velocities were measured on HARPS and ESPaDOnS spectra using the spectral region 5415 Å  < λ < 5465 Å, cross-correlating V2129 Oph spectra with spectra of the K7 weak T Tauri star (WTTS) V819 Tau (Fig. 1, panel b). The measured values from both spectrographs are in excellent agreement, and the values obtained for the HARPS spectra agree with those directly provided by the HARPS pipeline by cross-correlating the object’s spectrum with a K5 mask. The error in the HARPS and ESPaDOnS photospheric radial velocities is around 0.2  km s^-1 and 0.1  km s^-1, respectively, obtained by fitting the cross-correlation and taking into account the signal-to-noise ratio (S/N) of the spectra. However, looking at Fig. 1 (panel b) we can see that the HARPS data (filled symbols) are clearly less scattered than the ESPaDOnS ones (open symbols), which is probably due to the higher stability of the HARPS spectrograph.

We present in Figs. 2 and 3 the Hα and He i line profiles observed with HARPS, and folded in phase. The Hα line profile varies substantially in phase, having the same characteristics in similar phases, such as a red emission shoulder from phase 0.05 to 0.35 and a redshifted absorption component from phase 0.6 to 0.8.

The He i (5876 Å) line contains only a narrow component in the HARPS and ESPaDOnS spectra, which we fitted with a single Gaussian to measure the emission line parameters. The Gaussian central wavelength position was used to measure the He i radial velocity, which is clearly modulated over the rotational phase (Fig. 1, panel d). The HARPS and ESPaDOnS measurements are consistent and stable over a  ~33-day timescale. The error in the measured radial velocities is on the order of 0.5 to 2  km s^-1 and dominated by our fitting uncertainty in the emission line profile.

Fig. 1 V band photometry corrected from extinction (panel a)). SMARTS data are shown as squares, Palermo data as stars, and Maidanak data as triangles. Photospheric radial velocities (panel b)), veiling (panel c)), and He i (5876 Å) radial velocity (panel d)) phase variations. In panels b)–d), filled symbols are from HARPS data and open ones from ESPaDOnS. The error bars are shown for each observation in panel a), and are smaller than the symbol sizes in panel b). We show in panels c and d only the mean error bars to avoid crowding the figures. Different colors represent different rotational cycles identified as the integer part of the observational phase in Table 1: black – cycle 0, green – cycle 1, red – cycle 2, aqua – cycle 3, dark blue – cycle 4, pink – cycle 5.

3.2. Veiling

We measured veiling values in the spectral region 5415 Å  < λ < 5465 Å. This region contains data with a good S/N and many photospheric lines that are not blended, and is also devoid of emission lines. We used the projected rotational velocity of V2129 Oph obtained by Donati et al. (2007), vsini = 14.5 ± 0.3  km s^-1, which agrees with our determination. We initially used V819 Tau (K7) as a comparison star, but owing to its spectral type mismatch with V2129 Oph, we checked our measured veiling values with the WTTS V410 Tau (K5).

Fig. 2 Hα profiles obtained with the HARPS spectrograph. The profiles are normalized to the continuum level and the vertical dotted lines indicate the stellar rest velocity. Phases and JDs are given for each spectrum in the panels. The color code is the same as in Fig. 1.

Fig. 3 He i 5876 Å profiles obtained with the HARPS spectrograph. The profiles are normalized to the continuum level and the vertical dotted lines indicate the stellar rest velocity. Phases and JDs are given for each spectrum in the panels. The color code is the same as in Fig. 1.

While V819 Tau rotates slowly (vsini < 15  km s^-1, Herbig & Bell 1988) and is an excellent standard star to use in the determination of both vsini and veiling, V410 Tau is a rapid rotator (vsini = 74 ± 3  km s^-1, Skelly et al. 2010). We first calculated the veiling of V410 Tau by using V819 Tau as a standard, obtaining the value of 0.11 ± 0.04. This is not a true veiling, since V410 Tau is not accreting, but it gives us an idea of the veiling value created by the mismatch of spectral types between a K5 star and a K7 standard. We then rotationally broadened the spectra of V2129 Oph to the rotation velocity of V410 Tau and calculated the veiling of V2129 Oph with V410 Tau as a standard star. The values obtained are consistent with those determined using V819 Tau as a standard, but shifted to lower values by 0.15 ± 0.04, which is consistent with the veiling induced by the spectral type mismatch between V819 Tau and V2129 Oph.

The veiling error obtained using V819 Tau as a standard is on the order of 0.01 and 0.02 for the HARPS and ESPaDOnS data, respectively, while it reaches 0.04 when we use V410 Tau, owing to its high rotational velocity. We therefore decided to use the veiling values obtained with V819 Tau that had been decreased by a factor of 0.15 to take into account the spectral type difference between V819 Tau and V2129 Oph. Although the absolute values still have an error of about 0.04 due to the scaling process, the relative error between each measurement is on the order of 0.02. The final result is presented in Fig. 1 (panel c). The veiling that we determined for the ESPaDOnS data are compatible, within the errors, with the values calculated by Donati et al. (2011).

As shown in Sect. 4.1, a large cold spot is actually located in the stellar photosphere, which should affect the measured veiling values, as the system’s continuum flux is modulated by the spot. For Δ(V) ~ 0.21 mag, the relative veiling variation amounts to about 20%, which is, however, much lower than the veiling variations we observe.

Fig. 4 Rotational variation in the He i EW and the mass accretion rate. The mass accretion rate errors come from the relation between the He i (5876 Å) luminosity and the accretion luminosity by Fang et al. (2009). Filled symbols are from HARPS data and open ones from ESPaDOnS. Different colors represent different rotational cycles identified as the integer part of the observational phase in Table 1: black – cycle 0, green – cycle 1, red – cycle 2, aqua – cycle 3, dark blue – cycle 4, pink – cycle 5.

3.3. Mass accretion rate

We measured the He i 5876 Å line equivalent widths (EWs) (see Table 1 and Fig. 4, top panel) and computed line fluxes using the V band photometry. We do not have simultaneous V measurements for all the spectroscopic data, but the phase coverage of the V data is quite complete (Fig. 1, panel a). We therefore interpolated the V magnitude folded in phase at the rotational phase of the spectroscopic observations. Using the empirical correlations of Fang et al. (2009), we calculated the mass-accretion rates values from the He i line fluxes. The results are shown in Fig. 4, where we see that the mass accretion rate remains fairly stable most of the time at a value of (1.5 ± 0.6) × 10^-9 M_⊙ yr^-1, which agrees with the mass accretion rate determined by Donati et al. (2011). However, we note that, around phase 0.5, the mass accretion rate varied substantially over a period of days during the HARPS observations. This phase corresponds to the main accretion spot facing the observer. There is further evidence from the Hα and Hβ emission line profiles that the mass accretion increased during the HARPS observations (see Sect. 4.2).

4. Discussion

4.1. Hot and cold spots

The veiling clearly varies in phase during the HARPS observations (Fig. 1, panel c), with a maximum value of  ~0.30 close to phase φ = 0.5 and a minimum value generally above zero, suggesting that the hot spot remains in view during the whole rotational cycle. This is confirmed by the analysis of line profile variations (see Sect. 4.3). We developed a hot spot model with a circular spot to reproduce the veiling modulation. We used T_eff = 4500 K, T_spot = 8000 K and varied the spot radius, latitude, and location with rotational phase. We run models with two different values of the inclination of the system with respect to our line of sight, i = 45° and i = 60°, which correspond to values suggested in the literature for the system (Donati et al. 2011). For a system inclination of i = 60°, the best-fit model corresponds to a hot spot with a filling factor of 1.2%, located at phase 0.55 and high latitude (78°). If we use, instead, a system inclination of i = 45°  we obtain a hot spot with a filling factor of 0.3%, located at phase 0.55 and latitude of 70.5°. The two solutions, which are indicated as black lines in Fig. 5, are degenerate, having the same χ² and, although they do not constrain the system inclination, they both converge to a major hot spot with very similar characteristics. The low observed veiling contrast during the rotational cycle can only be reproduced with a very high-latitude spot, i.e., close to the rotational axis. However, none of the models perfectly reproduce the narrow veiling peak near phase 0.5, but instead predict a much smoother modulation. We have overplotted in red in Fig. 5 the veiling variations generated by the hot spot models that most closely reproduce the observed He i 5876 Å radial velocity curve (Fig. 6). The hot spot inferred from the He i radial velocity variations is consistent overall with the observed veiling changes. We also note that the veiling variability was not seen in all the observed rotational cycles (Fig. 1), which suggests that a hot spot has evolved and dimmed on a timescale of about a week. Very low veiling values are indeed derived from the ESPaDOnS spectra, which were acquired the week after the HARPS measurements.

Fig. 5 Best-fit hot spot model for the HARPS veiling variations (black lines). The red lines represent the veiling generated by the optimal hot spot fits to the He i radial velocity variations (Fig. 6). Models with i = 60° and i = 45° are shown as solid and dashed lines, respectively . The mean error bar in the observed veiling is shown on the left and the color code of the points is the same as in Fig. 1.

The He i 5876 Å line is thought to be produced in the post-shock region of the accretion shock. The modulation of its radial velocity (Fig. 1, panel d) is then likely due to the hot spot rotation. The phase of the He i radial velocity curve is indeed exactly the opposite of that of the photospheric radial velocity curve, the bluest velocity being reached at φ ~ 0.35 and the reddest at φ ~ 0.6. We assumed that all the He i line is produced at the location of the hot spot and fitted the measured radial velocities of He i with a hot spot model. We set the spot latitude and phase, the infall velocity, and the system inclination as free parameters. However, the system inclination was not well-constrained in the fitting procedure, so we decided to compute models with inclination values fixed at i = 45° and i = 60°, as we did in the veiling model calculations. The best-fit results are shown in Table 3 and in red in Fig. 6. The final values for each parameter obtained with i = 45° and i = 60° agree within the errors and correspond to a hot spot located at high latitude (80°) near phase 0.6 with an infall velocity of about 10  km s^-1. This is consistent with the assumption that the hot spot required to account for the veiling variations is also responsible for the modulation of the He i radial velocity. To illustrate this, we overplotted in Fig. 6, in black, the He i radial velocity model that corresponds to the best fit of the veiling curve. That the He i radial velocity changes during the whole rotational cycle indicates that the hot spot is always in view, which agrees with the veiling analysis. The small difference in spot parameters required to explain the veiling variations and the He i radial velocity curve are probably related to the basic assumption of circular spots at the stellar surface, an assumption that is probably wrong, but it would be beyond the scope of this paper to investigate more complex spot geometries.

The photometric variations appear to be dominated by a low-level modulation from a cool spot (Fig. 1, panel a) and the smooth shape of the light curve suggests that a high latitude spot is always in view. A cold spot model was run with a fixed inclination at either i = 45° or i = 60° to attempt to reproduce the observed light curves, but the model solutions are unable to efficently constrain all the free parameters, owing to the small number of data points available. The best-fit models with i = 45° correspond to a cold spot with temperature T_spot = (4200 ± 100) K located at phase 0.57 ± 0.07 and latitude of 60° ±  20°. Unfortunately, the cold spot radius is almost unconstrained (R_spot > 35°), which corresponds to a large incertitude in the spot filling factor (f > 9%). We show in Fig. 7 the χ² contours of the parameters adjusted in the cold spot models. In Fig. 8, we present the fits to the photometric data of the best-fit cold spot model with i = 45°. The best-fit models with i = 60° correspond to a cold spot with temperature T_spot = (3500 ± 500) K located at phase 0.57 ± 0.07 and latitude of 84° ±  5°. The cold spot radius is again unconstrained (R_spot > 50°, f > 18%), although a large spot is clearly needed in both models.

Fig. 6 HARPS (filled symbols) and ESPaDOnS (open symbols) He i radial velocity measurements, with respective error bars. The red curves represent the best-fit hot-spot model for the He i radial velocity data obtained with i = 60° (solid line) and i = 45° (dashed line). The black solid line corresponds to the best fit to the veiling variations with i = 60°, presented in Fig. 5.

Fig. 7 χ2 contours of the cold spot models with i = 45°. Each contour level corresponds to an increase of χ2 minimum. The red lines correspond to a 2χ2 minimum and represent the error in each parameter quoted in the text.

As pointed out by Petrov et al. (2011), hot and cold spots can distort in similar ways the photospheric line profiles, and photometry taken simultaneously with the spectra is needed to distinguish between the two effects. In the case of V2129 Oph, the photometric variations clearly indicate that there is a large cold spot in the stellar photosphere that would cause the photospheric radial velocity modulation seen in Fig. 1 (panel b). Moreover, the photospheric radial velocity variations are very stable over the observed period and clearly present in all rotational cycles, including those with very low veiling values (Fig. 1, panel b). This also points to a cold spot as the most likely cause of the photospheric radial velocity variations, instead of a hot spot. The reddest velocity is seen at φ ~ 0.35, the bluest at φ ~ 0.6, and the mean velocity of  ~−7  km s^-1 is crossed at φ ~ 0.5, indicating that the spot then faces the observer. This is consistent with the assumption that the cold spot required to account for the BVRI light curves is also responsible for the photospheric radial velocity modulation. We then tried to reproduce the observed HARPS photospheric radial-velocity variations using cold spot models with a system inclination of 60°. We varied the spot brightness, its latitude, the phase at which it occurs, and the systemic velocity of the star, but we could not find a unique solution to the data. We calculated 10 000 solutions exploring the parameter space of possible solutions and obtained probability density functions for each parameter. Some parameters are well-constrained, such as the phase where the spot faces the observer (φ = 0.475 ± 0.025), while others are clearly poorly constrained, like the spot filling factor (9%  < f < 50%), its brightness (0.0 <  brightness  < 0.8), and the spot latitude (0°  <  latitude  < 55°). The solution is degenerated for several combinations of the poorly constrained parameters. We present in Fig. 9 the radial velocity variations due to a cold spot calculated in each of the 10 000 models run at the phases of the HARPS observations and a system inclination of i = 60°.

Fig. 8 Best-fit cold-spot model with i = 45° to the photometric observations. The parameters of the cold-spot model used to generate the above light curves correspond to a χ2 minimum (Tspot = 4200 K, spot latitude of 60°, spot phase of 0.57, and spot filling factor of 40%).

The hot and cold spots in V2129 Oph are both located at high latitudes and appear at about the same rotational phases, in agreement with the Zeeman-Doppler results (Donati et al. 2011), but the hot spot, which produces the veiling and some of the emission lines, is most likely located at the chromospheric level, more or less spatially overlapping the major cool spot at the photospheric level.

Fig. 9 Radial velocity variations induced by a cold spot model obtained in 10 000 Monte Carlo runs varying the parameter space of possible solutions (small black dots). The triangles represent the median values of the theoretical solutions and the HARPS photospheric radial velocities are shown as large filled circles.

Fig. 10 Hα (top) and Hβ (bottom) profiles obtained with the HARPS (solid), ESPaDOnS (dashed), and SMARTS (dash-dotted, Hα only) spectrographs. The spectra correspond to a mean profile around the phases given in the panels. The profiles are normalized to the continuum level and the vertical dotted lines indicate the stellar rest velocity.

4.2. Emission lines

V2129 Oph displays Balmer and He i 5876 Å lines in emission, as usually observed in CTTSs (Figs. 2 and 3). The Balmer emission lines contain a major emission component with a rather triangular shape that is very difficult to decompose with simple Gaussian profiles. Blue and red emission bumps appear as the star rotates and a redshifted absorption component is seen in both Hα and Hβ, in many phases in the second half of the rotational cycle. The Balmer emission-line profiles vary substantially in phase but also exhibit a strong phase coherence, presenting the same characteristics in similar phases, as can be seen in Fig. 10 where we show the HARPS, ESPaDOnS, and SMARTS Hα and Hβ emission-line profiles averaged at a few chosen phases.

The Hα, Hβ, and He i lines display periodic variability across the profiles, at periods that range typically from 6.0 to 6.9 days, which are probably mostly driven by rotation, given the close match to the stellar rotation period. In Fig. 11, we show the Hα periodogram, where we can see a 6.5-day-period across the red wing, while the blue wing shows periodicities at both 6.00 days and about 8.3 days (Fig. 11). We do not have an obvious interpretation of the longer period in the blue wing of Hα but it could, for example, be related to a diskwind coming from a region outside the co-rotation radius.

Fig. 11 Hα periodogram obtained with the HARPS, ESPaDOnS, and SMARTS observations. The power scale ranges from 0 (white) to a maximum value of 18.7 (black).

The Hα profiles show a clear red-emission shoulder from phases 0.05 to 0.35 and the redshifted absorption appears from phases 0.6 to 0.8. The Hβ profiles have a redshifted absorption component in a more extended phase range, going from 0.5 to 0.9. This is consistent with the formation of Hβ in the accretion column, closer to the star than the source of Hα. In that way, the hot spot is seen through the accretion column projected along our line of sight during a more extended phase range in Hβ than in Hα, creating the redshifted absorption component. The redshifted absorption components are seen in both Hα and Hβ at phases that follow the passage of the hotspot. This indicates that the major accretion column must be trailing the hotspot passage. This could happen if the stellar magnetic field were at least partially anchored to the disk at radii beyond the co-rotation radius. Donati et al. (2011) obtained a co-rotation radius of r_cor ≃ 7.7 R_⋆ and a magnetospheric radius of 7.2 R_⋆ with the ESPaDOnS observations. The maximum value of the redshifted absorption velocity is  ~260  km s^-1 in Hα and  ~300  km s^-1 in Hβ. Although V2129 Oph was shown to have an octupolar magnetic-field component that is stronger than the dipolar component at the surface of the star, near the co-rotation radius the dipolar component is much stronger (by a factor of about 25) than the octupolar (Donati et al. 2011) and accretion is then expected to occur mostly through dipole field lines. Assuming that the system is seen at i = 60°, the hot spot is located at lat = 80°, and the accreting material free-falls along a dipolar accretion column, we obtain a truncation radius of 7.8 R _⋆ . This value is consistent with the possibility that part of the stellar magnetic field may interact with the disk beyond the co-rotation radius, creating a trailing accretion column in the 2009 observations.

The Hα and Hβ lines showed an increase in both intensity and width one day before (JD = 4998) the He i line, the mass accretion rate, and the veiling all increased significantly (Figs. 3, 4, 1). If Hα and Hβ are produced in a more extended region than He i, they could be affected by an increase in the mass accretion rate before the hot spot (veiling) and the He i line, which trace regions close to the accretion shock. Interestingly, the day Hα and Hβ increased, the He i line, which normally has only a slightly redshifted (v_peak ~ 4  km s^-1) narrow (FWHM ~ 20  km s^-1) component in V2129 Oph, showed a small broad component (BC), as can be seen in the two panels with phases 0.36 and JD = 4998 of Fig. 3. We decomposed these He i profiles using two Gaussians, and obtained BCs centered at  ~−20  km s^-1 and  ~−17  km s^-1 that have FWHM of 71  km s^-1 and 69  km s^-1, respectively. The BC is generally associated with hot winds (Beristain et al. 2001; Edwards 2003), hence this might indicate that an accretion burst has affected both the diskwind and the Balmer lines about one day before reaching the stellar surface. There was another burst episode in our data at JD = 4993, where the He i lines are again very strong. Unfortunately, we do not have, however, spectra taken one day before, to verify whether Hα and Hβ had also increased in advance of He i.

We have discussed the many characteristics of the Balmer emission lines, and we now attempt to reproduce them with theoretical line profiles, based on a 3D MHD circumstellar structure calculated for V2129 Oph.

4.3. Theoretical line profiles

Theoretical emission-line profiles of CTTSs have often been computed in the literature and compared to the observations to determine the main characteristics of the circumstellar accretion flow of these objects (Muzerolle et al. 1998, 2001; Kurosawa et al. 2006). However, the steady-state axisymmetric dipolar models that are more commonly used to describe the stellar magnetosphere cannot reproduce the variability observed in the emission profiles of CTTSs. These are mainly driven by the rotational modulation of non-axisymmetric multipolar fields and non-steady accretion processes, which are not straightforward to model. Progress has been made in the determination of the surface magnetic fields of CTTSs (Donati et al. 2007; Hussain et al. 2009; Gregory & Donati 2011) and the development of 3D MHD simulations that include non-axisymmetric multipolar components of the stellar magnetic field (Long et al. 2007, 2008; Romanova et al. 2008). It is now possible to compute emission line profiles from the funnel flow structure determined by 3D MHD simulations that are based on observed magnetic field configurations of CTTSs (Kurosawa et al. 2008, 2011).

The observed Hα and Hβ profiles of V2129 Oph cannot be easily decomposed with the simple Gaussian components typically used to represent the main accretion flow, the wind, and the accretion shock, in a similar way to the analyses of other stars, such as AA Tau (Bouvier et al. 2007). To understand the observed emission-line variabilities, we therefore decided to try to reproduce the Balmer line profiles with radiative transfer calculations based on the 3D MHD structure of the accretion flow. The Hα and Hβ profiles were calculated with the stellar parameters of V2129 Oph and two different magnetic field configurations: a pure dipolar field and the dipole+octupole configuration proposed by Donati et al. (2011) for the 2009 observations.

The 3D MHD code and model used here were developed and described earlier in Koldoba et al. (2002) and Romanova et al. (2003, 2004, 2008). The code has been modified by Long et al. (2007, 2008) to incorporate the quadrupole and octupole components, and 3D modeling of accretion onto BP Tau and V2129 Oph, with close to realistic fields, has been successfully performed (Long et al. 2011; Romanova et al. 2011). Here, we computed a special set of simulations for V2129 Oph, where the stellar magnetic field was approximated with a superposition of dipole and octupole fields, and the field values correspond to the more recent spectropolarimetric observations of this star (Donati et al. 2011). The parameter values used in the simulations are listed in Table 4. The inclinations of the dipole and the octupole with respect to the rotation axis were 15° and 25°, respectively, and the octupole is about 0.1, in rotational phase, ahead of the dipole, as proposed by Donati et al. (2011). The mass accretion rate used in the model, 6.3 × 10^-10 M_⊙ yr^-1, is slightly lower than the lowest values measured for the system (9 × 10^-10 M_⊙ yr^-1). This low accretion rate was necessary to generate velocity profiles that were compatible with the observed ones, since increasing the mass accretion rate would bring the disk closer to the star and produce narrower emission-line profiles. The simulations were done in dimentionless form, and the current simulations are applicable if we increase the dipole component of the field and accordingly increase the mass accretion rate. For example, the MHD simulation run would correspond to the observed mass accretion rate, if the dipole component of the field were increased by a factor of 1.2, such that the dipole field were 1.08 kG instead of 0.9 kG.

Fig. 12 Density slice in the xz plane for the 3D MHD simulations with dipole-only (top) and dipole+octupole (bottom) magnetic fields. The parameters used in the simulations are listed in Table 4. The length scales are in units of the stellar radius (R∗) and the intensity is shown on a logarithmic scale with an arbitrary unit.

The three-dimensional simulations show that, in a pure dipole case, matter flows in two ordered funnel streams towards the closest magnetic pole, as can be seen in Fig. 12 (top panel), which shows a slice of the density distribution in the xz plane and selected magnetic field lines. The disk is truncated by the magnetosphere at radius r_t ≈ 6−7 R_⋆, and matter hits the star below the magnetic pole to form two bean-shaped hot spots. In the dipole+octupole case (Fig. 12, bottom panel), the dipole component strongly dominates the field at large distances from the star, determining the truncation of the disk at similar distances as in a pure dipole case. However, close to the star, the octupole component dominates the flow and the octupolar field redirects the funnel to higher latitudes, where matter hits the star closer to the octupolar magnetic pole (see also Romanova et al. 2011). This has an effect in the hot spot location. In a pure dipole case, the hot spot forms at some distance (about 20°) from the magnetic pole. Hence, the hot spot and the cold spot, which is associated with the magnetic pole, should not coincide in phase. However, in the dipole+octupole case the hot spot is located much closer to the magnetic pole, because the octupolar component redirects the initial funnel flow to higher latitudes, which agrees with the observational results presented in Sect. 4.1 and in Donati et al. (2011).

Next, we calculated spectral lines from the MHD modeled flows. We followed the same approach described in Kurosawa et al. (2008) to obtain emission line profiles from the density, velocity, and temperature structures of the MHD simulations. These distributions were mapped on to the radiative transfer grid of the TORUS code (Harries 2000; Kurosawa et al. 2006, 2011), and the corresponding line source functions were calculated. The theoretical profiles of Hα and Hβ were computed using the Sobolev escape probability method, as described in Klein & Castor (1978) and Hartmann et al. (1994).

Our initial 3D MHD calculations were done with an adiabatic index of γ = 5/3 and, although the Hα and Hβ line profiles computed with the dipole-only magnetosphere were quite similar to the observations, the profiles obtained with the dipole+octupole magnetic field configuration were much narrower than the observed ones. The 3D simulations do not include radiative cooling, but only adiabatic cooling and heating. This leads to funnel temperatures that are too high, relative to observational values, and generates high pressure and pressure gradients in the funnel, which slow down the accretion. These effects were more pronounced in the presence of the octupole component, which, with its many loops close to the stellar surface (Fig. 12, bottom panel), acted as a wall against accretion. We then adjusted the adiabatic index of the gas in the equation of state to a lower value (γ = 1.2 instead or 5/3) to mimic a cooling effect. This helped to decrease the thermal pressure of the gas near the surface, and to produce a faster flow near the star. With the above adjustments, the flow geometry, on a large scale, for the dipole+octupole case (Figs. 15 and 16, top panels) is somewhat similar to that of the dipole-only case (Figs. 13 and 14, top panels).

Fig. 13 Hα model intensity maps (top panels) and the corresponding profiles (bottom panels, solid lines) computed at rotational phases of 0.0, 0.25, 0.5, and 0.75 (from left to right) for the 3D MHD simulation with a dipole-only magnetic field (see Fig. 12). The intensity maps are shown in the plane perpendicular to our line-of-sight. The parameters used in the simulations are listed in Table 4. The intensity is shown on a logarithmic scale with an arbitrary unit. The model profiles (solid lines) are compared with the observed mean Hα profiles (dotted lines) obtained with HARPS that are presented in Fig. 10. The profiles are normalized to the continuum level.

Fig. 14 Same as in Fig. 13, but for Hβ.

In the bottom panels of Figs. 13–16, we show the mean Hα and Hβ observed profiles at selected phases, overplotted on the mean theoretical line profiles at the same phases. We can see that the Balmer line profiles obtained with the models resemble overall the observed profiles for both Hα and Hβ and for both the dipole-only and dipole+octupole models. In Figs. 17 and 18, we show the theoretical and observed line EWs, the full width at half maxima (FWHM), and the width on the blue side at 10% peak intensity for both Hα and Hβ. We only use the blue side width at 10% peak intensity in our comparison, since the red side width at 10% peak intensity is strongly affected by the redshifted absorption component that is not very accurately modeled (see discussion later in this section).

Fig. 15 Hα model intensity maps (top panels) and the corresponding profiles (bottom panels, solid lines) computed at rotational phases 0.0, 0.25, 0.5, and 0.75 (from left to right) for the 3D MHD simulation with a dipole+octupole magnetic field (see Fig. 12). The intensity maps are shown in the plane perpendicular to our line-of-sight. The parameters used in the simulations are listed in Table 4. The intensity is shown on a logarithmic scale with an arbitrary unit. The model profiles (solid lines) are compared with the observed mean Hα profiles (dotted lines) obtained with HARPS and shown in Fig. 10. The profiles are normalized to the continuum level.

Fig. 16 Same as in Fig. 15, but for Hβ.

In Figs. 17 and 18 (top panels), we see that the Hα and Hβ theoretical EWs are of the same order of magnitude as the observed ones. However, the variability predicted by the models for the EW of both lines is quite different from the observations. On the other hand, the theoretical Hβ FWHM is very similar to the observed values, both qualitatively and quantitatively, while the theoretical FWHM of Hα are a bit overestimated when compared to the observations. The Hα and Hβ theoretical emission-line widths vary during the rotational cycle, being wider at phases 0.0 and 0.25 than at phases 0.5 and 0.75, as observed. However, as seen in the bottom panels of Figs. 17 and 18, although the qualitative behavior of the line width is similar between models and observations, the model profiles are much narrower than the observed ones. This is also clearly seen in the bottom panels of Figs. 13–16. In these figures, we can see that the intensity and shape of the theoretical and observed emission lines closely resemble each other, especially around phases 0.0 and 0.25. At phases 0.5 and 0.75, the match is not so good, the line strengths of the models being slightly weaker and stronger than those of the observations at phases 0.5 and 0.75, respectively. The differences in the line strengths at phases 0.5 and 0.75 may be attributed to differences between the flow geometry of the MHD models and the real system. The extension of the red wing is similar in models and observations and both the dipole-only and dipole+octupole models display clear redshifted absorption going below the continuum and extending up to 300–350  km s^-1 in Hβ at phase 0.5, in agreement with the observations. In the dipole-only models, we also see some hint of redshifted absorption in Hβ at phase 0.75, as observed, while no absorption is seen at these phases in the dipole+octupole models. Most of the line variability and the appearance of redshifted absorption in the Balmer lines can therefore be attributed to changes in the observer’s viewing angle of the system due to a combination of the rotational motion of the star and its inclined and non-aligned magnetosphere.

Fig. 17 Hα line EW (top), full width at half maxima (middle), and blue side width at 10% of the peak intensity (bottom) as a function of rotational phase. HARPS and ESPaDONs data are shown as filled and open dots, respectively. The values measured in the profiles of the dipole-only (solid line) and dipole+octupole (dashed line) models are overplotted.

Fig. 18 Same as in Fig. 17, but for Hβ.

A noticeable difference between models and observations is the strength and extension of the blue wings of the Balmer line profiles (especially in Hβ), as can be seen in the bottom panels of Figs. 13–16. The emission in the blue wings extends to  ~300   km   s^-1 in the observations, but it extends only to  ~150   km   s^-1 in most of the Hβ theoretical profiles. An exception is the dipole+octupole model at a phase of around 0.25, which shows a weak extension of the blue wing emission up to  ~300   km   s^-1. Possible causes of this discrepancy are (1) the uncertainty in the flow geometry; (2) the inclination angle; or (3) the presence of non-negligible wind emission. Although the observations do not show any clear sign of an outflow/wind, such as a blueshifted absorption component, the wind may contribute to the emission. Unfortunately, the wind is not included in our current MHD models.

Another difference is that the observations show clear redshifted absorption in both Hα and Hβ, while in the present model, there is clear redshifted absorption only in Hβ. The redshifted absorption in Hβ also appears at lower velocities in the simulations, starting at 50  km s^-1, than in the observations, starting at 100  km s^-1, which reflects the inflow speed of the gas in the adopted MHD solution. Another small difference between the models and observations is that, in the MHD simulation, the major funnel flow is centered around the hot spot, instead of mostly trailing it, as seen in the observations. In the dipole-only simulations, the Hβ redshifted absorption appears from phases 0.40 to 0.68. In the dipole+octupole simulations, it appears from phases 0.44 to 0.64, almost centered at phase 0.5, when the hot spot faces the observer, while in the observations it appears from phases 0.45 to 0.9 (see Sect. 4.2). This difference is probably caused by the magnetosphere being inside the co-rotation radius (R_c = 7.7 R_⋆), going from about 5.5 R_⋆ to 7.5 R_⋆ in the 3D MHD simulations (Fig. 12, top). The inner radius in the simulations was estimated from the location of the edge of the funnel flow, and the outer radius from the location of the largest closed magnetic field line. The outer magnetosphere radius is actually quite close to the magnetospheric radius derived by Donati et al. (2011) (7.2 R_⋆). However, since it is smaller than the co-rotation radius, it is inconsistent with a trailing funnel flow and it might instead represent a funnel flow that leads rather than trails the hot spot in the simulations. A way to increase the phase range where redshifted absorption appears in the theoretical profiles would be to either decrease the tilt of the dipole component or increase the mass accretion rate of the system. Our MHD simulations usually show that, decreasing the tilt of the dipole component, the magnetic funnel becomes wider in the azimuthal direction and the width of the funnel curtain also depends on the density level (Romanova et al. 2003). At lower density levels, the magnetic curtain covers a wider range of azimuthal angles, even at the present dipole tilt. Hence, at a somewhat higher mass accretion rate, we may have a wider curtain at the density level required to obtain the redshifted absorption in Hβ. This would make the redshifted absorption appear in a more extended phase range.

Fig. 19 Hα autocorrelation matrices. The gray scale represents the linear correlation coefficient (r). White corresponds to a positive correlation (r = 1) and black to an anti-correlation (r =  −1). Top left: observed data. Top right: dipole-only model with γ = 5/3. Bottom left: dipole-only model with γ = 1.2. Bottom right: dipole+octupole model with γ = 1.2.

The observed line profiles may come from different emitting/absorbing regions, which is specially true for Hα, and it is interesting to see how these different regions are correlated. To evaluate the correlation variability across the line profile, we computed autocorrelation matrices for Hα with the observed and theoretical line profiles. The results are presented in Fig. 19. Donati et al. (2007) found an anti-correlation between the red wing (v > 100  km s^-1) and the central emission peak in their ESPaDOnS data taken over seven days, covering one rotational cycle. They attributed the anti-correlation to the contributions at different velocities to the Hα profile from distinct physical regions (accretion pre-shock, winds, shock region). With the present data, we are unable to see such an anti-correlation between the red wing and the emission peak, and the two regions actually do not seem to be strongly correlated, which agrees with the matrices calculated with the theoretical profiles. This could indeed be due to the redshifted absorption component coming from the pre-shock region of the accretion funnel, while the central part of the emission profile might have contributions from either a much larger volume of the accretion column or a wind. They could then represent independent variations. The differences between the two observational sets (Donati et al. 2007 and ours) could be due to the larger number of rotational cycles covered with the present data set that would tend to smooth out any peculiarities in a single rotational cycle. The observed Hα autocorrelation matrix bears some resemblance to the dipole-only autocorrelation matrix (especially the γ = 5/3 model). Although the degree of correlation is not exactly the same across the line in both matrices, most of the correlated/anti-correlated regions appear in both the dipole-only and the observed matrices. Alternatively, the dipole+octupole matrix contains many correlated/anti-correlated regions that are not present in the matrix obtained with the observed data.

None of our theoretical models have included a wind component, which is expected to be present in these systems, even if it does not always contribute significantly to the line profile (Kurosawa et al. 2006; Lima et al. 2010). There are some blueshifted features, such as a shoulder, in the observed Hα profiles that could be due either to an absorption, coming from a wind, or to an emission excess, from different projections of the accretion funnel. To verify the importance of a diskwind to the observed Hα line, we calculated theoretical Hα profiles with the hybrid MHD model of Lima et al. (2010), which includes both magnetospheric and diskwind components. In this model, the magnetosphere is an axisymmetric dipole and we aim to reproduce the average observed profile of V2129 Oph, using the parameters presented in Table 4. The magnetosphere extends from 7.2 R_⋆ to 7.8 R_⋆, where we kept the maximum temperature of the magnetosphere around T_mag,max = 9000 K and used the photospheric temperature of T_phot = 4500 K. Slight changes in T_mag,max, as shown in Fig. 20, led to variations in the line intensity and details in the line profile, such as a more pronounced redshifted shoulder. The maximum wind temperature used in the models was T_wind,max = 9000 K and the temperature of the hot spot was varied from T_ring = 5500 K to 8000 K, without there being much change to the line profile. The diskwind starts at 7.81 R_⋆, immediately outside the outer magnetospheric radius and its extension was varied from 30 R_⋆ to 8 R_⋆, while the mass outflow rate remained constant (Ṁ_out = 0.1Ṁ_acc). We therefore increased the wind density as we decreased the outer diskwind radius, although neither significant blueshifted absorption nor emission were ever observed in the theoretical Hα profiles. This agrees with the results discussed in Lima et al. (2010), where it was shown that the wind contribution to the line profile is negligible for Ṁ_acc < 10^-9 M_⊙ yr^-1 owing to very low wind densities, which is the case in the present calculations. The mean observed Hα profile is well-reproduced by the model (Fig. 20, solid lines), and our diskwind component has no visible influence on the line profile at the low mass-accretion rate of V2129 Oph (dashed lines). With the parameters of Table 4 we were unable to reproduce the redshifted absorption component in Hα commonly observed in V2129 Oph. We then varied the system inclination to mimic a non-aligned magnetospheric and rotation axis. With a lower inclination, on the order of 45−50°, a redshifted absorption component was found to be visible, but the profile intensity was lower than usually observed.

5. Conclusions

The magnetospheric accretion scenario predicts that the stellar magnetic field will interact with the disk and disrupt it when the magnetic pressure overcomes the ram pressure due to the accretion process. Accretion columns will then be formed and a hot spot/ring will result at the stellar surface where the accreting gas hits the star. In the past few years, we have seen evidence from observations and numerical simulations that the star-disk interaction is dynamic and mediated by a stellar magnetic field that may vary in time and be composed of different multipoles that have to be taken into account in the magnetospheric accretion paradigm. The magnetic field moments are also generally found to be inclined with respect to the rotation axis, creating a non-axisymmetric circumstellar environment, whose observed characteristics change on a timescale of days, owing to the system rotation.

Fig. 20 Hα profiles calculated with the dipolar and axisymmetric hybrid model that includes both magnetospheric accretion and diskwind (solid lines). The dashed lines show only the diskwind contribution to each profile. We used the parameters presented in Table 4, a maximum wind temperature of Twind,max = 9000 K and a photosphetric temperatures of Tphot = 4500 K. We varied the maximum temperature of the magnetosphere from Tmag,max = 9000 K (left) to Tmag,max = 8950 K (right).

We have shown in this paper that the magnetic field configuration proposed by Donati et al. (2011) for the classical T Tauri star V2129 Oph, which includes dipole and octupole fields inclined with respect to the rotation axis, reproduces quite well the photometric and spectroscopic variability observed over several rotational cycles of the system.

We have used the magnetic field configuration obtained from observations in the calculations of 3D MHD magnetospheric simulations, whose density, velocity, and scaled temperature structures were mapped on to a radiative transfer grid to obtain line source functions from which we calculated the theoretical line profiles. The Hα and Hβ observed emission profiles vary due to rotation modulation and this variability is qualitatively well reproduced by the computed theoretical models. The observed profiles are however wider than the theoretical ones (especially for Hβ), containing an extended emission in the blue wing and broader-than-predicted full widths at half maximum. This indicates that additional modeling is needed to fully describe the accretion flows around CTTSs.

Spectral lines calculated for MHD model accretion flows with pure dipole and dipole+octupole fields are similar, since in both cases they form in funnel streams that are dominated by the dipole component at the truncation radius. However, the correlated variability across the emission line profiles has been most accurately reproduced by a dipole-only model. We also calculated emission line profiles with a dipolar and axisymmetric model that included a magnetosphere and a diskwind, but our diskwind model did not provide any significant contribution to the Hα line profile of V2129 Oph.

The hot and cold spots of V2129 Oph were found to be almost coincident in phase and at high latitudes but are expected to be located at different levels, the hot spot at the chromosphere and the cold spot at the photosphere. The inferred cold and hot spots are able to explain the radial velocity variations in the photospheric and He i (5876 Å) lines. The phase coincidence of the hot and cold spots is more accurately explained by the dipole+octupole models than the dipole-only model, since the octupolar component redirects the funnel flow towards the magnetic pole, while in the dipole-only case the magnetic pole and hot spots are separate.

The mass accretion rate of the system is generally around 1.5 × 10^-9 M_⊙ yr^-1, but accretion bursts can occur, as observed, on timescales of days, when the mass accretion rate is higher by as much as three times its quiescent value.

The spectroscopic and photometric variabilities observed in V2129 Oph are thus consistent with the general predictions of complex magnetospheric accretion models with non-axisymmetric, multipolar fields.

Acknowledgments

We gratefully acknowledge Christophe Lovis who obtained part of the HARPS spectra. S.H.P.A. acknowledges support from CAPES, CNPq and Fapemig.

References

All Tables

All Figures

Fig. 1 V band photometry corrected from extinction (panel a)). SMARTS data are shown as squares, Palermo data as stars, and Maidanak data as triangles. Photospheric radial velocities (panel b)), veiling (panel c)), and He i (5876 Å) radial velocity (panel d)) phase variations. In panels b)–d), filled symbols are from HARPS data and open ones from ESPaDOnS. The error bars are shown for each observation in panel a), and are smaller than the symbol sizes in panel b). We show in panels c and d only the mean error bars to avoid crowding the figures. Different colors represent different rotational cycles identified as the integer part of the observational phase in Table 1: black – cycle 0, green – cycle 1, red – cycle 2, aqua – cycle 3, dark blue – cycle 4, pink – cycle 5.

In the text

	Fig. 2 Hα profiles obtained with the HARPS spectrograph. The profiles are normalized to the continuum level and the vertical dotted lines indicate the stellar rest velocity. Phases and JDs are given for each spectrum in the panels. The color code is the same as in Fig. 1.
In the text

	Fig. 3 He i 5876 Å profiles obtained with the HARPS spectrograph. The profiles are normalized to the continuum level and the vertical dotted lines indicate the stellar rest velocity. Phases and JDs are given for each spectrum in the panels. The color code is the same as in Fig. 1.
In the text

Fig. 4 Rotational variation in the He i EW and the mass accretion rate. The mass accretion rate errors come from the relation between the He i (5876 Å) luminosity and the accretion luminosity by Fang et al. (2009). Filled symbols are from HARPS data and open ones from ESPaDOnS. Different colors represent different rotational cycles identified as the integer part of the observational phase in Table 1: black – cycle 0, green – cycle 1, red – cycle 2, aqua – cycle 3, dark blue – cycle 4, pink – cycle 5.

In the text

Fig. 5 Best-fit hot spot model for the HARPS veiling variations (black lines). The red lines represent the veiling generated by the optimal hot spot fits to the He i radial velocity variations (Fig. 6). Models with i = 60° and i = 45° are shown as solid and dashed lines, respectively . The mean error bar in the observed veiling is shown on the left and the color code of the points is the same as in Fig. 1.

In the text

	Fig. 6 HARPS (filled symbols) and ESPaDOnS (open symbols) He i radial velocity measurements, with respective error bars. The red curves represent the best-fit hot-spot model for the He i radial velocity data obtained with i = 60° (solid line) and i = 45° (dashed line). The black solid line corresponds to the best fit to the veiling variations with i = 60°, presented in Fig. 5.
In the text

	Fig. 7 χ2 contours of the cold spot models with i = 45°. Each contour level corresponds to an increase of χ2 minimum. The red lines correspond to a 2χ2 minimum and represent the error in each parameter quoted in the text.
In the text

	Fig. 8 Best-fit cold-spot model with i = 45° to the photometric observations. The parameters of the cold-spot model used to generate the above light curves correspond to a χ2 minimum (Tspot = 4200 K, spot latitude of 60°, spot phase of 0.57, and spot filling factor of 40%).
In the text

	Fig. 9 Radial velocity variations induced by a cold spot model obtained in 10 000 Monte Carlo runs varying the parameter space of possible solutions (small black dots). The triangles represent the median values of the theoretical solutions and the HARPS photospheric radial velocities are shown as large filled circles.
In the text

	Fig. 10 Hα (top) and Hβ (bottom) profiles obtained with the HARPS (solid), ESPaDOnS (dashed), and SMARTS (dash-dotted, Hα only) spectrographs. The spectra correspond to a mean profile around the phases given in the panels. The profiles are normalized to the continuum level and the vertical dotted lines indicate the stellar rest velocity.
In the text

	Fig. 11 Hα periodogram obtained with the HARPS, ESPaDOnS, and SMARTS observations. The power scale ranges from 0 (white) to a maximum value of 18.7 (black).
In the text

	Fig. 12 Density slice in the xz plane for the 3D MHD simulations with dipole-only (top) and dipole+octupole (bottom) magnetic fields. The parameters used in the simulations are listed in Table 4. The length scales are in units of the stellar radius (R∗) and the intensity is shown on a logarithmic scale with an arbitrary unit.
In the text

Fig. 13 Hα model intensity maps (top panels) and the corresponding profiles (bottom panels, solid lines) computed at rotational phases of 0.0, 0.25, 0.5, and 0.75 (from left to right) for the 3D MHD simulation with a dipole-only magnetic field (see Fig. 12). The intensity maps are shown in the plane perpendicular to our line-of-sight. The parameters used in the simulations are listed in Table 4. The intensity is shown on a logarithmic scale with an arbitrary unit. The model profiles (solid lines) are compared with the observed mean Hα profiles (dotted lines) obtained with HARPS that are presented in Fig. 10. The profiles are normalized to the continuum level.

In the text

	Fig. 14 Same as in Fig. 13, but for Hβ.
In the text

Fig. 15 Hα model intensity maps (top panels) and the corresponding profiles (bottom panels, solid lines) computed at rotational phases 0.0, 0.25, 0.5, and 0.75 (from left to right) for the 3D MHD simulation with a dipole+octupole magnetic field (see Fig. 12). The intensity maps are shown in the plane perpendicular to our line-of-sight. The parameters used in the simulations are listed in Table 4. The intensity is shown on a logarithmic scale with an arbitrary unit. The model profiles (solid lines) are compared with the observed mean Hα profiles (dotted lines) obtained with HARPS and shown in Fig. 10. The profiles are normalized to the continuum level.

In the text

	Fig. 16 Same as in Fig. 15, but for Hβ.
In the text

	Fig. 17 Hα line EW (top), full width at half maxima (middle), and blue side width at 10% of the peak intensity (bottom) as a function of rotational phase. HARPS and ESPaDONs data are shown as filled and open dots, respectively. The values measured in the profiles of the dipole-only (solid line) and dipole+octupole (dashed line) models are overplotted.
In the text

	Fig. 18 Same as in Fig. 17, but for Hβ.
In the text

	Fig. 19 Hα autocorrelation matrices. The gray scale represents the linear correlation coefficient (r). White corresponds to a positive correlation (r = 1) and black to an anti-correlation (r =  −1). Top left: observed data. Top right: dipole-only model with γ = 5/3. Bottom left: dipole-only model with γ = 1.2. Bottom right: dipole+octupole model with γ = 1.2.
In the text

Fig. 20 Hα profiles calculated with the dipolar and axisymmetric hybrid model that includes both magnetospheric accretion and diskwind (solid lines). The dashed lines show only the diskwind contribution to each profile. We used the parameters presented in Table 4, a maximum wind temperature of Twind,max = 9000 K and a photosphetric temperatures of Tphot = 4500 K. We varied the maximum temperature of the magnetosphere from Tmag,max = 9000 K (left) to Tmag,max = 8950 K (right).

In the text