© The Authors 2023

1. Introduction

The early evolution of low-mass stars (M_* < 2 M_⊙) during the classical T Tauri (CTT) phase depends on the interaction between the star and its accretion disc, on a distance of a few stellar radii. At the truncation radius, matter from the disc surface is channelled onto the stellar surface following the magnetic field lines and forming an accretion funnel or column (Ghosh et al. 1977; Zanni & Ferreira 2009; Romanova & Owocki 2016; Pantolmos et al. 2020). The star-disc interaction is responsible for accretion and ejection phenomena that have a strong impact on spectral lines formed in the close vicinity of the star’s surface.

Ghosh et al. (1977) developed an analytical model of magnetospheric accretion around a rotating neutron star with a dipolar magnetic field. Hartmann et al. (1994) applied this magnetospheric accretion model to the formation of emission lines in the spectrum of T Tauri stars. This fundamental paper set the general theoretical framework for the density and temperature distributions in aligned axisymmetric magnetospheres. The coupling between this representation of magnetospheric accretion in T Tauri systems with radiative transfer calculations has provided a crucial tool to interpret spectroscopic, photometric, and interferometric observations. The sensitivity of hydrogen lines to the parameters of the magnetospheric models was studied in detail by Muzerolle et al. (2001), improving the earlier calculations by Hartmann et al. (1994).

Near-infrared observations of the Brackett γ (Brγ) line with the Very Large Telescope Interferometer (VLTI) GRAVITY instrument (GRAVITY Collaboration 2017) also probe the inner part of the star-disc interaction region (GRAVITY Collaboration 2020; Bouvier et al. 2020a). However, it is still difficult to associate the characteristic sizes derived from interferometry with the actual size of the magnetospheric accretion region, which is a key parameter in our understanding of the star-disc interaction.

In this work, we aim to study the formation of the Brγ line and compute its spectroscopic and interferometric signatures for non-axisymmetric models of the inner star-disc interaction region, akin to state-of-art magnetohydrodynamics (MHD) simulations (Romanova & Owocki 2016). In particular, we want to clarify the meaning of the sizes inferred through near-infrared interferometric observations and how they compare with the overall size of the magnetospheric accretion region.

In Sects. 2 and 3, we describe the model used to compute the line formation in accreting magnetospheres. We discuss spectroscopic and interferometric signatures in Sects. 4 and 5, respectively.

2. Radiative transfer framework

We used the MCFOST code¹ (Pinte et al. 2006, 2009; Tessore et al. 2021) to compute emergent line fluxes from multi-dimensional models of magnetospheres for a 20-level hydrogen atom. The atomic model, with 19 bound levels and the ground state of H II, consists of 171 bound-bound transitions (atomic lines) and 19 bound-free transitions (continua). We focus here on the Brγ line at 2.1661 μm, although the Balmer lines Hα and Hβ and the Paschen β line (Paβ) are modelled as well. These specific hydrogen lines are commonly used to characterise accretion and ejection phenomena in young systems (Folha & Emerson 2001; Alencar et al. 2012; Bouvier et al. 2020a; Pouilly et al. 2020; Sousa et al. 2021). The method to solve for the non-local thermodynamic equilibrium (non-LTE) populations of hydrogen and the microphysics are the same as in Tessore et al. (2021). The updated version of the code we use now simultaneously solves the charge equation and the statistical equilibrium equations, which has been proven to increase the convergence in chromospheric conditions (Leenaarts et al. 2007). We tested our code for different magnetospheric models taken as benchmarks in Muzerolle et al. (2001) and Kurosawa et al. (2006). The results of this comparison are presented and discussed in Appendix A.

3. Magnetospheric accretion model

Matter from the circumstellar disc is channelled onto the stellar surface along the dipolar magnetic field lines. The stellar magnetic field truncates the disc at a distance R_t from the star, the truncation radius. In practice, the interaction between the stellar magnetic field and the disc takes place over a small region between R_t and R_t + δr. Both R_t and δr are used to define the size of the disc region magnetically connected to the star. As the gas approaches the stellar surface, it decelerates in a shock and is heated at coronal temperatures. Theoretical models of accretion shocks by Calvet & Gullbring (1998) show that the optically thin emission of the pre-shock and post-shock dominates below the Balmer jump and that the optically thick emission of the heated photosphere contributes to the total continuum emission at larger wavelengths. In the following, we only consider the contribution of the heated photosphere to the shock radiation. The shock² temperature was computed from the energy of the gas infalling onto the stellar surface following the prescription of Romanova et al. (2004) unless stated otherwise. This approach assumes energy conservation and that the shock radiates as a black body, meaning that its temperature is determined by the specific kinetic energy and enthalpy of the gas deposited at the stellar surface. The shock temperature hence derived is of the order of 4500 K–6000 K.

3.1. The stellar surface

The stellar surface is considered as the inner boundary of the model and emits as a blackbody whose temperature is determined by the stellar parameters. Throughout the paper, the stellar parameters are T_* = 4000 K, M_* = 0.5 M_⊙, and R_* = 2 R_⊙. We set the distance to the star at 140 pc, which is typical of the nearest star-forming regions such as Upper Scorpius (≈146 pc Galli et al. 2018a) or Taurus (≈130 pc Galli et al. 2018b).

3.2. Geometry of the accretion funnels

We considered 3D non-axisymmetric models of the magnetospheric accretion region. These models are parametrised by the same set of parameters as the axisymmetric magnetospheric model of Hartmann et al. (1994; see also Muzerolle et al. 1998, 2001; Kurosawa et al. 2006, 2011; Lima et al. 2010; Dmitriev et al. 2019). The density and the velocity fields of the accretion columns are fully described with a set of independent parameters: the mass accretion rate Ṁ, the rotation period P_rot, R_t, and δ_r.

For our study, the value of Ṁ, R_t and δ_r, and of the temperature of the magnetosphere are fixed. The impact of these parameters on the line formation has been discussed thoroughly in Muzerolle et al. (1998, 2001, see also Appendix A). The line’s response to the mass accretion rate and to the temperature is an essential proxy for understanding the physics of the star-disc interaction region. We use a mass accretion rate Ṁ = 10⁻⁸ M_⊙ yr⁻¹, a truncation radius R_t = 4R_*, and δ_r = 1R_*. The value of the rotation period is deduced from the maximum truncation radius (R_t + δ_r), imposing that stable accretion occurs at 90% of the corotation radius, consistent with the work of Blinova et al. (2016). The rotation period is therefore fixed at P_rot = 6 days, corresponding to slowly rotating T Tauri stars (see Herbst et al. 2007 and Bouvier et al. 2014 for a review). The rotational velocity for that period is thus of the order of 80 km s⁻¹ at the outer edge of the magnetosphere.

When the magnetic field axis (μ) is misaligned with respect to the rotational axis (Ω) of the star, the geometry of the accretion flow changes dramatically. The equations for the magnetic field components of a non-axisymmetric dipole, that is with a non-zero obliquity, are provided in Mahdavi & Kenyon (1998). The parameter β_ma describes the angle between the dipole moment and the star’s rotational axis, which is the magnetic obliquity.

We approximate the density, ρ, along the non-axisymmetric magnetic field lines with

$$\begin{array}{c}\hfill \rho =\alpha \frac{B}{v}=\alpha B\frac{{\rho }_{\mathrm{axi}}}{B_{\mathrm{axi}}},\end{array}$$(1)

where α is a constant along a given field line and B is the analytic misaligned dipolar field. v, ρ_axi, and B_axi denote the velocity field, density, and dipolar magnetic field, respectively, and they are taken from the axisymmetric model of Hartmann et al. (1994). In other words, the 3D density structure is computed from Eq. (1) under the assumption that the infalling gas has a velocity field on the poloidal plane. The value of α is computed from the numerical integration³ of the mass flux over the shock area,

$$\begin{array}{c}\hfill {\displaystyle \dot{M}={\displaystyle \int }\rho \mathbf{v}·\mathrm{d}\mathbf{S},}\end{array}$$(2)

where dS is the surface element and v is the velocity field. In our model, the value of Ṁ is an input parameter and is held constant. Therefore, α is obtained to ensure consistency between Eqs. (1) and (2). We compute five models with an obliquity β_ma ranging from five to forty degrees in steps of ten degrees, which is representative of what has been measured for T Tauri stars with spectroscopy (McGinnis et al. 2020) and spectropolarimetry (Donati et al. 2008, 2010, 2013; Johnstone et al. 2014; Pouilly et al. 2020). For these non-axisymmetric models, the shortest field lines – defining the main accretion columns⁴ – carry most of the gas density. We remove the longest field lines – the secondary columns – from our modelling, as was done in Esau et al. (2014). This yields models with one crescent-shaped accretion spot per stellar hemisphere reminiscent of numerical simulations of misaligned dipoles (Romanova et al. 2003).

Figure 1 shows the density of a non-axisymmetric magnetosphere with an obliquity of 10°.

Fig. 1. Density distribution of non-axisymmetric model with an obliquity of 10°. The rotation axis of the star, Ω, is shown with a white arrow and the dipole axis, μ, is shown with a red arrow. The density is computed from Eqs. (1) and (2). The colour map scales with the density.

3.3. Temperature of the funnels

The temperature of the magnetospheric accretion region is not well constrained. The determination of the temperature by Martin (1996) from first principles was not able to reproduce the observations. A self-consistent calculation of the temperature of the magnetosphere is beyond the scope of this paper. Instead, we adopted the temperature profile of Hartmann et al. (1994), which has been extensively used in the past to model line fluxes from accreting T Tauri stars. The temperature is computed using a volumetric heating rate (∝r⁻³) and balancing the energy input with the radiative cooling rates of Hartmann et al. (1982). The exact balance between the heating and cooling mechanisms is unknown. Instead, the temperature profile is normalised to a free parameter, T_max, that sets the value of the maximum temperature in the funnel flow. In the following, we have set the temperature maximum to T_max = 8000 K.

4. Spectroscopic signatures

Thanks to the Doppler shift of the funnel flow, it is possible to reconstruct the origin of the emission line by looking at the brightness maps in various velocity channels. Figure 2 shows the contribution of the different parts of the magnetosphere to the total integrated Brγ line flux at a given velocity for an inclination of 30°, matching the model illustrated in Fig. 1. At that density and temperature, the continuum emission comes from the stellar surface (I_surf/I_mag > 100). Locally, the continuum emission from the shock is three times larger than the emission from the star. Overall, given the small covering area of the accretion shock (around 1%), the total continuum emission at the frequency of the Brγ line is dominated by the star’s radiation: F_shock/F_* = 3%.

Fig. 2. Origin of emission seen across the Brackett γ line. The contribution of individual images to the total line flux is indicated on the central image showing the line profile. The brightness maps are in units of the maximum emission. The emission of the stellar surface is saturated. Orange to red colours indicate the regions of maximum emission. The system is seen at an inclination of 30° and an rotational phase of ∼0.25, similarly to Fig. 1.

The low-velocity components (< 50 km s⁻¹) of the line form in the regions where the projected velocity along the line-of-sight is close to zero and near the disc.

The geometry of the non-axisymmetric model, defined in Sect. 3, is responsible for a rotational modulation of the integrated line flux. Many classical T Tauri stars show modulated photometric variability (e.g., Cody et al. 2014), and, more directly related to the magnetospheric region, many also show rotational modulation of the longitudinal component of the stellar magnetic field (e.g., Donati et al. 2020). Indeed, the periodic variability of optical and emission line profiles has been reported in various systems (e.g., Sousa et al. 2016; Alencar et al. 2018; Bouvier et al. 2020a), which indicates that the emission region is stable on a timescale of several rotation periods.

Figure 3 shows the variability of the Brγ line at different phases of rotation at an inclination of 60° for different obliquities. The origin of the rotational phase is defined such that at phase 0.5, the accretion shock is facing the observer. The red-shifted absorption seen for the Brγ line at phases 0.250, 0.47, and 0.69 results from a lower source function of the gas above the shock (see Appendix A). From observations, redshifted absorption in the Paβ and Brγ lines is seen in less than 34% and 20% of the line profiles, respectively (Folha & Emerson 2001). The inverse P Cygni profile disappears when the shock, or a significant fraction of it, is hidden on the opposite side of the star. The line, with either a double-peaked profile or a moderately redshifted absorption, is reminiscent of cases II and IV in Reipurth et al. (1996). While the profiles with redshifted absorption agree with observations, those that display an M shape are usually not observed in young stellar objects. This suggests that magnetospheric accretion is not the only contribution to the profile, which can also be impacted by various types of outflows (e.g., stellar, interface, and disc winds Lima et al. 2010; Kurosawa et al. 2011). The optically thick accretion disc is not included in our models. The effect of the disc emission and absorption on the spectroscopic and interferometric observables will be discussed in a subsequent paper.

Fig. 3. Brackett γ line variability along the rotational cycle. Each column corresponds to a specific rotational phase. At phase 0, the shock area is unseen on the stellar surface, while at a phase of 0.5, the shock is fully seen on the visible hemisphere. The colours correspond to different values of the obliquity. All fluxes are computed with an inclination of 60°.

We also observe a decrease in the line flux as the obliquity increases. Figure 4 shows the radius encompassing 90% of the total line flux, R₉₀, at an inclination of 60° for non-axisymmetric models with different obliquities for the Hα, Hβ, Paβ, and Brγ lines.

Fig. 4. Radius encompassing 90% of total flux (R90) for each line as a function of the obliquity, βma. Hydrogen lines are labelled with different colours.

As β_ma increases, the volume of the magnetospheric accretion region decreases because the arc length of the accreting field lines shortens. Therefore, the total flux, for all lines, decreases accordingly, independently of the viewing angle of the system. However, we also note a dependence of R₉₀ on the line. The value of R₉₀ represents the size of the emitting region on a given line, which is a function of density and temperature, and that of the viewing angle.

5. Interferometric signatures

In this section, we compute the size of the Brγ line-emitting region inferred from interferometric observations, and we compare it to model flux radii (see Sect. 4).

5.1. Interferometric observables

The interferometric observables are derived from the radiative transfer (RT) model using the ASPRO2⁵ software developed by the Jean-Marie Mariotti Center (JMMC). These observables represent what would be observed with GRAVITY in the near-infrared. We considered the configuration obtained with the Very Large Telescope (i.e. 4×8m telescopes), encompassing a range of baselines from 35 to 135 m. With a typical eight-hour night, we computed one observing point per hour for the six baselines of the VLTI to increase the Fourier sampling, namely uv coverage, which is crucial for the fitting part of our approach. As described in Bourgès & Duvert (2016), we derived the observables from the RT images (see Fig. 2) by computing the complex visibility in each spectral channel around the Brγ line and interpolating them to match GRAVITY’s spectral resolution (R = 4000). Specifically, we simulated a total of 37 spectral channels (from 2.161 to 2.171 μm with a step of 2.8 10⁻⁴ μm) for the six projected baselines repeated eight times. Within this range, 31 spectral channels are used to measure the K-band continuum and six channels sample the Brγ line emitting region.

Figure B.1 illustrates the resulting uv plane projected on-sky for a typical object observed at the VLTI with a declination of −34° (e.g., TW Hydrae). Figure 5 shows the interferometric observables along the rotational cycle for a model with an inclination of 60° and an obliquity of 10°. The two main observables are the modulus of the complex visibility – the visibility amplitude – and the differential phase – its argument – dispersed in wavelength. The phase is normalised to zero in the continuum. The visibility amplitude can then be used to estimate the object’s size, while the phase measures the photo-centre shifts between the line-emitting region and the continuum. The phase can only be used as a relative measurement (e.g., between the line and the continuum), the absolute phase being lost due to a combination of atmospheric and instrumental effects. We repeated the simulated observations and computed nine datasets over a rotational cycle sampled every 40 degrees (0.11 in phase). In this study, we are interested in the line’s emitting region only. Therefore, we used pure line quantities, instead of total visibilities and phases, to remove the contribution from the stellar surface (see Appendix B for the derivation of the pure line interferometric quantities).

Fig. 5. Synthetic interferometric measurements and modelling across the rotational phase of the system. Top: integrated images over the Brγ line. Green (uniform disc) and yellow (Gaussian disc) ellipses are the characteristic sizes measured with a GRAVITY-like instrument. Middle-top: pure line visibility amplitude observables associated with the corresponding models. The visibility variation (for a given baseline) as the uv plane rotates is the specific signature of an elongated object. Middle-bottom: pure phase visibility across the line profile for the six baselines of the VLTI. Colours encode the observing time. Bottom: velocity map of the radiative transfer model. The coloured dots represent the measurement of the photo-centre derived from the phase visibility in each available spectral channel (see Appendix B). In each figure, the magnetosphere and the stellar surface have been normalised independently for display purposes.

5.2. Physical characteristics and sizes

Once the interferometric observables are computed, we applied standard modelling methods to interpret the data (Berger 2003). Firstly, we averaged the visibility amplitude of the six spectral channels available within the Brγ line⁶. We used the average visibilities to recover the global size of the Brγ emitting region, where the different velocities probe specific parts of the moving material within the magnetosphere. Then, we fitted the averaged visibility amplitude using elongated Gaussian or uniform disc models. Such models are typically used in interferometry to estimate the system’s characteristic size and on-sky orientation. The source’s brightness distribution is defined by an elongation factor and a position angle, and its half-flux radius or radius in the case of a Gaussian disc or a uniform disc model, respectively. In the following, we adopt the definition of ‘radius’ for both models; this corresponds to the half-flux semi-major axis for the Gaussian model and the semi-major axis for the uniform disc model. The recovered sizes and orientations are represented in the top panel of Fig. 5. While neither model can fully account for the size of the magnetosphere, the uniform disc probes a larger area of the magnetosphere, while the Gaussian disc seems limited to the most luminous parts. We note that the fit of the visibility is equally good for both models and, thus, does not allow us to discriminate between the models from the synthetic visibilities only (middle-top, Fig. 5).

In order to quantify the physical meaning of the interferometric measurements, we compared the interferometric sizes with reference flux radii of the RT models. We set these radii to represent 50, 80, 90, and 99% of the total flux emitted by the magnetospheric accretion region. Figure 6 compares the sizes derived with interferometry to the characteristic radii of the RT models.

Fig. 6. Interferometric radii as function of rotational phase. Uniform and Gaussian disc models are shown with green and yellow markers, respectively. Blue lines correspond to the radii encompassing 50, 80, 90, and 99% of the total RT model’s flux. The blue shaded areas represent the standard deviation of these radii across the rotational phase. The red shaded area indicates the inner (Rt) and outer radius (Rt + δr) of the RT model.

We find that the size derived from the uniform disc model is modulated around an average value of 3.5R_* corresponding to 90% of the Brγ emitting region. The size obtained by interferometry appears to be modulated by the position of the funnel flows close to the star, with a minimum located around phase 0.8. The Gaussian model exhibits the same modulation but with a lower amplitude (2.1 ± 0.4R_*) and appears sensitive to the magnetosphere’s innermost region, close to the 50% flux radius. The size derived from the uniform disc model emerges as being the most appropriate to recover the reference model size, accounting for at least 80% of the total flux emitted by the magnetosphere.

The derived orientations obtained from interferometry seem to be particularly representative of the position of the accretion funnel flow and the on-sky orientation of the Brγ emitting region (Fig. 5). The measured position angle agrees with the magnetosphere’s orientation, particularly when the shock faces the observer (phase = 0.5). Nevertheless, it appears somewhat hazardous to decipher the shape and orientation of the emitting region across the rotational cycle from this observable only, as different magnetospheric configurations can be described by a very similar interferometric model (e.g., phases 0.03 and 0.25). A stronger constraint on the orientation of the funnel flows arises from differential phase measurements.

5.3. Differential phases and photo-centre shifts

From the differential phases, we can derive the photo-centre shift between the continuum and the Brγ line emitting region. In the regime of marginally resolved sources, there is a direct relationship between the projected photo-centre displacement vector (P) and the phase along each baseline (Lachaume 2003):

$$\begin{array}{c}\hfill {\varphi }_i=-2\pi \frac{B_i}{\lambda }\mathbf{P},\end{array}$$(3)

where ϕ_i is the differential phase measured for the ith baseline, B_i is the length of the corresponding baseline, and λ is the effective wavelength of the spectral channel. A four-telescope beam-combiner such as GRAVITY gives access to six projected baselines that enable us to accurately retrieve the value and orientation of the photo-centre shifts in each spectral channel (Le Bouquin et al. 2009; Waisberg et al. 2017). Such a measurement results in a position-velocity plot of the displacement of the photo-centre across the Brγ line relatively to the continuum. This is illustrated in the bottom panels of Fig. 5.

The photo-centre shifts trace the accretion funnel flow’s direction and follow the stellar rotation. For instance, when the northern accretion shock (N-shock) is located behind the star (phase = 0), the accreting material falls onto the stellar surface in the direction of the observer. Accordingly, the photo-centre measured in the blueshifted part of the line profile (≃−75 km s⁻¹) lies on the blueshifted part of the velocity map, corresponding to the approaching funnel flow. Equivalently, the photo-centre measured in positive velocity channels of the line profile (≃+75 km s⁻¹) is shifted towards the receding funnel flow. In contrast, when the shock faces the observer (phase = 0.5), the velocity map goes from blue to red in the east-west direction and the photo-centre shifts recover this trend, as demonstrated at phase 0.47.

We can thus identify three main directions and shapes of the photo-centre shifts: linear north-south at phase ≃0 (N-shock behind), S shape at phases 0.25 and 0.69, and linear east-west at phase ≃0.5 (N-shock in front). The differential phase is, therefore, a key element in the recovery of the geometry and orientation of the line-emitting region, tracing the moving material along a rotational cycle.

5.4. Signal-to-noise considerations

As a proof of concept, the results presented above assume infinite signal-to-noise ratio. The goal is to predict the spectroscopic and interferometric signatures of the magnetospheric accretion process. Thus, the models predict typical visibility amplitudes ranging from 1 down to 0.97 (see Fig. 5). Such a modest interferometric signal requires a measurement accuracy of about 1% to be securely detected. Similarly, the models predict a deviation of the differential phases by 1–2 degrees (Fig. 5), which requires a level of accuracy of approximately a fraction of a degree to yield a robust detection. Recent interferometric studies performed with GRAVITY in the K band demonstrate that these levels of accuracy can be routinely obtained with reasonable exposure times on young stellar objects (e.g., Bouvier et al. 2020b; GRAVITY Collaboration 2020, 2023) or active galactic nuclei (GRAVITY Collaboration 2018).

6. Summary and conclusion

We present non-LTE radiative transfer modelling of the Brackett γ line emission for non-axisymmetric models of accreting magnetospheres. We used the equations of a misaligned dipolar magnetic field to derive the geometry of the magnetospheric accretion region for different obliquities of the magnetic dipole. We used MCFOST to compute radiative signatures of the Brγ line along a full stellar rotational cycle. Furthermore, we derived near-infrared interferometric observables for the line, comparable to what the GRAVITY instrument has already measured for T Tauri stars.

The main conclusions of this study are the following:

• 1)

  The total flux in the line and the line-to-continuum ratio depend on the obliquity of the dipole. As the obliquity increases, the size of the emitting region decreases, leading to a lower integrated flux. Also, projection effects make the emission region of lines forming close to the stellar surface appearing narrower.

• 2)

  The Brγ line total flux varies with the rotational phase due to the non-axisymmetry of the models induced by the magnetic obliquity. The line profiles exhibit a red-shifted absorption, that is, an inverse P Cygni profile, when a significant fraction of the accretion shock is aligned with the observer’s line of sight. When the shock is hidden on the opposite side of the star, the line profiles exhibit a double-peaked shape, reminiscent of the lines formed in a rotating envelope. The latter is due to the relatively large rotational velocity of the magnetospheric model (∼80 km s⁻¹).

• 3)

  Near-infrared interferometric observations in the Brγ line directly probe the size of the magnetospheric accretion region. The Gaussian disc model is sensitive to the brightest parts of the magnetosphere, up to 50% of the truncation radius, while a uniform disc model grasps 90% of the magnetosphere. It is of prime importance to consider this aspect when estimating magnetospheric radius from interferometric measurements. In both cases, the measured radius varies with the rotational phase (due to the non-axisymmetry of the dipole). A robust interferometric estimate of the magnetospheric radius therefore requires monitoring the system over a full rotational cycle.

• 4)

  The combined knowledge of the differential phase and the associated photo-centre shifts gives hints on the object orientation and geometry. More specifically, the relative direction of the photo-centre shifts indicates the changing orientation of the accreting material along the rotational cycle in the non-axisymmetric case.

Near-infrared interferometry of the Brackett γ line is used to characterise the inner star-disc interaction region and offers a good estima te of the size of the line’s forming region at sub-astronomical unit levels of precision. Comparing this size with reference model radii, such as the truncation radius, allows us to distinguish between multiple origins of the Brγ line, within or beyond these radii (e.g., magnetosphere, stellar and disc winds, jets). Furthermore, simultaneous spectroscopic and interferometric observations along a rotational cycle have the potential to unveil the geometry and orientation of the line’s emitting region. The variability of the line associated with the photo-centre shifts provides a unique and unambiguous proxy of the physical processes occurring in the magnetosphere of young accreting systems, within a few hundredths of an astronomical unit around the central star.

Acknowledgments

The authors thank Claudio Zanni, Lucas Labadie, Catherine Dougados, and Alexander Wojtczak for fruitful discussions. 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). B. Tessore thanks the french minister of Europe and foreign affairs and the minister of superior education, research and innovation (MEAE and MESRI) for research funding through FASIC partnership. C. Pinte acknowledges funding from the Australian Research Council via FT170100040 and DP180104235. The numerical simulations presented in this paper were performed with the Dahu supercomputer of the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities.

References

Appendix A: Benchmark

We present the comparison between line profiles obtained with MCFOST and previous studies. The magnetospheric model corresponds to the axisymmetric and compact configuration of Muzerolle et al. (2001) with a fixed shock temperature at 8 000 K, a rotation period of five days, and the following canonical T Tauri parameters: M_* = 0.8 M_⊙, R_* = 2 R_⊙ and T_* = 4 000 K. The inclination of the system is 60 degrees. The continuum emission of the stellar surface (shock and photosphere) is constant for all models. Figures A.1, A.2, A.3, and A.4 show the Hα, Hβ, Paβ, and Brγ lines profiles for different values of T_max and Ṁ. An inverse P Cygni profile with a redshifted absorption is seen for all lines, although it is dependent on the value of the mass accretion rate and of the maximum temperature. For a given mass accretion rate, an increase of the maximum temperature results in a higher line emission peak and a shallower redshifted absorption. As the temperature increases, the line source function increases, which is the cause of a higher emission above the continuum emission. The appearance of the red-shifted absorption component is caused by absorption from the gas above the stellar surface. It is controlled by the ratio between the source function of the line in the accretion funnel and that of the underlying continuum from the stellar surface, especially at low mass accretion rates and temperatures. Eventually, for the highest mass accretion rate and temperature, the lines become so optically thick that the redshifted absorption is washed out by the large wings of the line. The red-shifted absorption is more pronounced for lines forming closer to the accretion shock such as the Hβ line. At a temperature above 8,000 K and a mass accretion rate above 10⁻⁸ M_⊙ yr⁻¹, the continuum emission from the magnetosphere becomes significant and the line-to-continuum ratio decreases. This effect is seen, for instance, in the Hα line (see Fig. A.1). When the mass accretion rate increases for a given temperature, the density of the magnetosphere increases. As a consequence, the line source function increases. At a high temperature and high density, the background continuum emission of the magnetosphere dominates for certain wavelengths, and absorption occurs. The latter effect is seen in the Paβ (Fig. A.3) and Brγ (Fig. A.4) lines where the strong continuum contribution at the disc surface leads to absorption at low velocities; there, the line’s source function is small. These results are consistent with the previous studies and demonstrate the robustness of our code for modelling the close environment of T Tauri stars (Tessore et al. 2021).

Fig. A.1. Dependence of Hα line flux with mass accretion rates Ṁ and maximum temperatures Tmax. The inclination of the system is 60°.

Fig. A.2. Same as Fig. A.1, but for Hβ.

Fig. A.3. Same as Fig. A.1, but for Paβ.

Fig. A.4. Same as Fig. A.1, but for Brγ.

Appendix B: Derivation of the interferometric pure-line phase and visibility

We focused on the magnetospheric emission probed by the Brγ line and, therefore, aimed to remove any additional contributions (stellar photosphere, the accretion shocks, dusty disc, etc...). Following Kraus et al. (2008) and Bouvier et al. (2020b), we computed the continuum-subtracted observables, which are the so-called pure-line visibility and phase, using the emission line profiles computed in §3. This is of prime importance in the case of the Brγ line as the magnetospheric emission is quite faint in the infrared (≈1.3 excess flux compared to the continuum, Fig. 3). The derivation of the pure line quantities is only possible if the source is marginally resolved (i.e, size < λ/2B).

In this case, the pure line visibility V_line and phase Φ_line are given by

$$\begin{array}{c}\hfill V_{\mathit{Line}}=\frac{F_{L/C}{V}_{\mathit{Tot}}-V_{\mathit{Cont}}}{F_{L/C}-1},\end{array}$$(B.1)

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathit{Line}}=arcsin(\frac{F_{L/C}}{F_{L/C}-1}\frac{V_{\mathit{Tot}}}{V_{\mathit{Line}}}sin{\mathrm{\Phi }}_{\mathit{Tot}}),\end{array}$$(B.2)

where F_L/C denotes the line-to-continuum flux ratio as taken from the normalised spectrum (Fig. 3), V_Cont is the visibility computed in the continuum (star+shock only), and V_Tot, Φ_Tot are the total complex quantities measured by GRAVITY. In Eq. (B.1), we note that in the case when F_L/C is close to one, the derived V_Line cannot exist (converges to infinity). Such non-ideal profiles appear if the red absorption becomes too significant. Therefore, we decided to discard the affected spectral channels for phases 0.25 and 0.47, where F_L/C is too close to one: one point (v = 53 km s⁻¹) at phases 0.25 and two points (v = 15, 53 km s⁻¹) at phase 0.47.

Fig. B.1. Fourier sampling (uv coverage) of the simulated data. The different colours correspond to the six different baselines of the VLTI. The eight points per baseline represent a typical observational sequence with one data point per hour.

All Figures

	Fig. 1. Density distribution of non-axisymmetric model with an obliquity of 10°. The rotation axis of the star, Ω, is shown with a white arrow and the dipole axis, μ, is shown with a red arrow. The density is computed from Eqs. (1) and (2). The colour map scales with the density.
In the text

Fig. 2. Origin of emission seen across the Brackett γ line. The contribution of individual images to the total line flux is indicated on the central image showing the line profile. The brightness maps are in units of the maximum emission. The emission of the stellar surface is saturated. Orange to red colours indicate the regions of maximum emission. The system is seen at an inclination of 30° and an rotational phase of ∼0.25, similarly to Fig. 1.

In the text

	Fig. 3. Brackett γ line variability along the rotational cycle. Each column corresponds to a specific rotational phase. At phase 0, the shock area is unseen on the stellar surface, while at a phase of 0.5, the shock is fully seen on the visible hemisphere. The colours correspond to different values of the obliquity. All fluxes are computed with an inclination of 60°.
In the text

	Fig. 4. Radius encompassing 90% of total flux (R90) for each line as a function of the obliquity, βma. Hydrogen lines are labelled with different colours.
In the text

Fig. 5. Synthetic interferometric measurements and modelling across the rotational phase of the system. Top: integrated images over the Brγ line. Green (uniform disc) and yellow (Gaussian disc) ellipses are the characteristic sizes measured with a GRAVITY-like instrument. Middle-top: pure line visibility amplitude observables associated with the corresponding models. The visibility variation (for a given baseline) as the uv plane rotates is the specific signature of an elongated object. Middle-bottom: pure phase visibility across the line profile for the six baselines of the VLTI. Colours encode the observing time. Bottom: velocity map of the radiative transfer model. The coloured dots represent the measurement of the photo-centre derived from the phase visibility in each available spectral channel (see Appendix B). In each figure, the magnetosphere and the stellar surface have been normalised independently for display purposes.

In the text

Fig. 6. Interferometric radii as function of rotational phase. Uniform and Gaussian disc models are shown with green and yellow markers, respectively. Blue lines correspond to the radii encompassing 50, 80, 90, and 99% of the total RT model’s flux. The blue shaded areas represent the standard deviation of these radii across the rotational phase. The red shaded area indicates the inner (Rt) and outer radius (Rt + δr) of the RT model.

In the text

	Fig. A.1. Dependence of Hα line flux with mass accretion rates Ṁ and maximum temperatures Tmax. The inclination of the system is 60°.
In the text

	Fig. A.2. Same as Fig. A.1, but for Hβ.
In the text

	Fig. A.3. Same as Fig. A.1, but for Paβ.
In the text

	Fig. A.4. Same as Fig. A.1, but for Brγ.
In the text

	Fig. B.1. Fourier sampling (uv coverage) of the simulated data. The different colours correspond to the six different baselines of the VLTI. The eight points per baseline represent a typical observational sequence with one data point per hour.
In the text