11institutetext: I. Physikalisches Institut, Universität zu Köln, Zülpicher
Str. 77, 50937, Köln, Germany
11email<EMAIL_ADDRESS>22institutetext: Univ. Grenoble Alpes,
CNRS, IPAG, 38000 Grenoble, France 33institutetext: Max Planck Institute for
Astronomy, Königstuhl 17, 69117 Heidelberg, Germany 44institutetext: Max
Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, 85741
Garching bei München, Germany 55institutetext: Max-Planck-Institute for Radio
Astronomy, Auf dem Hügel 69, 53121 Bonn, Germany 66institutetext: Dublin
Institute for Advanced Studies, 31 Fitzwilliam Place, D02,XF86 Dublin, Ireland
77institutetext: INAF – Osservatorio Astronomico di Capodimonte, via
Moiariello 16, 80131 Napoli, Italy 88institutetext: Instituto de Astronomía,
Universidad Nacional Autónoma de México, Apdo. Postal 70264, Ciudad de México
04510, Mexico 99institutetext: CENTRA, Centro de Astrofísica e Gravitação,
Instituto Superior Técnico, Avenida Rovisco Pais 1, 1049 Lisboa, Portugal
1010institutetext: Universidade do Porto, Faculdade de Engenharia, Rua Dr.
Roberto Frias, 4200-465 Porto, Portugal 1111institutetext: School of Physics,
University College Dublin, Belfield, Dublin 4, Ireland 1212institutetext:
Universidade de Lisboa, Faculdade de Ciências, Campo Grande, 1749-016 Lisboa,
Portugal 1313institutetext: LESIA, Observatoire de Paris, PSL Research
University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Paris
Diderot, Sorbonne Paris Cité, France 1414institutetext: Sterrewacht Leiden,
Leiden University, Postbus 9513, 2300 RA Leiden, The Netherlands
1515institutetext: Dipartimento di Fisica & Chimica, Università di Palermo,
Piazza del Parlamento 1, I-90134 Palermo, Italy

# The GRAVITY Young Stellar Object survey

IX. Spatially resolved kinematics of hot hydrogen gas in the star/disk
interaction region of T Tauri stars

J. A. Wojtczak 11 L. Labadie 11 K. Perraut 22 B. Tessore 22 A. Soulain 22 V.
Ganci 115 5  J. Bouvier 22 C. Dougados 22 E. Alécian 22 H. Nowacki 22 G. Cozzo
221515 W. Brandner 33 A. Caratti o Garatti 336677 P. Garcia 991010 R. Garcia
Lopez 33661111 J. Sanchez-Bermudez 3388 A. Amorim 991212 M. Benisty 22 J.-P.
Berger 22 G. Bourdarot 44 P. Caselli 44 Y. Clénet 1313 P. T. de Zeeuw 441414
R. Davies 44 A. Drescher 44 G. Duvert 22 A. Eckart 1155 F. Eisenhauer 44 F.
Eupen 11 N. M. Förster-Schreiber 44 E. Gendron 1313 S. Gillessen 44 S. Grant
44 R. Grellmann 11 G. Heißel 1313 Th. Henning 33 S. Hippler 33 M. Horrobin 11
Z. Hubert 22 L. Jocou 22 P. Kervella 1313 S. Lacour 1313 V. Lapeyrère 1313
J.-B. Le Bouquin 22 P. Léna 1313 D. Lutz 44 F. Mang 44 T. Ott 44 T. Paumard
1313 G. Perrin 1313 S. Scheithauer 33 J. Shangguan 44 T. Shimizu 44 S.
Spezzano 44 O. Straub 44 C. Straubmeier 11 E. Sturm 44 E. van Dishoeck 441414
F. Vincent 1313 F. Widmann 44

(Received -; accepted -)

###### Abstract

Context. Hot atomic hydrogen emission lines in pre-main sequence stars serve
as tracers for physical processes in the innermost regions of circumstellar
accretion disks, where the interaction between star and disk is the dominant
influence on the formation of infalls and outflows. In the highly magnetically
active T Tauri stars, this interaction region is particularly shaped by the
stellar magnetic field and the associated magnetosphere, covering the inner
five stellar radii around the central star. Even for the closest T Tauri
stars, a region as compact as this is only observed on the sky plane at sub-
mas scales. To resolve it spatially, the capabilities of optical long baseline
interferometry are required.

Aims. We aim to spatially and spectrally resolve the Br$\gamma$ hydrogen
emission line with the methods of interferometry in order to examine the
kinematics of the hydrogen gas emission region in the inner accretion disk of
a sample of solar-like young stellar objects. The goal is to identify trends
and categories among the sources of our sample and to discuss whether or not
they can be tied to different origin mechanisms associated with Br$\gamma$
emission in T Tauri stars, chiefly and most prominently magnetospheric
accretion.

Methods. We observed a sample of seven T Tauri stars for the first time with
VLTI GRAVITY, recording spectra and spectrally dispersed interferometric
quantities across the Br$\gamma$ line at 2.16 $\mu$m in the NIR K-band. We use
the visibilities and differential phases to extract the size of the Br$\gamma$
emission region and the photocenter shifts on a channel by channel basis,
probing the variation of spatial extent at different radial velocities. To
assist in the interpretation, we also make use of radiative transfer models of
magnetospheric accretion to establish a baseline of expected interferometric
signatures if accretion is the primary driver of Br$\gamma$ emission.

Results. From among our sample, we find that five of the seven T Tauri stars
show an emission region with a half-flux radius in the four to seven stellar
radii range that is broadly expected for magnetospheric truncation. Two of the
five objects also show Br$\gamma$ emission primarily originating from within
the corotation radius, which is an important criterion for magnetospheric
accretion. Two objects exhibit extended emission on a scale beyond 10 R∗, one
of them even beyond the K band continuum half-flux radius of 11.3 R∗. The
observed photocenter shifts across the line can be either similar to what is
expected for disks in rotation, or show patterns of higher complexity.

Conclusions. Based on the observational findings and the comparison with the
radiative transfer models, we find strong evidence to suggest that for the two
weakest accretors in the sample, magnetospheric accretion is the primary
driver of Br$\gamma$ radiation. The results for the remaining sources imply
either partial or strong contributions coming from additional, spatially
extended emission components in the form of outflows, such as stellar or disk
winds. We expect that in actively accreting T Tauri stars these phenomena
typically occur simultaneously on different spatial scales. Through more
advanced modelling, interferometry will be a key factor in disentangling their
distinct contributions to the the total Br$\gamma$ flux arising from the
innermost disk regions.

###### Key Words.:

stars: formation - stars: variables: T Tauri, Herbig Ae/Be - accretion,
accretion disks - techniques: high angular resolution - techniques:
interferometric -

##  1 Introduction

Figure 1: A schematic depiction of the innermost parts of the circumstellar
disk of a T Tauri star.Br$\gamma$ emission in T Tauri stars is associated with
the concept of magnetospheric accretion, but outflows caused by certain types
of disk or stellar winds (dashed lines) can heat the hydrogen gas to the
necessary temperatures to produce Br$\gamma$ emission as well. Different
origin mechanisms can be tied to different spatial scales. Magnetospheric
accretion flows (depicted in red) can only stably form from within the
corotation radius R${}_{\textrm{co}}$, where the gaseous disk rotates at the
same angular velocity as the star. Br$\gamma$ radiation from within
R${}_{\textrm{co}}$ is likely to trace processes close to the truncation
radius R${}_{\textrm{tr}}$, whereas emission from regions on the scale of the
sublimation radius R${}_{\textrm{sub}}$ is tracing a disk/wind phenomenon.
This figure shows specifically the common case of a non-axisymmetric
magnetosphere, featuring a tilt between the stellar rotational axis
$\Omega_{*}$ and the magnetic dipole $\mu$. The obliquity of the magnetosphere
leads to the formation of two accretion columns, which will funnel gas from a
wide arc at the truncation onto one primary hotspot per hemisphere, as first
suggested in Bertout et al. (1988), Calvet & Hartmann (1992). The equations
and physical quantities relevant for the determination of the corotation
radius $R_{\rm co}$ and the magnetic truncation radius $R_{\rm tr}$ are
discussed in Section 5. For the equations used to estimate the sublimation
radius $R_{\rm sub}$, which indicates the area at which the disk becomes
sufficiently hot to destroy dust grains, we refer to Gravity Collaboration et
al. (2021b).

As the interplay between a young stellar object (YSO) and its circumstellar
disk has a major influence on the different phases of its development, the
study of such a star-disk interaction is of great importance in providing
insight into many aspects of star formation. Active research into this topic
is connected to a variety of astrophysical concepts, ranging from high level
models of stellar evolution down to the environmental conditions in the
cradles of planet formation and the mechanics that govern mass accretion
during the pre-main sequence stage of a YSO. A number of ideas regarding the
processes behind accretion and ejection of matter in the inner disk have been
extensively explored in theoretical models and are relatively well established
on a conceptual level. However, the exact nature of the physical interaction
in the innermost circumstellar region inside 0.1 au from the central star is
still debated because the observational evidence for many of these concepts
remains mostly indirect. In this context, a distinction arises between
intermediate-mass Herbig stars and low-mass T Tauri stars. The latter have a
lower effective temperature and weaker mass accretion rates, but also exhibit
significant magnetic field strengths of 1-2 kG due their convective nature
(e.g. Johns-Krull 2007). This substantial difference in the level of magnetic
activity has given rise to the paradigm of magnetospheric accretion. As a
theory concerning accretion in low mass YSOs, it has attained preeminence due
to its ability to reconcile accretion properties with the presence of stellar
hot spots and the observations of hydrogen line profiles (Camenzind, 1990;
Calvet & Hartmann, 1992). It has superseded the earlier theory of boundary
layer accretion, which is still thought to govern the processes of mass
accretion in magnetically weak intermediate mass Herbig Ae/Be type YSOs
(Mendigutía, 2020).
In the model of magnetospheric accretion, it is proposed that the inner
gaseous disk is truncated by the strong stellar magnetic field because of the
magnetic pressure exerted perpendicular to the magnetic field lines threading
the disk. Thus, a region of a few stellar radii surrounding the star is formed
where the disk plane is almost completely evacuated from gas particles
(Bouvier et al., 2007). Instead of further spiralling inwards, the gas
particles can only move freely along the field lines, which effectively serve
as funnels that channel the gas onto the star. As the gas impacts the stellar
surface at high latitudes, it creates heated shock regions close to the poles
(Hartmann et al., 2016).

Hydrogen emission lines serve as tracers for the kinematics of the hot gas in
the inner star-disk system. In particular the near-infrared Br$\gamma$ line at
2.16 $\mu$m, directly relevant for this work, requires the gas to be at
temperatures of 8000 to 10000 K (Muzerolle et al. (1998a), Tambovtseva et al.
(2016)) in order to produce strong emission. A heated gaseous disk in
Keplerian rotation can in principle be the cause for the Br$\gamma$ emission
in YSOs, in particular for intermediate-mass Herbig Ae/Be stars (Kraus et al.,
2008). However, even in those cases, theoretical considerations indicate that
the total line flux contribution coming from the hot disk layers is faint
compared to the contribution from outflows and infalls, see again Tambovtseva
et al. (2016). In T Tauri stars, with their much lower effective temperatures
and comparatively lower mass-accretion rates, the radiation field strength
alone is insufficient to heat up the gas disk to such temperatures at scales
of 0.1 au (Natta, 1993). It is then reasonable to assume that only specific
heating mechanisms acting much closer to the star in the accretion funnel
flows or in the stellar/disk wind launching regions will result in the
formation of Br$\gamma$ emission lines. Consequently, Br$\gamma$ emission in T
Tauri stars is an excellent tracer of processes associated to such gas heating
mechanisms, in particular for the magnetospheric accretion model which is
believed to play a central role in the evolution of low-mass pre-main sequence
stars.
Fig. 1 illustrates how different possible mechanisms producing Br$\gamma$
emission are tied to different spatial scales. Br$\gamma$ emission coming from
extended parts of the gaseous disk closer to the sublimation radius is likely
to trace outflows in the form of a disk or even a stellar wind, launched from
either the circumstellar disk or the polar region of the star, respectively
(Ferreira, 2013). For Herbig stars, a wind contribution in the Br$\gamma$
emission has been proposed in previous interferometric works (Garcia Lopez et
al., 2015; Kurosawa et al., 2016). On the contrary, Br$\gamma$ emission
arising from a more compact region close to the star is considered, in
particular in magnetically active T Tauri stars, to originate in the accretion
columns located inside the magnestopheric truncation radius $R_{tr}$, at which
point the ram pressure of the infalling material is compensated by the
magnetic pressure exerted perpendicular to the field lines (Bessolaz et al.,
2008). A useful parameter to differentiate compact versus extended emission in
the context of magnetospheric accretion can be the corotation radius $R_{\rm
co}$, which is the radius at which the angular velocities of the gaseous disk
and of the star equate (Bouvier, 2013). Since the formation of stable funnels
is only possible if the disk is truncated within the corotation radius
(Romanova & Owocki, 2016), Br$\gamma$ emission detected inside $R_{\rm co}$
indicates the probable presence of magnetospheric accretion flows. Following
Fig. 1, it is clear that the mechanisms involved in the accretion process of
matter in low-mass pre-main sequence stars are acting on a scale of a few (for
the magnetic accretion) to a few tens (in the case of winds) of stellar radii.
Prior to the advent of NIR long baseline interferometry, observational
evidence for the magnetospheric accretion scenario has mainly been discussed
in the context of high spectral resolution spectroscopy studies by Muzerolle
et al. (1998b), who have shown that the magnetospheric accretion model could
predict the Br$\gamma$ profiles. The process of magnetospheric accretion would
influence the shape and strength of hydrogen emission lines, with a
characteristic inverse P Cygni profile caused by infalling optically thick gas
( e.g. Calvet & Hartmann 1992, Edwards 1994,Petrov et al. 2014), typically
observed at high inclination of the system. Kurosawa et al. (2011)
investigated through theoretical simulations how the combination of the
magnetospheric model with a (stellar/disk) wind model can qualitatively
reproduce the observed profiles. However, without spatially resolving the
star-disk interacting regions, such a picture will remain incomplete since
such spectroscopic features with their overall line shapes will not be unique
to magnetospheric accretion.
With the advent of long baseline interferometry delivering sub-milliarcsecond
resolution and the concomitant improvement in sensitivity of the technique, a
larger reservoir of Br$\gamma$ emitting T Tauri stars at a few hundred parsecs
distance has become accessible. This has been exploited in the past in the
works of Eisner et al. (2010), Eisner et al. (2014), who examined Br$\gamma$
emission in a number of T Tauri stars with the Keck Interferometer, as well as
in recent pioneering studies with the VLTI by Gravity Collaboration et al.
(2017b) for S CrA, Gravity Collaboration et al. (2020) for TW Hya and Bouvier
et al. (2020) for DoAr 44.
In this paper, we report on observations made with GRAVITY (Gravity
Collaboration et al., 2017a), the second generation beam combiner instrument
of the VLTI, of a small sample of seven accreting T Tauri stars. We aim to
improve upon the results obtained in previous studies by combining the high
angular and spectral resolution offered by GRAVITY to present spectrally
dispersed measurements of the Br$\gamma$ emission region across the sample,
capable of identifying broader trends and categories between the different
objects. We contextualise our findings by drawing upon the latest results in
radiative transfer modelling of magnetospheric accretion in Br$\gamma$ ,
deriving sets of synthetic observables to compare our observations to.
Sections 2 and 3 describe the sample, observations and data reduction.
Sections 4 to 5 describe the modelling and the results. Comparisons with
synthetic data obtained from the models can be found in Section 6, along with
a wider discussion on magnetospheric accretion and the origin of Br$\gamma$ in
T Tauri stars

##  2 Scientific sample and observations

Table 1: Properties of the stars of the T Tauri sample.

Name | Sp. Type | $\text{T}_{\rm eff}$ | $d$ | $M_{*}$ | $L_{*}$ | $R_{*}$ | $\dot{M}_{acc}$ | Age | References
---|---|---|---|---|---|---|---|---|---
|  | [K] | [pc] | $[M_{\odot}]$ |  [$L_{\odot}$] |  [$R_{\odot}$] |  [10-8 $M_{\odot}$yr-1] | [Ma] |
AS 353 | K5 | 4450 |  399.6 $\pm$ 4.2 |  1 $\pm$ 0.2 |  3.9 $\pm$ 0.67 |  3.32 $\pm$ 0.29 | 7.9 - 45 |  $\leq$ 1 |  Rei et al. (2018),
|  |  |  |  |  |  |  |  | Prato et al. (2003)
RU Lup | K7 | 4050 |  157.5 $\pm$ 1.3 |  0.6 $\pm$ 0.2 |  1.64 $\pm$ 0.22 |  2.60 $\pm$ 0.17 | 1.8 - 30 | 1-2 |  Herczeg & Hillenbrand (2008),
|  |  |  |  |  |  |  |  | Siwak et al. (2016)
VV CrA SW | K7 | 4050 |  156.6 $\pm$1.2 |  0.6 $\pm$ 0.2 |  1.73 $\pm$ 0.24 |  2.67 $\pm$ 0.19 | 3.33 - 22 | 1-2 | Sullivan et al. (2019)
TW Hya | K6 | 4200 |  60.14 $\pm$ 0.07 |  0.9 $\pm$ 0.1 |  0.41 $\pm$ 0.05 |  1.21 $\pm$ 0.07 | 0.13 - 0.23 | 4-12 |  Donati et al. (2011),
|  |  |  |  |  |  |  |  | Venuti et al. (2019)
S CrA N | K5-K6 | 4300 |  160.5 $\pm$ 2.2 |  0.9 $\pm$ 0.2 |  3.82 $\pm$ 0.54 |  3.52 $\pm$ 0.25 | 7.8 - 50 |  $\leq$ 1 | Sullivan et al. (2019)
|  |  |  |  |  |  |  |  |  Gahm et al. (2018),
DG Tau | K6 | 4200 |  125.2 $\pm$ 2.3 |  0.8 $\pm$ 0.2 |  1.58 $\pm$ 0.37 |  2.37 $\pm$ 0.28 | 4.6 - 74 | 1-2 |  White & Ghez (2001),
|  |  |  |  |  |  |  |  | White & Hillenbrand (2004)
DoAr 44 | K2-K3 | 4840 |  146.3 $\pm$ 0.6 |  1.5 $\pm$ 0.2 |  1.61 $\pm$ 0.21 |  1.80 $\pm$ 0.12 | 0.63 - 0.9 | 3-6 |  Manara (2014),
|  |  |  |  |  |  |  |  | Espaillat et al. (2010)

111References are referring to the lower and upper limit $\dot{M}_{acc}$
literature values. $R_{*}$ is computed from effective Temperature and stellar
Luminosity. Other properties are adopted from Gravity Collaboration et al.
(2021b).

###  2.1 Sample

We have studied a total of seven T Tauri sources, which in the context of this
work will be referred to as the GRAVITY T Tauri sample. The targets in the
sample consist of a number of mostly strong accretors with mass accretion
rates on the order of $10^{-8}M_{\odot}\ yr^{-1}$, masses of less than 1.5
$M_{\odot}$, and stellar luminosities ranging up to 4 $L_{\odot}$. The
majority are mid K type stars with surface temperatures in the span between
4000 K and 4900 K and ages of 1 to 2 million years or younger. Their
individual properties are summarised in Table 1.

###  2.2 Observations

The individual sources of the sample were observed between 2016 and 2021 with
VLTI GRAVITY, the latest generation K-band beam combiner based at the ESO
facilities of the Very Large Telescope Inteferometer at Cerro Paranal, Chile.
The analysis of the Br$\gamma$ emission region relies exclusively on
observations taken with the four Unit Telescopes (UTs) with their 8.2 m
primary mirrors. The Auxilliary Telescopes (ATs), while providing better
coverage in the uv plane and potentially higher spatial resolution at the
longest baselines, are usually not sensitive enough to detect the relatively
faint emission signal of Br$\gamma$ originating from T Tauri type stars.
The standard UT configuration combines the 4 telescopes to a total of six
baselines, spanning up to 130 m in length on the ground and providing a
spatial resolution of up to $\frac{\lambda}{2B}$, which translates into 1.6 to
1.7 mas when observing in the near infrared between 2 and 2.4 micrometers.
With the exception of S CrA N, all of the sources were observed in single
field mode, meaning the fringe data were simultaneously recorded on the
separate detectors of the instrument’s fringe tracker (FT) and science channel
(SC).
The FT data were recorded at low spectral resolution (R=22) in 6 spectral
channels over the range of the entire K-band. While the fringe tracker data in
its own right can be utilised in the analysis of the K-band continuum, its
main purpose with respect to the investigation of the hydrogen emission lines
is to allow for longer integration times in the SC than would otherwise be
feasible by stabilising the fringes against small pertubations of the
atmosphere (Lacour et al., 2019). We made use of the FT capability to take the
SC observations in GRAVITY’s high resolution mode (R=4000) with a width of
$\sim$ 5 $\AA$, or 75 km/s, per channel, enabling us to spectrally resolve the
Br$\gamma$ feature and trace the gas kinematics over a range of up to $\pm$ 6
channels from the center of the line, depending on the strength of the feature
in the individual sources.
The observations of the different targets are subdivided into a set of short
term runs of about 5 minutes in length, each of which was saved as a separate
file, stretching over a total observation time of around 1 1/2 hours for most
of the objects, although the observation of DoAr 44 went on for notably longer
with a total time of slightly over 4 hours. Measurements of suitable
calibrator stars in a sufficiently close part of the sky were predominantly
taken either up to 45 minutes before the beginning or after the end of the
telescope run and thus potentially probing the sky under close but not
entirely identical atmospheric conditions. For the 2021 epoch of the RU Lup
observation, no calibrator close to the target was available for the night, so
that an object in a significantly different part of the sky had to be used for
data calibration.
Weather conditions for the relevant nights were mostly stable with typical
seeing values of around 0.5 to 0.7 arcseconds, however, due to the relatively
long obervation of DoAr 44, we note a larger fluctuation of the seeing
throughout the night and deteriorating atmospheric conditions towards the end
of the run.
The GRAVITY data on TW Hya, S CrA and DoAr 44 was previously published in
Gravity Collaboration et al. (2020), Gravity Collaboration et al. (2017b) and
Bouvier et al. (2020), respectively. A log of the observations of the T Tauri
sample is given in Table 2.

##  3 Data reduction and spectral calibration

Table 2: Log of observations.

Source | Date | Time (UTC) | Configuration | Calibrator | Airmass | Seeing [arcsec] | Files
---|---|---|---|---|---|---|---
DG Tau | 21.01.2019 | 01:31 - 03:06 | U1-U2-U3-U4 | HD 37491 | 1.59 | 0.65 | 10
TW Hya | 21.01.2019 | 06:59 - 08:24 | U1-U2-U3-U4 | HD 95470 | 1.03 | 0.51 | 13
AS 353 | 21.04.2019 | 08:44 - 10:22 | U1-U2-U3-U4 | HD 183442 | 1.31 | 0.34 | 12
RU Lup 2018 | 27.04.2018 | 03:23 - 04:30 | U1-U2-U3-U4 | HD 142448 | 1.27 | 0.57 | 8
RU Lup 2021 | 30.05.2021 | 01:58 - 03:19 | U1-U2-U3-U4 | HD 99264 | 1.31 | 0.64 | 8
S CrA N | 20.07.2016 | 05:30 - 06:13 | U1-U2-U3-U4 | HD 188787 | 1.10 | 0.56 | 6
VV CrA | 20.06.2019 | 02:53 - 04:04 | U1-U2-U3-U4 | HD 162926 | 1.31 | 0.77 | 3
DoAr 44 | 22.06.2019 | 02:04 - 05:46 | U1-U2-U3-U4 | HD 149562, | 1.04 | 0.68 | 28
|  |  |  | HD 147701, |  |  |
|  |  |  | HD 147578 |  |  |

222Time denotes the beginning of the first and the last file of the
observation.

###  3.1 Pipeline reduction and post-processing

The data were reduced and calibrated through the use of the standard GRAVITY
pipeline software provided by ESO (Lapeyrere et al. 2014). A full account of
the calibrated SC data from the observations is given in Appendix D.
In the context of this work we focus on the analysis of the SC data to extract
information about the Br$\gamma$ gas emission region, whereas FT data is only
used in a supplemental manner as described further below. For an overview of
the treatment of the K-band continuum, based on the same data sets obtained
with GRAVITY, we refer to the work of Gravity Collaboration et al. (2021b).
Each of the individual files obtained for all of our targets contains 4 sets
of flux data coming from the 4 Unit Telescopes, 6 sets of normalised
visibility amplitude data coming from the 6 UT baselines, 6 corresponding sets
of visiblity phase data, and 4 sets of closure phase data from the 4 unique
closed baseline triangles that can be formed by the UTs.
The individual files were merged before further analysis using the post-
processing routine of the GRAVITY pipeline, effectively time-averaging our
observables with a weighted mean over the 1 1/2 hours of observation. The data
were globally shifted to the local standard of rest (LSR) to eliminate
observational effects due to the earth’s rotation around the sun and the
radial velocity of the source.
The flux data were treated by first merging the 4 individual telescope spectra
and then normalising the Br$\gamma$ peak to the surrounding continuum via a
polynomial fit in the region between 2.1 $\mathrm{\mu}$m to 2.2
$\mathrm{\mu}$m . The normalised mean spectrum was further fitted with a
Gaussian function in a region (2.15 $\mathrm{\mu}$m to 2.18 $\mathrm{\mu}$m )
closer to the Br$\gamma$ peak to smooth out local fluctuations when
determining the line-to-continuum ratio on a channel to channel basis. For
those objects in the sample which exhibit asymmetric line shapes, a
combination of two single Gaussians was used to fit the spectrum. We then
consider the function value for a given wavelength as the true line to
continuum flux ratio for that wavelength and use this value in our
computations concerning the size and photocenter shift as detailed in Appendix
A. In order to ensure the interpretation of the data remains meaningful, we
implemented a selection criterion based on the strength of the flux ratio in
each channel relative to the continuum dispersion so that only those channels
that lie above a 2$\sigma$ threshold are taken into account for further
analysis.
For the visibility amplitudes, which we take as the root of the $V^{2}$, we
used FT data to normalise the science channel continuum around the Br$\gamma$
line to the continuum value measured by the fringe tracker for the
corresponding baseline. To this end, we took an average continuum value from
the FT channels, discarding the channel at the lowest wavelength as it is
affected more strongly by atmospheric conditions and the meteorology laser
operating at 1.908 $\mathrm{\mu}$m . To determine the visibility value in each
of considered the spectral channels, a single peaked Gaussian function was
then fitted to the data, similar to the treatment of the flux.
The visibility phase data received an analogous treatment after a polynomial
function was first fitted to and then subtracted from the data in order to
obtain the differential phase between line and continuum channels.
We attempted to correct for the telluric effects from atmospheric absorption
lines by using the calibrator spectra, but found that many of the calibrators
were not suitable for this method due to large photospheric absorption
features broad enough to cover the telluric features. An alternative attempt,
using an atmospheric model spectrum, to remove the telluric lines proved also
to be problematic. Instead, we decided against a full correction and opted to
restrain our analysis to a range of channels that falls inside the two
telluric features at 2.163446 $\mathrm{\mu}$m and 2.168645 $\mathrm{\mu}$m .
We proceeded to remove the channels affected by those two features from the
spectrum before attempting to fit the Gaussian profile as described above.
Whilst we recognize that the flux ratios obtained close to the telluric
features obtained in this manner are potentially incorrect, we propose that
the overall impact of this discrepancy on our results is negligible as low
signal strengths in those areas means we do not consider those channels in our
analysis for the majority of our targets.

###  3.2 High-precision spectral calibration

A more detailed analysis of the flux recorded by the 4 UTs has shown a
systematic spectral shift of the spectra with respect to known telluric
absorption lines in the K-band. This shift, which is on average on the order
of 5 Å, is not affecting the individual UT spectra equally, leading to a
misalignment between the 4 Br$\gamma$ lines recorded with the UTs. Such a
misalignment not only affects the flux values in the merged spectrum, but also
the computation of the complex visibilities by the pipeline during the
reduction process. This leads to a distortion of the Br$\gamma$ line signature
in both the visibility amplitudes and differential phases for the different
baselines, depending on the UTs involved. As a practical consequence of this
effect we noticed that the spectrally resolved line features in the visibility
amplitudes, which are generally expected to be single peaked, were not well
aligned with the peak of the mean spectrum and instead were for some baselines
shifted to the left or right, thus affecting the results of the channel by
channel analysis of those observables. The line feature in the differential
phases is equally affected, as S-shape signatures can be effectively
suppressed and appear single peaked.
We adjust for this effect by introducing a correction based on a fit of the 24
individual spectra of each observational run to the known telluric lines close
to Br$\gamma$ . If the telluric lines are too weak or the data too noisy to
produce a proper fit, a master correction based on an average shift is
applied. An example of this, showing the individual spectra before and after
the correction, is shown in Fig. 2.
We found the overall effect of the correction to be more pronounced in the
photocenter shift profiles derived from the differential phases, whereas the
visibilitiy data appeared to be relatively less affected. The general impact
is highest on the pure line quantities (see Appendix A), which is shown in
Fig. 3. While the marginally changed shape of spectra and visibility signals
after the correction will also affect the pure line observables, the dominant
effect here stems from the relative shift between the mean spectrum and the
interferometric observables. When computing $\mathrm{V_{Line}}$ from Eq. 5, a
different combination of $\mathrm{F_{L/C}}$ and $\mathrm{V_{tot}}$ values in
each channel is used after the correction, as now the visibility peaks are
better aligned with the mean spectral peak. The difference is immediately
obvious when comparing the pure line visibilities in Fig. 3 between top and
bottom. Before the correction, the peak of the visibility was shifted to the
right with respect to the four-telescope mean spectral peak, resulting in
comparatively large values for $\mathrm{V_{tot}}$ at much weaker flux ratios
in the blue wing. Consequently, the pure line visibilites are increasing as we
move to higher wavelengths. After the correction, the decreasing values for
$\mathrm{F_{L/C}}$ in the blue wing are paired with decreasing values in in
the visibilities, leading to a peak in pure line visibilites that is roughly
aligned with the peak in total visibilities.
The telluric shift has been detected in all of the data sets except for the
2021 epoch of RU Lup. The effect is currently understood to be connected to
the old grism of the science spectrometer and was consequently resolved with
the grism upgrade in October 2019, meaning all GRAVITY data taken after this
month should no longer be affected.

Figure 2: Top: The Br$\gamma$ spectra of the AS 353 observation, taken with
the four unit telescopes of the VLTI before the application of the spectral
calibration. The vertical line indicates the mean position of the line peak.
Before the correction, the individual spectra for most of our objects were
misaligned both with respect to the telluric lines and also between each
other. Bottom: The same spectra after they were brought into alignment through
the spectral calibration. The corrected spectra and recomputed observables
were then further shifted to bring them to the local standard of rest (LSR).

Figure 3: Top: Total and pure line visibilities for an exemplary baseline from
the AS 353 observation before the spectral correction was applied. Blue and
grey shaded regions indicate the uncertainties on the total and pure line
visibilities, respectively. Bottom: The same baseline after the correction.
Especially with respect to the computation of the pure line visibilities, the
change appears significant in the red wing of the signal.

##  4 Data and methodology overview

Through our GRAVITY observation we obtained K-band spectra between 2 and 2.4
microns, showing clear Br$\gamma$ features at around 2.16 microns with peak
line-to-continuum flux ratios spanning a range of values from 1.91 in the
extreme case of AS 353 down to only 1.19 for the faintest line we observed for
DoAr 44. The uncertainties on the flux ratios, derived from the instrumental
error bars given by the pipeline, are on the order of 0.01 for the majority of
the T Tauri sample, although in the case of VV CrA we note a severe increase
to around 0.03, even though the dispersion of the continuum flux values
remains comparable to the dispersion found for the other sources.
The continuum visibility data in the vicinity of the Br$\gamma$ feature for
the different objects in the sample ranges from visibilities of about 0.7 to
amplitudes close to 1, indicating that the total line emission region and the
surrounding continuum were partially resolved over the entire K-band at
baseline lengths between 40 m and 130 m for all of the targets. The Br$\gamma$
visibilities signals are single peaked and in most cases well aligned with the
position of the feature in the corresponding spectrum, see e.g. Fig. 22 in
Appendix D.
There are clear distinctions in terms of visibility signal strength at the
Br$\gamma$ line between the different sources. Particularly clear and strong
features were observed in both epochs of RU Lup, where the difference between
Br$\gamma$ peak and surrounding continuum can be as high as 0.12 (Fig. 23),
compared to the relatively faint signals detected at the position of
Br$\gamma$ for DoAr 44 with a line to continuum difference in visibility of
only about 0.01 at even the longest baselines.
Noise levels can equally vary greatly, with low continuum dispersion in
wavelength regions adjacent to the Br$\gamma$ line again for RU Lup, where we
measure a continuum standard deviation of 0.0023, compared to the very noisy
data obtained for VV CrA (Fig. 26) where the continuum standard deviation can
go up to 0.035.
We found non-zero differential phases at the position of Br$\gamma$ at several
baselines for each target in the sample, with typical absolute phase values of
around 0.5∘ to 1∘, although the largest measured differential phase goes up to
6.3∘. S-shaped features, as expected for systems in rotation, can be observed
for some sources at some baselines (e.g. RU Lup 2021, UT4-UT1, see Fig. 23),
but are rather exceptional and notably absent from many other targets. The
features appear single peaked in many cases, but the peak is not aligned with
the line center, suggesting that in such cases a potential S-shape might be
masked by one of the two constituent peaks being somehow attenuated.
The differential phase data is significantly more noisy when compared to the
visibilities. Standard deviations in the continuum adjacent to Br$\gamma$ ,
even for sources with clear differential phase Br$\gamma$ peaks, are on
average the order of 0.4∘. For the very noisy data of VV CrA the continuum
dispersion even goes up to 2∘ in standard deviation.
The closure phases at the position of Br$\gamma$ are consistent with the
closure phases of the nearby continuum throughout the sample with no clear
Br$\gamma$ signal visible for any of the sources. The continuum closure phases
are for most of our objects are close to 0, although small offsets of around
1∘ can be detected at certain baseline triangles for AS 353, RU Lup and DoAr
44 (Fig. 21, 23 and 27, respectively). For VV CrA (Fig. 26) we found again a
comparatively large dispersion of closure phase values, with standard
deviations across the different triangles reaching up to 4∘.
The data were treated by using the total line to continuum flux ratio and
continuum quantities to remove the influence of the continuum on the
interferometric observables. This way a ”pure line” Br$\gamma$ visibility and
differential phase was extracted from the total observables in each of the
chosen spectral channels and then further used to determine the size and
photocenter of the Br$\gamma$ emission region. Based on a simple geometric
Gaussian disk model (Berger & Segransan, 2007), we extracted the half-flux
radius (or half width half maximum, HWHM) in each channel, as well as the
relative offset of the photocenter with respect to the position of the
continuum photocenter. This treatment was applied to a number of spectral
channels based on the previously explained flux selection criterion to analyse
the change of those derived properties across the different velocity
components of the Br$\gamma$ line, see e.g. Fig. 6. A full explanation of the
process is given in Appendix A.

##  5 Results

###  5.1 General results

Figure 4: Top: The fitted HWHM of the Br$\gamma$ pure line emission region of
the central channel versus the mass accretion rate. The latter is depicted in
two forms: Once as a grey bar which denotes the range of $\dot{M}$ values
found in the literature on the respective source, then as the instantaneous
mass accretion rate, which was derived from the equivalent width of the
GRAVITY Br$\gamma$ spectra. The colour coding represents the classification of
the sample objects as either weak (green) or intermediate to strong accretors
(red). Most of the T Tauri stars cluster in a 4-7 R∗ range, which is on the
order of typical magnetospheric sizes.
Bottom: Br$\gamma$ HWHM versus accretion luminosities. The colour coding
corresponds to the upper plot. Table 3: Results for equivalent widths, mass
accretion rates and characteristic sizes of the Br$\gamma$ emission region.

| AS353 | RU Lup 2018 | RU Lup 2021 | TW Hya | DG Tau | DoAr 44 | S CrA | VV CrA
---|---|---|---|---|---|---|---|---
$\mathrm{F_{L/C}}$ | 1.91 | 1.86 | 1.83 | 1.61 | 1.29 | 1.19 | 1.53 | 1.65
$\mathrm{FWHM\ [km\ s^{-1}]}$ | 233 | 183 | 187 | 121 | 224 | 137 | 175 | 150
$\mathrm{W10\%\ [km\ s^{-1}]}$ | 594 | 607 | 420 | 241 | 445 | 274 | 524 | 416
$\mathrm{\mu_{HV}\ [km\ s^{-1}]}$ | -97 | -55 | -50 | -183 | - | - | -273 | -75.49
$\mathrm{\mu_{LV}\ [km\ s^{-1}]}$ | 3 | 2 | 6 | 1 | - | - | 3 | 0.84
$\mathrm{FWHM_{HV}\ [km\ s^{-1}]}$ | 520 | 586 | 295 | 482 | - | - | 549 | 598.54
$\mathrm{FWHM_{LV}\ [km\ s^{-1}]}$ | 204 | 154 | 154 | 116 | - | - | 170 | 141.32
EW [$\AA$] |  -19.62 $\pm$ 0.91 |  -16.82 $\pm$ 1.76 |  -13.47 $\pm$ 0.43 |  -7.11 $\pm$ 0.76 |  -5.80 $\pm$ 0.35 |  -2.65 $\pm$ 0.38 |  -8.71 $\pm$ 0.66 |  -10.09 $\pm$ 0.99
$\mathrm{L_{acc}\ [L_{\odot}]}$ |  3.07 $\pm$ 0.17 |  0.98 $\pm$ 0.12 |  0.76 $\pm$ 0.03 |  0.03 $\pm$ 0.00 |  0.33 $\pm$ 0.02 |  0.07 $\pm$ 0.01 |  1.74 $\pm$ 0.15 |  1.64 $\pm$ 0.19
$\mathrm{log\dot{M}_{acc,inst}\ [M_{\odot}yr^{-1}]}$ |  -6.38 $\pm$ 0.44 |  -6.76 $\pm$ 0.48 |  -6.87 $\pm$ 0.48 |  -8.74 $\pm$ 0.49 |  -7.40 $\pm$ 0.47 |  -8.46 $\pm$ 0.49 |  -6.56 $\pm$ 0.44 |  -6.52 $\pm$ 0.47
$\mathrm{\dot{M}_{acc,inst}\ [10^{-8}\ M_{\odot}yr^{-1}]}$ |  41.48${}_{-20.16}^{+41.35}$ |  17.31${}_{-8.54}^{+20.78}$ |  13.38${}_{-6.66}^{+15.79}$ |  0.18${}_{-0.1}^{+0.32}$ |  3.95${}_{-2}^{+4.44}$ |  0.35${}_{-0.19}^{+0.38}$ |  27.72${}^{+28.96}_{-13.57}$ |  29.77${}^{+34.39}_{-14.59}$
$\mathrm{PA\ [^{\circ}]}$ |  173 $\pm$ 3 |  99 $\pm$ 31 |  101 $\pm$ 31 |  130 $\pm$ 32 |  143 $\pm$ 12 |  137 $\pm$ 4 |  1 $\pm$ 6 |  91 $\pm$ 6
$\mathrm{i\ [^{\circ}]}$ |  41 $\pm$ 2 |  16${}^{+6}_{-8}$ |  20${}^{+6}_{-8}$ |  14${}^{+6}_{-14}$ |  49 $\pm$ 4 |  32 $\pm$ 4 |  27 $\pm$ 3 |  32 $\pm$ 3
$\mathrm{Continuum\ HWHM\ [au]}$ |  0.28 $\pm$ 0.05 |  0.21 $\pm$ 0.06 |  0.21 $\pm$ 0.06 |  0.042 $\pm$ 0.003 |  0.13 $\pm$ 0.01 |  0.16 $\pm$ 0.02 |  0.17 $\pm$ 0.02 |  0.16 $\pm$ 0.01
$\mathrm{Br\textrm{$\gamma$}\ HWHM\ [au]}$ |  0.206 $\pm$ 0.003 |  0.076 $\pm$ 0.001 |  0.061 $\pm$ 0.001 |  0.029 $\pm$ 0.001 |  0.143 $\pm$ 0.004 |  0.050 $\pm$ 0.002 |  0.071 $\pm$ 0.011 |  0.081 $\pm$ 0.005
$\mathrm{R_{tr}\ [au]}$ | 0.04 - 0.08 | 0.04 - 0.07 | 0.04 - 0.08 | 0.05 - 0.07 | 0.05 - 0.09 | 0.09 - 0.11 | 0.05 - 0.10 | 0.04 - 0.07
$\mathrm{R_{co}\ [au]}$ | 0.10 | 0.04 | 0.04 | 0.05 | 0.06 | 0.05 | 0.04 | 0.10
$\mathrm{Continuum\ HWHM\ [R_{*}]}$ |  18.13 $\pm$ 3.59 |  17.36 $\pm$ 5.10 |  17.36 $\pm$ 5.10 |  7.47 $\pm$ 0.70 |  11.32 $\pm$ 1.61 |  19.07 $\pm$ 2.69 |  10.38 $\pm$ 1.42 |  12.64 $\pm$ 1.04
$\mathrm{Br\textrm{$\gamma$}\ HWHM\ [R_{*}]}$ |  13.32 $\pm$ 1.16 |  6.32 $\pm$ 0.43 |  5.01 $\pm$ 0.35 |  5.13 $\pm$ 0.33 |  12.93 $\pm$ 1.56 |  5.93 $\pm$ 0.47 |  4.36 $\pm$ 0.75 |  6.56 $\pm$ 0.59
$\mathrm{R_{tr}\ [R_{*}]}$ | 2.83 - 5.30 | 3.28 - 6.14 | 3.53 - 6.61 | 8.31 - 12.35 | 4.50 - 8.42 | 10.62 - 12.68 | 3.36 - 6.29 | 2.86 - 5.36
$\mathrm{R_{co}\ [R_{*}]}$ | 6.47 | 3.31 | 3.31 | 8.89 | 5.44 | 5.96 | 2.44 | 7.25

333FL/C is the total line to continuum flux ratio at the center of the line.
FWHM and W10$\%$ are the line width at 50% and 10% of the peak flux,
respectively. $\mu$ is the position of the high (HV) and low velocity (LV)
component from the double Gaussian fit to the spectrum.
The Br$\gamma$ half-flux radius (HWHM) is given for the central channel. Rco,
continuum HWHM, PA, and i were adopted from Gravity Collaboration et al.
(2021b). Accretion luminosities were estimated from the line luminosities
based on the empirical relationship in Eq. 1. The instantaneous mass accretion
rates were then computed as in Eq 2. Rtr is given based on Eq. 4 for magnetic
field strengths between 1 kG to 3 kG and using the mass accretion rates
derived from the equivalent widths (EW). The asymmetric error bars on the
absolute values of $\mathrm{\dot{M}_{acc,inst}}$ represent the 25th and 75th
percentiles of the distribution. Other error bars denote 1$\mathrm{\sigma}$
uncertainties.

In this section we will first present a breakdown of the general results
obtained from the GRAVITY T Tauri data, introduce some additional concepts
relevant to the interpretation and give an overview over the global trends we
were able to identify across the entire sample. In the following sections we
will then proceed to go over each source individually in detail and
contextualise our findings by drawing upon relevant references, if and where
available. Table 3 summarises key results concerning the characteristics of
the Br$\gamma$ emission region and accretion properties, respectively. Figures
4 and 5 serve as a visual overview of the sample results with regards to the
most important criteria we use to interpret our results.
The sample consists of mostly strong Br$\gamma$ emitters, with normalised peak
line to continuum flux ratios between 1.91 and 1.53 for 6 out of our 8 data
sets and ranging down to 1.19 for the weakest emitter DoAr 44. Many of the
objects exhibit line shapes that are asymmetric, featuring an excess of
blueshifted emission. We approximate this asymmetry by fitting a superposition
of two Gaussian functions to the data, basically dividing the feature into a
narrow centered low velocity component and a broader offset high velocity
component. Whether these Gaussian components can necessarily be directly
attributed to specific physical phenomena, for example such as additional
outflows on top of an accreting magnetosphere, is not certain. Inverse P Cygni
features, as expected to be seen when observing magnetospheric accretion
columns at low inclination, are generally absent from the line profiles of the
sample, even for those objects observed at a close to pole-on configuration,
such as TW Hya. The width and offset of the broad component, and thus of the
total emission line, can vary significantly throughout the sample, even
between two epochs of the same target (RU Lup), although the overall strength
of the asymmetry correlates with the peak flux ratio and is also reflected in
the equivalent widths. We computed the line luminosities for each object from
the equivalent widths and K band magnitudes from 2MASS (Cutri et al., 2003).
We then used the empirical relationship from Alcalá et al. (2014)

$\displaystyle log\left(\frac{L_{acc}}{L_{\odot}}\right)$
$\displaystyle=a\cdot log\left(\frac{L_{Line}}{L_{\odot}}\right)+b$ (1)

with a=1.16 $\pm$ 0.07 and b=3.6$\pm$ 0.38 to estimate the accretion
luminosities and subsequently determined the corresponding mass accretion
rates via (Hartmann et al., 1998)

$\displaystyle\dot{M}_{acc}$
$\displaystyle=\left(1-\frac{R_{*}}{R_{in}}\right)^{-1}L_{acc}\frac{R_{*}}{GM_{*}}\approx
1.25L_{acc}\frac{R_{*}}{GM_{*}}$ (2)

which describes the relationship between the gravitational energy released and
radiated away by the mass falling in from the truncation radius. For the right
hand side approximation we follow the common convention of using Rin $\approx$
5 R∗ as a typical size for truncation radii.
These results are summarised in Table 3. They are also shown in Fig. 4, where
we present both the instantaneous mass accretion rates as well as a range
based on the lowest and highest previous estimates found in the literature.
Based on this comparison, we found that our mass accretion rates tend to fall
into the range established by past observations, even if mostly towards the
higher end for the stronger accretors. The uncertainties on the coefficients
a,b of Eq. 1 translate into relatively large and asymmetric error bars on the
absolute value of $\mathrm{L_{acc}}$ and consequently also
$\mathrm{\dot{M}_{acc}}$. We used a Monte-Carlo approach to determine the
final uncertainties on the accretion rate rather than Gaussian error
propagation not applicable here.
We fitted the visibility data with a 2D Gaussian in order to find the half-
flux radii (HWHM) of the pure Br$\gamma$ emission region, using the
inclination and position angles of the K-band continuum emission region
presented in Gravity Collaboration et al. (2021b) to constrain those two
parameters in our fit. The HWHMs we obtained through this process are, across
the entire sample, on the order of less than 0.21 au at their greatest extent,
ranging down to regions as compact as 0.03 au. In terms of stellar radii, the
central channel sizes vary from up to 13 R∗ to less than 5 $R_{*}$, although
the majority of objects show clustering (see Figures 5, 4) in a range from 4
to 7 R∗, which is fairly typical for magnetospheric radii (Bouvier et al.,
2007). This clustering can be observed across the entire range of mass
accretion rates of our small sample, so that a correlation between those sizes
and accretion rates cannot be established here. Equally, within our sample we
do not detect a clear correlation of those Br$\gamma$ sizes with the
luminosity of the objects, as opposed to the NIR continuum sizes, which were
shown to correlate as per R $\propto$ L1/2 (Gravity Collaboration et al.,
2021b). However, it is apparent from Fig. 5 that the two weakest accretors in
the sample (TW Hya, DoAr 44) also show most clearly an emission region smaller
than the corotation radius, whereas most of the other objects show Br$\gamma$
emission coming from beyond the corotation radius.
Changes in HWHM in the spectral channel across the Br$\gamma$ line can be
significant, with largest to smallest size-in-channel ratios being around 1.2
to 1.5 for many of the objects, although in extreme cases this can go up to
almost 2, whilst others show a relatively flat profile relative to their error
bars. Under a magnetospheric accretion scenario, a decrease of size towards
the edges of the line (i.e. higher velocities) would be naively expected, as
the gas at the highest velocities should be in free fall close to the star. In
our sample, only one source (TW Hya) clearly shows this decrease, whereas most
of the other objects show either a straight increase with a minimum at the
center channel or even a mixed profile with a decreasing size at low
velocities and then an increase as we move to higher velocities. These mixed
profiles could be indicative of a multi-component Br$\gamma$ emission region,
although in many cases the extent of these more complex signatures is within
the uncertainties.
An important point of reference in the search for the origin of Br$\gamma$
emission is the corotation radius, since we know that the inner disk must be
truncated inside the Rco so that stable magnetospheric accretion columns can
form. If the dominant driver of Br$\gamma$ emission is magnetospheric
accretion, we expect that the Br$\gamma$ emission region would be on the order
of the corotation radius or even slightly more compact. The corotation radius
is defined as the radius at which the angular velocity of the rotating disk
matches the angular velocity of the star $\Omega_{*}$ (Bouvier, 2013):

$\displaystyle R_{co}$
$\displaystyle=\left(\frac{GM_{*}}{\Omega^{2}_{*}}\right)^{1/3}$ (3)

Stellar angular velocities are determined from measurements of the rotational
period of the star or the v sin(i). We refer to Gravity Collaboration et al.
(2021b) for the relevant values for our sample objects. The corotation radii
are given in Table 3. Of course, the best way to determine whether the size of
the Br$\gamma$ emission region corresponds to the size of the magnetosphere in
any object would be to compare the region size against the magnetospheric
truncation radius directly. The truncation radius can be computed from stellar
parameters, if the strength of the magnetic dipole field is known:

$\displaystyle R_{tr}$
$\displaystyle=12.6\frac{B_{3}^{4/7}R_{2}^{12/7}}{M_{0.5}^{1/7}\dot{M}_{-8}^{2/7}}$
(4)

In this notation we follow the convention used in e.g. Hartmann et al. (2016).
$M_{0.5}$ is the stellar mass in units of 0.5 M⊙, R2 the stellar radius in
units of 2 R⊙, B3 the surface field strength of the dipolar magnetic field at
the stellar equator in kG and $\dot{M}_{-8}$ the mass accretion rate in units
of $\mathrm{10^{-8}M_{\odot}yr^{-1}}$ . The resulting truncation radius is
given in units of R⊙. The derivation is conceptually discussed in Bessolaz et
al. (2008).
Measurements of magnetic field strengths in T Tauri stars are scarce, so that
we cannot rely on the truncation radius as a primary criterion to identify
strong magnetospheric accretion cases. Even so, we can still use general
estimates of Rtr as complementary information. Building upon the approach laid
out in Gravity Collaboration et al. (2021b), we set a range of 1 - 3 kG for B
and compute the truncation radii for this range as presented in Table 3 from
the relevant stellar properties.
An overview of the emission region sizes in the central line channels,
compared to the respective corotation radii and the NIR K-band continuum sizes
taken from Gravity Collaboration et al. (2021b), can be found in Fig. 5. The
sample can generally be divided into targets with emission regions on the
order of the corotation radius (DoAr 44, TW Hya), those with more extended
Br$\gamma$ emission than the corotation radius, but still close to typical
magnetospheric radii around 5 R∗ (RU Lup, S CrA N, VV CrA SW), and those with
highly spatially extended Br$\gamma$ regions outside the typical range of
magnetospheric sizes (AS 353, DG Tau). This separation partially corresponds
to the the distinction between weak accretors ( $\dot{M}_{acc}$ ¡ 1$\cdot$
10${}^{-8}M_{\odot}yr^{-1}$, again DoAr 44 and TW Hya) and the intermediate to
strong accretors, which we indicated in Fig.4 via the different colours. While
the weak accretors do appear to conform to the corotation criterion with
compact emission regions, the situation is more complex for objects with
higher accretion rates, as there is no clear correlation between the relative
region sizes and the mass accretion rates for those targets.
Magnitudes for the observed photocenter displacements are on the order of 0.07
au or below. For multiple targets the entire photocenter profile across the
line is sufficiently compact to be within or close to 1 R∗, although such a
compact concentration of photocenters does not necessarily correlate with the
size of the emission region. The photocenter profiles are complex and not
easily classified, although a broad distinction can be made between targets
that show quasi-rotational profiles, featuring distinct blue and red arms
aligned along a common axis, targets with compact profiles that lack any clear
structures, and complex profiles which again appear potentially as a
superposition of multiple components. In Fig. 5 we give an overview over the
possible alignment of the photocenter shifts across the line with different
axes. To account for the offset of the profiles from the plot center, we
compute the difference vectors between each point and the center channel
(shown as a white square) and then compute the minimum angular difference
between the angle of the difference vector and the NIR continuum disk axis
position angle. The resulting value is a measure of how well aligned the
photocenter offsets are relative to the disk axis. The plot can be interpreted
in a meaningful way by considering the degree of clustering, and at which
angles the clustering occurs. A concentration of points in a narrow angular
range (AS 353) indicates that the photocenter shift is aligned along a single
axis. If the clustering range is close to zero, the photocenters are aligned
with the disk semi-major axis, which is something that would be expected from
a disk in Keplerian rotation. If the points cluster at 90∘, the photocenters
might be aligned with some form of polar outflow (red wing in DG Tau). An even
spread of points indicates an absence of recognisable structures and a random
distribution of photocenters (TW Hya). The figure clearly shows the previously
introduced distinction between simple alignments and complex patterns, where
we see multiple clusters (DG Tau, RU Lup).
It is important to remember that the differential phases track the relative
offset of the Br$\gamma$ emission region from the overall photocenter of the
continuum. In e.g. Fig. 6 the center of the plot is the combined photocenter
of the star and other continuum emitting components such as the dusty disk. If
the dusty disk is perfectly centrosymmteric, its own photocenter will coincide
with that of the star and the overall continuum photocenter will match the
position of the star. If the continuum photocenter however is not
centrosymmetric, either offset from the star or with an asymmetric brightness
distribution, then the zero position in those plots can significantly deviate
from the stellar position. We therefore make use of the closure phases to
extract information about the continuum asymmetry. For the majority of our
sample objects the closure phases themselves are close to zero and thus
suggest a centrosymmetric continuum. We then instead use the uncertainty given
by the phase dispersion to derive an uncertainty on the stellar position. To
do this we follow Gravity Collaboration et al. (2021b) by fitting an
azimuthally modulated ring model to the continuum closure phase data. The
results are depicted as red circles on the photocenter profile plots, showing
that indeed in the majority of the cases the center channel photocenter could
conceivably be centered on the stellar position.

Figure 5: Top: Br$\gamma$ half-flux radii (yellow), compared to K-band NIR
continuum half-flux radii (red) and corotation radii (green) across the
sample. Note that Br$\gamma$ region and continuum HWHMs are based on a
Gaussian disk and a Gaussian ring model, respectively, which should
approximately reflect the differences in the morphology of both regions. The
objects are given in order of descending mass accretion rate.
Bottom: The distribution of photocenter shift angles across the line for the
entire sample. The x-axis shows the minimum difference angle between the the
relative shift vector and the position angle of the continuum disk. The white
channel serves as the point of reference for the relative shift vector. An
angle of 0∘ means the photocenter shift vector is aligned with the disk axis,
as would for example be expected for a disk in rotation. A difference angle of
90∘ indicates that the shift vector is perpendicular to the disk axis. A
clustering of points close to a certain difference angle indicates an
alignment of the photocenter shifts along a preferential axis. VV CrA shows a
monoaxial distribution of photocenters, while for other targets more complex
profiles with multiple alignments are indicated. The number stated in brackets
after the object name gives the average error on the angle for the five
central channels.

###  5.2 Individual results

####  5.2.1 AS 353

Figure 6: Spectrum (top), size (middle) and photocenter shift (bottom)
profiles for AS 353. Shaded regions indicate uncertainties.

With a peak line to continuum flux ratio of 1.91 and an equivalent width of
-19.6 Å, AS 353 (Fig. 6) is the strongest Br$\gamma$ emitter in the GRAVITY T
Tauri sample. These measurements are consistent with the results of a
spectroscopic survey of young binaries done in 1996 by Prato et al. (2003),
who found an equivalent width of -21.1 Å for the Br$\gamma$ emission line of
this object. We do not find a significant change in size across the Br$\gamma$
line for the emission region, with all channels appearing consistent within
their respective uncertainties with a half flux radius of 0.206 au at the
center of the line. This puts the emission region at about twice the size of
the corotation radius at 0.1 au, but still well within the NIR continuum HWHM
of 0.28 au. However, when comparing these quantities, it is important to point
out a certain amount of confusion surrounding the distance estimate to the AS
353 system. Since the system is thought to be connected with the Aquila star
forming region (Rice et al., 2006) and thus to the complex of dark
interstellar clouds of the Aquila rift, accurate measurements of the distance
have historically been difficult to obtain. Previous works such as Edwards &
Snell (1982) refer to a distance of 150 pc, which was also adopted by Prato et
al. (2003) with a estimated error of 50 pc based on the idea that AS 353 is in
front of the Aquila rift at 200 $\pm$ 100 pc (Dame & Thaddeus, 1985). More
recent high accuracy parallax measurements of AS 353 taken by GAIA and
released as part of DR2 and DR3 put the binary at a comparatively remote 400
pc (Gaia Collaboration, 2020), which is the value we adopt in this paper. This
is supported by other recent studies covering parts of the Aquila rift, such
as the work of Ortiz-León et al. (2017), which presents astrometric
measurements that put the WS40/Serpens cloud complex at 436 $\pm$ 9.2 pc.
Still, even relatively recent publications like the interferometric study of
Eisner et al. (2014) use the 150 pc estimate to convert their angular scales
into physical sizes. For the Br$\gamma$ emission region size Eisner et al.
(2014) report an angular ring diameter of 1.06 mas and a subsequent ring
radius of 0.08 au, which, if converted to the more recent distance
measurements, would translate to a radius of 0.21 au. This, within a 1$\sigma$
error bar, agrees with our findings, despite the fact that their model does
not take into account the relatively high inclination of 41∘ of the inner disk
of AS 353. Conversely, if our results were to be converted into a physical
size based on a 150 pc distance, we would obtain a half flux radius of 0.077
au, i.e. something more compact than the corotation radius. However, based on
the fact that the most recent distance measurements for AS 353 independently
came to similar results, we consider the 400 pc estimate to be robust. As
such, the very large extent of the emission region indicates that the
Br$\gamma$ emission is likely produced through some form of disk wind outflow
launched from the inner gaseous disk.
This could be supported by the spatial profile of the Br$\gamma$ line
photocenters, which shows similarity to a rotational profile, albeit not
aligned with the disk and severely truncated in the blue wing. It is possible
that such an asymmetric structure, with a maximum photocenter displacement of
0.039 au, or 2.5 $R_{*}$, in the red wing, is caused by the interplay of disk
and jet. The AS 353 system is connected to the Herbig-Haro object HH32 and
drives a large scale outflow at a very high inclination of 70∘ (Curiel et al.,
1997) with a fainter blue component and multiple more prominent redshifted
knots. At smaller scale, H$\alpha$ emission originating from a blueshifted jet
was reported at a position angle of (111 $\pm$ 18)∘ by (Takami et al., 2003).
This is, within the respective uncertainties, close to perpendicular to the
NIR continuum disk axis of (173 $\pm$ 3)∘ (Gravity Collaboration et al.,
2021b).
Assuming that the weaker emission in the blue knot observed in HH32 translates
to a similar asymmetry for the small scale jet, a superposition of two
distinct photocenter shift profiles, one aligned with the jet axis and the
other coming from the rotating base of a disk wind and thus aligned with the
disk semi-major axis, could lead to the observed distribution of photocenters.
This is also supported by the closure phase data, which indicates that the
position of the star is consistent with the centroid of the central velocity
channel within an uncertainty of 0.026 au.

####  5.2.2 RU Lup

Figure 7: Spectrum (top), size (middle) and photocenter shift (bottom)
profiles for the two epochs of RU Lup in 2019 and 2021. Shaded regions
indicate uncertainties.

RU Lup was observed in two distinct epochs in 2018 and 2021 (Fig. 7). It
appears as the second strongest Br$\gamma$ emitter of our sample of T Tauri
YSOs, exhibiting a peak flux line to continuum ratio of 1.86 and 1.83,
respectively. From the Br$\gamma$ line we compute an equivalent width of 16.8
$\AA$ in 2018 and 13.5 $\AA$ in 2021, which reflects the narrower shape of the
feature in the more recent epoch. Using the empirical relationship presented
before, we use the Br$\gamma$ EW to estimate the accretion luminosity to be
(0.98 $\pm$ 0.12) L⊙ and (0.76 $\pm$ 0.03) L⊙, respectively. By comparison,
the X-Shooter measurements from 2012 presented in Alcalá et al. (2014) put the
equivalent width at 9.6 $\AA$ and the accretion luminosity of the RU Lup
system at 0.5 L⊙. Using the accretion luminosity, stellar radius and stellar
mass of RU Lup, we compute the instantaneous mass accretion rate to be
17.31${}^{+20.78}_{-8.54}$ and 13.38${}^{+15.79}_{-6.66}$ $\cdot\ 10^{-8}\
\mathrm{M_{\odot}yr^{-1}}$, while they arrive at a more moderate accretion
rate of 4.3 $\cdot\ 10^{-8}\ \mathrm{M_{\odot}yr^{-1}}$ largely due to their
significantly higher estimate of the stellar mass. Similar accretion rates can
be found in Manara (2014), although higher rates similar to the ones we
derived were also presented in the past (Siwak et al., 2016).
The Br$\gamma$ emission region was marginally resolved across all channels at
even the shortest baselines, giving us an extremely clear size profile across
different velocity components of the line. For 2018 we detect a minimum half
flux radius of 0.076 au (6.3 R∗) at the line center, which increases to 0.1 au
(8.3 R∗) at the edges of our selected channel range. A similar size profile
presents itself in 2021, where the central minimum is at 0.061 au (5 R∗),
while in the extreme blue channel the half flux radius increases again to 0.1
au, but only 0.088 au (7.3 R∗) au in the red wing. Such changes on the order
of 30% to 50% are particularly noteworthy when taking into account the high
signal to noise in both visibilities and flux data, and the resulting very
small uncertainties we obtain for the half flux radii. While the region in
both epochs is significantly more compact than the NIR continuum at 0.21 au
(17.4 R∗), Br$\gamma$ emission does appear to extend beyond the corotation
radius of 0.04 au (3.3 R∗) at both epochs. Both the absolute size, the
relative size variation and the pronounced blueshifted excess emission are not
characteristic for a scenario where Br$\gamma$ emission is exclusively, or
even dominantly, driven by magnetospheric accretion itself, but rather imply
the existence of a more extended outflow component, possibly in the form of a
disk or stellar wind.
The Br$\gamma$ signals in the differential phases are relatively weak and
translate into a distribution of different velocity photocenters that is quite
complex and not easily interpreted. While the maximum shift is on the order of
less than 1 R∗, there seem to be substructures in the profile that align
themselves along 2-3 different axes. One for the low velocity channels which
are, within a 1$\sigma$ uncertainty, broadly aligned with the semi-major axis
of the inner NIR continuum disk at (99 $\pm$ 31)∘ and also the large scale
outer disk as measured with ALMA at 121∘ (Huang et al., 2018). Then there are
two parallel axes which run almost perpendicular to the low velocity axis,
although between the two epochs, the red and blue arm are not equally well
aligned. The overall profile gives the appearance of a ”crescent” type shape
which seems to rotate by about 45∘ towards North from 2018 to 2021. We find
that in 2018 the blue arm of the profile is almost perfectly aligned with the
average PA of the blueshifted jet at (229 $\pm$ 10)∘ (Whelan et al., 2021),
while in 2021 this is true for the red arm. It is noticeable that the arm
closest to the jet axis in each case is relatively straight while the other
arm appears somewhat distorted in shape, although the effect is within the
uncertainty on the photocenter displacements. While it is tempting to conclude
that some form of polar outflow is responsible for this behaviour, a stellar
wind would not explain why the blue arm runs anti-parallel to the jet
direction in 2018. The complexity of this case would require a better
understanding of what a superposition of different photocenter distributions
caused by different Br$\gamma$ origin mechanisms would look like in order to
comment on this behaviour in a more meaningful way.

####  5.2.3 VV CrA SW

Figure 8: Spectrum (top), size (middle) and photocenter shift (bottom)
profiles for VV CrA SW. Shaded regions indicate uncertainties.

The GRAVITY K-band spectrum for the VV CrA (Fig.8) primary shows a Br$\gamma$
emission feature that lacks the pronounced blueshifted excess detected in many
of the other objects in the sample. The line shape appears more Lorentzian,
featuring broader wings which contain more of the flux than in a Gaussian
profile. Flux ratios vary from 1.65 at the center of the line to 1.05 at the
edges of our selected channel range. We computed a Br$\gamma$ equivalent width
of 10.1 $\AA$, which appears significantly weaker than than the 30 $\AA$
reported by Sullivan et al. (2019) based on archival K-band NIRSPEC data from
2008. Even by visual comparison, the GRAVITY and NIRSPEC spectra around the
Br$\gamma$ line appear quite distinct, as the NIRSPEC data shows both a peak
flux ratio of about 2 and a higher excess of blueshifted emission, leading to
an asymmetric shape similar to those observed for other objects of the GRAVITY
sample. By contrast, Scicluna et al. (2016) present Br$\gamma$ spectra taken
with CRIRES in 2013 which are visually closer to ours, also displaying quasi-
Lorentzian line shape that is slightly skewed towards the blue wing. Their
equivalent width of 9 $\AA$ is also consistent with our findings, although due
to their very different stellar parameters, their estimate of the mass
accretion rate is significantly lower (4 $\cdot\ 10^{-8}\ M_{\odot}$yr-1) than
the 29.77${}^{+34.39}_{-14.59}$ $\cdot\ 10^{-8}\ \mathrm{M_{\odot}yr^{-1}}$)
we derived from our data.
At the center channel we obtain a size of 0.081 au (6.5 R∗) for the emission
region, while the change across the line appears to be mostly within the
significant uncertainties, which are amplified by the low flux and high
visibilities in the wings of the feature. Even so, there is seemingly an
increase in size towards the line edges and especially the over blueshifted
channels, where the half flux radius goes up to 0.127 au (10.2 R∗). In the
channels close to the line center, Br$\gamma$ emission originates from within
the corotation radius at 0.09 au (7.3 R∗). If we were to interpret the two
Gaussian components we fitted to the spectrum physically, this could signify
the presence of an additional emission source at scales beyond the corotation
radius, which is responsible for the broad component, while the narrow
component, which is dominant in the center, traces the emission coming from
the magnetosphere. This is further supported by the fact that within the range
of the narrow component, we first see a decrease in size as would be expected
for the case of magnetospheric accretion. A combination of a disk wind
launched from the inner disk and a system of accretion columns would be a
suitable candidate to explain the radius profile, although the large
uncertainties do not allow as to conclude this definitively.
A contributing wind could also explain the distribution of photocenters across
the line, as the separation between red and blue arm is similar to a
rotational profile which is aligned with the position angle of the NIR
continuum disk at (91 $\pm$ 6)∘. Such an alignment is typical for a disk in
rotation, and a wind launched from such a rotating disk might be expected to
display a similar displacement of photocenters, see e.g. Goto et al. (2012).
Again, the low signal to noise ratio in our differential phase data results in
very large error bars, making any kind of commentary ultimately speculative,
but it is notable that the shift (3.8 R∗) we detect for a photocenter is
significantly larger both in absolute terms and as a fraction of the size of
the Br$\gamma$ emission region than what we find for any of the other targets
in the sample. This implies that there is some influence which distorts the
symmetry of the Br$\gamma$ brightness distribution at larger scales than what
we usually see.

####  5.2.4 DG Tau

Figure 9: Spectrum (top), size (middle) and photocenter shift (bottom)
profiles for DG Tau. Shaded regions indicate uncertainties.

GRAVITY observations of the DG Tau (Fig. 9) Br$\gamma$ flux show a weak
spectral line, with a peak flux ratio of only 1.29, the second lowest in the
sample, and a FWHM of about $\pm$ 224 km/s. The original assumption of DG Tau
as one of the strongest accretors in the sample appears at odds with the small
equivalent width of only 5.79 $\AA$, from which we estimate a mass accretion
rate of just $\mathrm{3.95^{+4.44}_{-2}}$ $\cdot\ 10^{-8}\
\mathrm{M_{\odot}yr^{-1}}$. This puts this observation at the very low end of
a range of previously reported accretion rates, some of which can surpass the
current result by an order of magnitude (see e.g. Iguchi & Itoh (2016). Such a
degree of variation among past results indicates either high variability in
accretion itself, some form of temporary obscuration of the central star and
the Br$\gamma$ emission region, or a combination of similar factors leading to
an attenuated Br$\gamma$ feature.
Fits of the continuum data with a Gaussian ring model have constrained the
inner dust disk PA at (143 $\pm$ 12)∘ (Gravity Collaboration et al., 2021b),
so that we were not able to detect a misalignment between the K band continuum
and the outer disk at 135.4∘ position angle as measured by ALMA at 870 $\mu m$
(Güdel et al., 2018). The semi-major axis of the inner dusty disk appears
almost perpendicular to the reported jet PA (Solf & Boehm, 1993) of 226∘. By
contrast, there appears to be a significant divergence in terms of inclination
between the outer disk at 37∘ and the inner disk at (49 $\pm$ 4)∘. The known
jet inclination of 46∘ (Eislöffel & Mundt, 1998) agrees with the GRAVITY
result, which suggests that the inner disk is aligned with the magnetic field.
For DG Tau, Br$\gamma$ visibility signals tend to be faint, with multiple
baselines showing no recognizable Br$\gamma$ feature at all. The absence of a
clear detection for some of the baselines cannot be directly attributed to a
lack of spatial resolution as there seems to be overall no clear correlation
between strength of signal and baseline length. At a central channel half-flux
radius of 0.143 au (12.9 R∗), the emission region stretches far beyond the
corotation radius and even slightly beyond the NIR continuum HWHM of 0.125 au,
which makes it a unique case within our sample.
Our results appear more extended when compared to those of Eisner et al.
(2014), who examined the Br$\gamma$ feature of DG Tau with the Keck
Interferometer (KI) and found a uniform ring radius of 0.12 au in Br$\gamma$
versus 0.17 au in the NIR continuum for an assumed distance of 140 pc. If
updated to the more recent distance estimate of 125 pc from GAIA DR3, this
would translate into 0.10 au and 0.15 au, respectively. As the Br$\gamma$ line
recorded with KI shows an equivalent width of 6.9 $\AA$ and is thus broadly
similar in strength and width to the one measured with GRAVITY, this
discrepancy is likely the result of differing methods in determining the
region size. Due to their reliance on a single baseline, they fit a uniform
ring model without taking into account possible inclination effects. While the
methodology section of their paper does reference a ”circumstellar emission”
which seems similar in approach to the pure line emission used in our
analysis, it is not clear whether the reported Br$\gamma$ radius is based on
this or the total visibilities. If the latter is the case, this could
certainly cause the emission region to appear significantly more compact.
Differential phase signals are also relatively faint and translate into
photocenter shifts that, even at a maximum magnitude of 0.0125 au, do not
extend significantly beyond the stellar radius of 0.0123 au. The photocenter
profile shows a significant displacement between the central channel
photocenter position and the continuum photocenter at the zero position (i.e.
the center of Fig. 9). When taking into account the closure phase information,
this displacement can be explained in terms of an uncertainty on the relative
location of the continuum photocenter with respect to the star. The position
of the star, which in Fig. 9 is located somewhere within the red circle
determined from the closure phases, is still consistent with the shift
measured for some of the photocenters on a 1$\sigma$ level.
The photocenter profile features a blue arm that is extended along the semi-
major axis of the continuum disk, whereas the displacements across the red
wing appear to be more in line with the axis of the redshifted jet component.
However, the scale of the error bars on the individual points in both wings
prevent us from definitively associating the profile with either of these
alignments. Given the large uncertainties, it is certainly possible to
interpret the profile as essentially Keplerian, possibly tracing the
Br$\gamma$ emission coming from a wind that was launched from the rotating
base of the disk near the inner rim. Not only would such a case be consistent
with the large spatial extent of the Br$\gamma$ emission region, but also with
the velocity profiles obtained with ALMA by Podio et al. (2020) and Güdel et
al. (2018) from the large scale CO and other molecular line emission, which
also indicate blushifted emission along the south-east axis.

####  5.2.5 TW Hya

Figure 10: Spectrum (top), size (middle) and photocenter shift (bottom)
profiles for TW Hya. Shaded regions indicate uncertainties.

The Br$\gamma$ emission line in the GRAVITY spectrum of TW Hya (Fig. 10) shows
an excess of blueshifted emission, which is also seen in a spectrum recorded
with SINFONI (Goto et al., 2012) with a similar (R=4000) spectral resolution
as GRAVITY and further confirmed in high resolution (R=45000) (Sokal et al.,
2018) observations made with IGRINS between 2015 and 2017.
At the center of the feature we find that the Br$\gamma$ emission region
originates in area with 0.029 au (5.16 R∗) half flux radius, which falls off
to 0.016 au (2.85 R∗) in the extreme blue wing and 0.022 au (3.91 R∗) at the
lowest point in the red wing. As such TW Hya is the only source in the sample
that clearly shows a decline in size towards the edges of the line beyond the
uncertainty on the individual points.
We also find that the emission region in all channels is significantly more
compact than the corotation radius of 0.05 au and slightly smaller than the
magnetospheric radius of 0.037 au estimated from stellar parameters and an
assumed magnetic field strength of 1 kG. As most of the observational evidence
actually suggests the magnetic field to be on the order of up to 1.5 kG (e.g,
Bouvier et al. (2007), Donati et al. (2011), Johnstone et al. (2014)), we can
adjust the estimate of the truncation radius to 0.046 au and still find it to
be less extended, but close to the corotation radius, which is expected for
magnetospheric accretion. This agreement and the compact nature of the
Br$\gamma$ half flux region when compared to the corotation radius are making
a compelling case for the idea that Br$\gamma$ emission in TW Hya indeed
orginates from the magnetosphere. This is further supported by the size
profile. The decrease in HWHM would be expected under a magnetospheric
accretion scenario as line components with the highest velocity are likely
tracing gas in free fall towards the stellar surface.
This very important result, presenting the Br$\gamma$ emission region as the
first direct measurement of a spatially resolved magnetosphere, was previously
reported in Gravity Collaboration et al. (2020). Based on the same data, they
found a Br$\gamma$ region size of 0.021 au and thus something even more
compact, although the deviation ultimately comes down to methodology and data
correction. Even though a Gaussian disk was used as a geometric model in both
cases, they performed their analysis on an average of the three central
channels of the line and used previous ALMA measurements to fix their
inclination at the level of the outer disk (7∘, see e.g. Andrews et al.
(2016)), while we refer to the GRAVITY measurements of the NIR continuum to
perform our fit for a slightly more inclined (14∘) disk. Another key
difference is the spectral correction that we applied to align telescope
spectra with each other and with the telluric lines. As the interferometric
quantities are recomputed based on the corrected spectra, the realignment can
have a significant enough impact on the visibilities and differential phases
to account for deviations such as this.
It is notable that the magnetospheric radius at B=1.5 kG extends beyond the
K-band continuum half flux radius, whereas the dust sublimation radius was
estimated to fall between 0.02 and 0.04 au (Gravity Collaboration et al.,
2021b). This might indicate that dust grains exist up to the truncation radius
and that the magnetosphere is playing a part in truncating not just the
gaseous disk, but also the dusty disk, which was also suggested by Eisner et
al. (2006) based on their Keck Interferometer observations of this target.
The photocenter profile across the line is highly compact, concentrating in an
area within 0.01 au (1.78 R∗) around the continuum photocenter. Since the
dispersion of closure phases at Br$\gamma$ wavelengths is small, it is
reasonable to assert that, given the maximum shift of 0.009 au and an average
shift of less than 0.003 for the central channels compared to a stellar radius
of 0.006 au, the Br$\gamma$ photocenter across the line coincides with the
position of the central star. This suggests that most of the different
photocenter velocity components of Br$\gamma$ emission, observed almost pole-
on, are well centered on the stellar position. Little can be inferred directly
from such a distribution. At very low inclination, any centrosymmetric
brightness distribution is likely to lead to a similarly compact profile
centered on zero, i.e. only in combination with more conclusive evidence, as
presented above, does it support the magnetospheric accretion scenario. It is
still noticeable that the distribution, and indeed the differential phase
signals at the individual baselines, appear to show slightly larger shifts in
the red wing of the line. In the case of an inclined magnetosphere with two
accretion columns, the opposite could be naively expected in a pole-on
configuration, as infalling gas at the far side of the star at the highest
blueshifted velocities would be blocked by the star itself, which should lead
to a shift away from the center for the blue components. The absence of a
stronger shift in the blue arm could indicate that the configuration of the
magnetic field in TW Hya is close to an axisymmetric scenario with a small
obliquity angle between dipole and rotational axes. Such a case would be in
line with the relatively small difference in inclination between the inner NIR
continuum disk and the outer disk and also the lack of blueshifted absorption.
However, a previous study by Johnstone et al. (2014) indicates a substantially
tilted dipole, with an obliquity between 10∘ and 40∘, respectively, and it is
unclear whether this still would be consistent with such a scenario.

####  5.2.6 S CrA N

Figure 11: Spectrum (top) and size (bottom) profile for S CrA N. Due to the
poor quality of the phase data we chose not to include an analysis of the
photocenter shift distribution for this object. Shaded regions indicate
uncertainties.

We observe the a similarly asymmetric spectral line as for most of the
previous sources, although for S CrA (Fig. 11) it appears at a higher
blueshift, with a broad high velocity component centered at -273 km/s.
The equivalent width of -8.7 $\AA$ is consistent with measurements of the
Br$\gamma$ line taken in 2008 with NIRSPEC, as presented in Sullivan et al.
(2019), who reported -7.2 $\AA$. However, they also report a significantly
weaker mass accretion rate of 12.8 $\cdot$ 10-8 $\mathrm{M_{\odot}yr^{-1}}$
when compared to the estimate of 27.72${}^{+28.96}_{-13.57}$ $\cdot$ 10-8
$\mathrm{M_{\odot}yr^{-1}}$ we derived from our equivalent width. The NIRSPEC
K band spectrum around the Br$\gamma$ feature is visually consistent in shape
with our GRAVITY observations and shows the same asymmetry, arguably even more
so than the previously published work (Gravity Collaboration et al., 2017b) on
the same GRAVITY data, which was at the time not corrected for the
misalignment between the four telescopes. S CrA was also investigated by Gahm
et al. (2018) with the UVES spectrograph of the VLT, who fitted a
magnetospheric accretion model to their spectroscopic data and found an even
higher mass accretion rate of 50 $\cdot$ 10-8 $\mathrm{M_{\odot}yr^{-1}}$.
The spectral correction we applied also affects the shape and position of the
Br$\gamma$ signals in the visibilities. In Gravity Collaboration et al.
(2017b) their unusual shape and the misalignment between the visibility peaks
at different baselines with the position of the line in the spectrum was
presented as a possible indicator of the presence of different kinetic
signatures in the data. In our spectrally corrected data, we do no longer find
this misalignment for most of the baselines to any significant degree. Except
for the shortest baseline, the peaks are sufficiently aligned with the
spectral line and are well represented by a single peaked Gaussian shape. For
the UT3-UT2 baseline with a length of 41 m, we still find a visibility shape
that is strongly offset towards the blue wing and decreases much more quickly
in the red, imitating the asymmetry in the spectrum. This could be indicative
of a large scale component in the system that is overresolved at the longer
baselines, although this in itself would not explain the lack of signal in the
red wing at this baseline.
At a number of combinations of channels and baselines, the pure line emission
region is completely unresolved, as indicated by the fact that the computation
of the pure line visibility yields results that are consistent with 1 within
the uncertainties for those wavelengths at those baseline lengths. For the
purpose of determining the region size, such points were discarded so that the
fits of the ellipsoid are based purely on those baselines and channels for
which the region is marginally resolved. While this correction is reflected in
the error bars on the region size in each channel, the results in those cases
will also overestimate the true size and should thus be considered as an upper
limit.
With this in mind it is difficult to comment on the trend of the size across
the different velocity channels since for S CrA we find that it is mainly the
outer parts of the wing, and the red wing in particular, which are affected by
this issue. While the size profile shows an increase from 0.071 au (4.4 R∗) at
0 km/s to 0.104 au (6.4 R∗) in the extreme red and 0.094 (5.7 R∗) in the
extreme blue channel, the uncertainties and upper limit nature of those
results must put some doubt on any possible interpretation.
Even at the lowest extent in the central channels, the region extends outside
the range of the corotation radius at 0.04 au (2.4 R∗). While the large
uncertainties in the line wings again suggest that it is problematic to
overstate this case for the higher velocity ranges, we consider it unlikely
that the degree of overestimation is sufficiently large to mask an emission
point of origin within the corotation radius.
Due to the poor quality of the differential phases, we cannot present a
photocenter shift profile as for our other sample objects and have to rely
largely on the visibilities for interpretation. As such, even when taking into
account the degrading quality of the fits in the edges of the line, neither
the absolute size of the emission region nor the size profile across the line
appear characteristic for magnetospheric accretion. It is likely that
additional Br$\gamma$ emitting components are indeed present in the system, as
the peculiar shape of the UT3-UT2 visibility implied. This was in principle
also the argument put forward in Gravity Collaboration et al. (2017b), where
they found a size of the semi-major axis of a Br$\gamma$ Gaussian disk of
(0.54 $\pm$ 0.07) mas at the center of the line, which would translate to
(0.087 $\pm$ 0.011) au. Given the obvious impact of the spectral correction on
the visibilities, this is still very much consistent with the current results,
although in the previous work the size of the region declined throughout the
line rather than reach a minimum close to the line center. S CrA is associated
with a large scale jet (Peterson et al., 2011) and could also be driving an
outflow on smaller scales which contributes to the Br$\gamma$ feature of the
system, which might contribute to the blueshifted excess in the spectrum. The
-273 km/s at the center of the broad, blueshifted component fit the typical
velocities found in small scale jets, see e.g. Takami et al. (2003). There
they also found a positional displacement in H$\alpha$ in S CrA N as part of
their spectro-astrometic survey. While they argue that this displacement at a
position angle of 153∘ is likely caused by the companion, we do note that
baseline UT4-UT3 at a similar PA of 132∘ consistently yields the largest size
for the Br$\gamma$ emission region across all of the spectral channels. Given
the seperation between S CrA N and S CrA S of 1.3 arcseconds, the visibility
measurements would not be contaminated by this visual companion.
However, differentiating between such stellar outflow or a wind being launched
from the rotating inner disk as primary drivers of Br$\gamma$ in S CrA N is
challenging given the lack of robust differential phase data to trace the
kinematics of the gas.

####  5.2.7 DoAr 44

Figure 12: Spectrum (top), size (middle) and photocenter shift (bottom)
profiles for DoAr 44. Shaded regions indicate uncertainties.

DoAr 44 (Fig. 12) exhibits a symmetric line shape with peak line-to-continuum
flux ratio of only 1.19, making it the weakest Br$\gamma$ emitter in the
sample. Similar to the case of S CrA N, the pure line emission region is
completely unresolved at a number of channels for certain baselines, with pure
line visibilities reaching values of 1 within the uncertainties. The so
affected channel-baseline combinations were disregarded in the ellipse fitting
process, as with a geometric Gaussian disk model these can only be associated
with FWHMs of 0. The two outmost channels that meet the selection criterion
are most affected by the reduced number of data points to fit the ellipse to,
leading to disproportionately large error bars compared to the inner 4
channels. Even then, only for 3 of those 4 channels do we resolve the emission
region under all 6 baselines simultaneously, which consequently makes those 3
channels the most robust estimators of the true region size that should be
compared against the corotation radius of 0.05 au.
The emission region has a half flux radius of 0.05 au, which remains
relatively constant between the center channels within the uncertainties. With
Br$\gamma$ emission coming from a region similar in extent to the corotation
radius, it is likely that magnetospheric accretion is a driving mechanism
behind a large part of the total line emission and that a significant portion
of the emission will originate in the accretion columns of the magnetosphere.
At about 5.9 R∗, this size falls into the range of values typically associated
with magnetospheres (Bouvier et al., 2007).
Our results are in agreement with those previously reported by Bouvier et al.
(2020), who published a dedicated paper on a multi-instrument campaign on DoAr
44, including these GRAVITY interferometric observations. There they reported
an upper limit to the emission region size of 0.047 au, which, when taking
into account the impact of our additional spectral correction, is consistent
with what we found.
The differential phase signals are faint, the resulting distribution of
Br$\gamma$ component photocenters hints at a separation of blue and red arms
similar to a possible rotational profile, but the uncertainties prevent us
from showing this conclusively. The maximum shift in the most blueshifted
channel is of a magnitude of 0.016 au (1.9 R∗), while the red and center
channels cluster around an area at about 0.005 au (0.59 R∗) from the continuum
photocenter. The closure phases show a weak offset in the continuum of around
0.9∘ which is still consistent with a shifted continuum photocenter of around
that magnitude, meaning the observed Br$\gamma$ emission could in truth be
centered on the position of the star. This is also in agreement with the
findings of Bouvier et al. (2020), who reported a continuum offset of 50
$\mu$as, which almost exactly corresponds to the offset of the white channel
in Fig. 12.
Despite the relatively compact emission region, the photocenter profile
appears less concentrated than the similarly compact case of TW Hya, which is
likely rooted in the relatively higher inclination (32 $\pm$ 4)∘ of the inner
disk when compared to the pole-on case. The profile appears also to be almost
perfectly aligned with the semi-major axis of the inner disk (137 $\pm$ 4)∘,
which would be typical for a disk in rotation, but is not necessarily
associated with a rotating magnetosphere. This could be explained by an
additional wind component launched from the disk-magnetosphere interface, but
in that case a blue excess in the line shape similar to what was observed in
e.g. AS 353 would be expected. As it stands, the alignment is either
coincidental or the contribution from the wind is too weak to clearly modulate
the line shape.

##  6 Discussion

Based on the results presented in Section 5, the objects of our sample can be
broken down into three distinct categories: For two of them (TW Hya, DoAr 44)
we found an emission region size that is on the order of the corotation radius
or smaller, and also either on the order of, or smaller than, the truncation
radius we estimated from the equivalent width.
Two of the other objects (AS 353, DG Tau) showed emission regions that extend
beyond 10 R∗, far exceeding both the corotation radius and also the range of
truncation radii at magnetic field strengths of 1-3 kG.
The remaining T Tauri stars (both epochs of RU Lup, S CrA N, VV CrA) do not
fit clearly into either of those two groups, featuring Br$\gamma$ emission
regions with HWHMs either larger than the corotation radius, but consistent
with the estimated range of truncation radii, or smaller than the corotation
radius, but larger than even the upper limit of the truncation radius computed
for a 3 kG field.
In this section, we will discuss these results in the context of different
Br$\gamma$ origin mechanisms in T Tauri stars. We begin by first introducing
our work on a radiative transfer model of magnetospheric accretion in order to
establish a baseline for the expected spectral, size and photocenter profiles
for Br$\gamma$ emission originating exclusively from the magnetosphere. We
will then proceed to consider how deviations from these baseline profiles in
the form of additional components of extended Br$\gamma$ emission would be
able to account for those results which are not well explained by the
magnetospheric accretion scenario.

###  6.1 Radiative Transfer Models of Magnetospheric Accretion

Figure 13: An image of the magnetospheric accretion region of an RU Lup-like
system at the Br$\gamma$ peak wavelength, obtained from radiative transfer
modelling with MCFOST. Depicted here are the accretion columns at the line
center, the star and the shock region at high latitude, as well as a
contribution coming from a faint magnetospheric continuum. The model is
axisymmetric, i.e. the magnetic axis is aligned with the stellar rotational
axis, i.e. the vertical axis in the image above.

Figure 14: The Br$\gamma$ spectral line produced by the radiative transfer
models at different inclinations of the magnetosphere. As we move to higher
inclinations, the peak line-to-continuum flux ratio decreases and the line
broadens. The inverse P Cygni profile with a redshifted absorption feature
disappears beyond inclinations of 54∘, as the accretion columns cover less of
the stellar surface and the shock region. The blue- and redshifted velocity
components of the rotating magnetosphere are increasingly separated, leading
to a double peaked profile at high inclinations. Note that the synthetic data
were convolved with a Gaussian and brought to the spectral resolution of
GRAVITY.

Figure 15:  The three columns on the left show the continuum-subtracted
Br$\gamma$ emission region image at three different velocities, while the two
columns on the right show the the photocenter shift profile and the size
profile across all velocity components of the Br$\gamma$ line, for a range of
inclinations between 5∘ and 67∘. The model images correspond to the blue,
white and red channels marked with a bold black outline, and the photocenters
and region sizes in those particular channels are defined by the red dots and
the semi-major axes of the red ellipses, respectively. Only the semi-major
axis of the ellipse is left as a free parameter of the fit, which emulates the
approach taken with the observational data. Fig. 16 shows the corresponding
profiles if inclination and position angle are also free parameters of the
fit.

Figure 16: Model images, photocenter and size profiles at different
inclinations as in Fig. 15, but in this version the position angle and
inclination of the red ellipse are free parameters of the fit. This approach
does not match the data treatment and suffers from a lack of constraints along
the axis of the uv coverage gap, but will find better fitting solutions at
high inclinations. Note that while the size profile is different between the
two approaches, the photocenter profile is not affected by this and identical
to Fig. 15.

Estimates of the corotation and the truncation radii serve as simple criteria
to differentiate between objects for which Br$\gamma$ emission originates on
spatially compact scales, consistent with the magnetospheric accretion
scenario, and those showing an extended emission region, where other
mechanisms are likely to contribute to a significant degree. Beyond this, it
is challenging to assess the implications of the complete set of
interferometric observables without being overly speculative. To provide us
with a more complete idea of the expected signatures, and to aid us in the
interpretation of the line shapes, the size profiles and the photocenter shift
profiles, we additionally investigated the predictions of radiative transfer
models of the magnetospheric accretion region (Tessore et al., in prep). We
explored a number of model images, see Fig. 13, showing the magnetospheric
emission region at different inclinations and different radial velocities
across the Br$\gamma$ line. The model features an axisymmetric magnetosphere,
meaning the magnetic dipole is aligned with the rotational axis of the star.
While it has been shown in the past that a dipole tilt is common in T Tauris
(e.g. Johnstone et al. (2014)), we consider the axisymmetric model still as a
first step to establish broad interferometric trends and their dependencies on
the model parameters.
Each of the image cubes represents the magnetosphere at a specific set of
stellar and magnetospheric properties, including the stellar radius, mass
accretion rate, magnetospheric temperature and the temperature of the shock
region. We set the parameters to emulate one of our sample targets in order to
produce a set of synthetic observables analogous to those obtained from the
GRAVITY observations. We subsequently applied the methodology laid out in
Section A to determine the half-flux radius of a Gaussian disk and photocenter
shift in each channel, so that the synthetic data were treated in the same way
as the observational data.
To complement this, we chose the physical model of the magnetosphere to
represent an RU Lup-like system at R∗=2.5 R⊙ and an intermediate accretion
rate of $\dot{M}$ = 5 $\cdot$ 10-8 $M_{\odot}yr^{-1}$ at a mass of 0.8
$M_{\odot}$ and a rotation period of 3.7 days. The magnetospheric temperature
was set to 7500 K, with a shock region heated to 8000 K and a total
magnetospheric radius of 4.57 R∗ (0.053 au), consisting of a truncation radius
of 3.37 R∗ (0.039 au) and a magnetospheric width of 1.2 R∗ (0.014 au).
The model intensity distribution (Fig. 13) features Br$\gamma$ radiation
coming from the emission line itself as well as a continuum contribution
consisting of a photospheric component, the shock region and a faint
magnetospheric continuum. It does not include an extended continuum component,
such as the inner regions of a dusty disk. While this will affect the direct
comparison between synthetic and observational total observables, the pure
line quantities remain unchanged. As detailed in Appendix C, in depth
knowledge about the exact composition of the continuum is not required in
order to recover the pure Br$\gamma$ line visibilities and differential
phases.
A model spectrum was computed by integrating the image fluxes at each
wavelength. At low inclinations it shows a single peaked profile with a
shallow absorption feature in the red wing, caused by infalling optically
thick gas obscuring part of the stellar surface and the shock region. This
inverse P Cygni profile, which is considered to be a characteristic feature of
hot hydrogen emission lines produced in accretion columns observed at face-on
configurations (Bouvier et al., 2007), is visible in the model up to
inclinations of 54∘ at velocities of 150 km/s up to 200 km/s, see Fig. 14. At
the lowest modelled inclination of 5∘, the emission line extends over a
velocity range of about $\pm$ 150 km/s, whereas beyond 54∘, the profile
becomes double peaked and the range expands to $\pm$ 220 km/s.
When observed close to pole-on, the image of centrosymmetric accretion flows
forms a disk-like structure with an intensity distribution that features rings
of variable brightness, decreasing from the stellar surface to the edge of the
magnetosphere. At such low inclination, it is reasonably well approximated in
principle by a geometric Gaussian disk model, given that the ”hole” left by
the star in the continuum-subtracted images (Fig. 16 and 15, left side panel)
will be unresolved at typical GRAVITY baselines. As we move to higher
inclinations, centrosymmetry is lost and the internal geometry of the
intensity distribution becomes complex, although the outline is still well
demarcated by an inclined ellipse in the center channel. In the blue and red
wings of the line, Br$\gamma$ emission is concentrated in regions away from
the center, so that the synthetic visibilities and consequently the fitted
ellipse are tracing the corresponding substructures within the total
magnetosphere.
For an inclination of 5∘, the model system exhibits a decrease in size from
the center channel, with a half-flux radius of 0.032 au (2.75 R∗), to higher
velocities, where the region can be as compact as 0.016 au (1.38 R∗). At
higher inclinations, this trend becomes essentially inverted, with the center
channel forming a minimum and an increasing HWHM towards the line edges. In
this case illustrated in Fig. 15, the ellipse position angle and inclination
were fixed to be in line with the magnetosphere, following the approach we
took with the GRAVITY data when fixing the ellipse fit at the PA and
inclination of the NIR continuum disk. We revoked this assumption in Fig. 16.
The results indicate that the exact choice of methodology potentially can have
large impact on the recovered trends across the line at the various
inclinations. When all the parameters of the ellipse are left free for
fitting, the size profile at higher inclinations broadly retains the original
decrease in HWHM seen at 5∘ inclination, even if some of the individual
channels show some fluctuations.
The best fit result of the ellipse for the 5∘ inclination suggests that the
half-flux radius of the Gaussian disk underestimates the actual spatial extent
of the full Br$\gamma$ emission region by a factor of about 2/3 at the center
of the line. There the 0.032 au derived from the geometric model fall short in
comparison to a magnetosphere with a defined size of 0.053 au. This
discrepancy appears less pronounced at higher inclinations, although the
comparison is not necessarily meaningful, as the intensity distribution in
those cases is no longer well approximated by a simple geometric model. So
while the ellipse fits based on the Gaussian disk do yield size estimates on
roughly the appropriate order of magnitude, it is clear that a geometric
Gaussian disk or ring model is not well suited to represent the realities of a
system of accretion columns at higher inclinations, even for an axisymmetric
magnetosphere.
The photocenter in each velocity channel was determined from a fit of the
synthetic differential phases to again ensure consistency with the approach
taken with the GRAVITY data. As a consistency check, we also computed the
photocenter directly by calculating the barycenter of the intensity
distribution, and found excellent agreement between the results of the two
methods. The so obtained photocenter shift profiles are highly compact at
inclinations of 5∘ (Fig. 15 and 16, center panel), as at such a low angle the
system of funnel flows is essentially centrosymmetric for the model of an
axisymmetric magnetosphere. The maximum magnitude of shift in such a case is
close to 0.007 au at 5∘ , which corresponds to 0.6 R∗. As the inclination
increases, we see a more extended profile, similar to a disk in Keplerian
rotation, which begins to bend inwards at the edges at inclinations above 54∘.
The global maximum photocenter shift is thus found at moderate inclinations
and remains relatively stable from 37∘ to 54∘ at 0.034 au or 2.9 R∗. At low to
intermediate inclinations, the quasi-rotational profile is aligned at an angle
of about 60∘ relative to the ”North” (i.e. the vertical) axis. In this the
photocenter profile of the rotating magnetosphere is distinct from a rotating
disk, as in the latter case the profile would be aligned with the disk semi-
major axis.

###  6.2 Comparison of model to observational data

Figure 17: Synthetic observables (solid lines) computed from the radiative
transfer images of the magnetosphere, compared to the observational data of RU
Lup from 2021. A Br$\gamma$ spectrum was computed from integrating the image
fluxes at different wavelengths. Synthetic visibility amplitudes and phases
were obtained from Fourier transformation of the images. In order to ensure
comparability with the total observables recorded by GRAVITY, an additional
continuum disk contribution, taking into account parameters like disk size,
disk flux and halo flux, was included through a form of post processing.
The figure shows the synthetic data for disk parameters chosen in accordance
with the results reported in Gravity Collaboration et al. (2021b).

Figure 18: Left: A continuum subtracted model image of the magnetosphere at
26∘ inclination, representing the pure line emission region at the central
channel of the Br$\gamma$ line. The size of the region and its photocenter,
depicted here as ellipse and red dot, were both determined from the synthetic
observables in a process analogous to the way the GRAVITY data were treated.
Center: The model photocenter shift profile across the line, compared to the
results for RU Lup 2021. The model profile is largely aligned along a single
axis, while the observational profile shows a significantly higher degree
complexity.
Right: Variation of emission region size across the line. While the size of
the magnetosphere is an input parameter of the model, the change of the half
flux radii at different velocities still provides a useful baseline to which
the trends in the sample can be compared. The observational profile shows a
region more than twice the size of the model HWHM at the central channel,
growing larger at the line edges, whereas the model predicts a decline in size
in the wings of the line.

Whilst the model cube was set up to be comparable to an RU Lup-like system, it
is still useful to compare the broader trends against the rest of the sample.
The range of inclinations covered by the 7 objects starts at the almost face-
on case of TW Hya and goes up to (49 $\pm$ 4)∘ for DG Tau, with most of the
objects being observed at low inclinations between 20∘ and 30∘. In such cases,
the spectral line is single peaked with a slight asymmetry caused by
redshifted absorption at low inclinations, and a blueshifted line peak at
intermediate inclinations. This behaviour is not consistent with the observed
emission lines from the GRAVITY sample, as clearly visible absorption features
in the red wing or shifts of the peak emission away from the center are never
detected in the observational data. Conversely, the asymmetric broadening of
the base of the emission line, observed in 4 of the 7 targets, is not
reproduced by the model.
The size profile of the model shows that, at those low to moderate
inclinations represented by the majority of sample objects, the HWHM either
decreases towards the line edges, or remains relatively flat. The photocenter
profile is either very clustered within a radius of 1 R∗, if observed face-on,
or rotational and well aligned along a specific axis at an angle to the disk
PA. Both of these trends together are roughly consistent with only the two
weakest accretors. For the other objects we measure, within the error bars,
either flat or increasing sizes at higher velocities. Only the two epochs of
RU Lup exhibit a clear U-shape, such as can be seen in the model. However, in
the model this feature only appears at significantly higher inclinations than
those determined for the RU Lup system.
Most of the objects show photocenter profiles with higher orders of complexity
than a simple alignment along one axis. One notable exception here is AS 353,
for which the photocenter profile agrees with the model prediction in that it
shows both a single axis rotational profile, but also a misalignment between
that photocenter axis and the disk PA. However, for this target the emission
region size is highly extended far beyond typical truncation radii, indicating
that these photocenter profiles are not unique to the case of a rotating
magnetosphere.
Overall the model trends are most accurately replicated by the data for TW
Hya. Both the photocenter clustering and the smaller half-flux radii at higher
velocities correspond very well to what is depicted in Figures 15 for a 5∘
inclination, further strengthening the argument for a Br$\gamma$ emission
region dominated by magnetospheric accretion (Gravity Collaboration et al.,
2020). Still, even for this object, the inverse P Cygni profile that is
predicted by the models is curiously absent from the observation, as is the
case across the entire sample. It is possible that in TW Hya the asymmetric
line shape is rather an expression of the absorption in the red wing rather
than of an excess in the blue wing. The relatively narrow emission line in TW
Hya, which matches up well with the range of $\pm$200 km/s predicted by the
model, could support this notion. We detect a notably larger range of
velocities in emission for the other objects with strong asymmetric features,
where the emission region size indicated spatially extended Br$\gamma$
emission on a scale far beyond the magnetosphere.
DoAr 44 is another candidate for a Br$\gamma$ emission region that is strongly
shaped by magnetospheric accretion, although in this case the comparison
against the model profiles is not as conclusive as for TW Hya. At the higher
inclination of 32∘ of DoAr 44, the size profile derived from the model is more
impacted by the methodology behind the emission region fit, although Fig. 15
shows profiles which are relatively flat for the center channels at 26∘ and
37∘. This would be consistent with what is seen in Fig. 12. The photocenter
profile of DoAr 44 can be interpreted as quasi-rotational, although the low
number of selected channels with sufficient flux and the relatively large
uncertainties on the photocenter positions prevent this from being conclusive.
Whether the GRAVITY spectrum of DoAr 44 shows the expected redshifted
absorption is open to interpretation due to the relatively high dispersion
among the continuum-level data points. At this level of inclination, the
absorption peak would already be significantly less pronounced than in the
pole on case of TW Hya, see Fig. 14, so that the discrepancy between observed
spectral line and model line shape is not large. The model thus seems to again
support the idea that Br$\gamma$ emission in DoAr 44 is indeed largely
stemming from magnetospheric accretion processes. (Bouvier et al., 2020)
VV CrA presents itself as something of an edge case here due to the very large
uncertainties associated with the HWHMs and photocenter offsets. It appears
that at least for the central channels the region is more compact (6.56 R∗)
than the corotation radius (7.25 R∗), yet we also find a magnetospheric radius
between 2.86 and 5.36 R∗ for 1 and 3 kG magnetic fields, respectively. As the
emission region extends further than what we would consider a reasonable upper
limit on the truncation radius, the Br$\gamma$ emission does not seem to trace
the magnetosphere. For the remaining objects, comparisons between images and
observations indicate that Br$\gamma$ emission is not well described by only
the magnetospheric accretion model. We see strong asymmetries in the line
shapes, with moderate to high amounts of excess blueshifted emission, which
are not explained by magnetospheric accretion alone. This is also true for the
size profiles across the line, which appears either flat or U-shaped, with a
minimum at the center. Such profiles can be seen in the model when fixing
inclination and position angle for the ellipse fit (Fig. 15), but generally at
higher inclinations than those measured for our sample objects. The
photocenter profiles equally show a high degree of complexity, with either
multiple alignments in different parts of the line or shift magnitudes far
beyond the magnetosphere. We also must account for the possibility that some
of the discrepancies between models and observations certainly will arise from
the axisymmetric nature of the model, given that non-aligned cases are common
and obliquity angles between stellar dipole and rotational axis are typically
on orders of 10∘ (McGinnis et al., 2020). This would likely affect the
photocenter shift profiles the most, whereas the impact on the size profile is
not trivial to evaluate intuitively.
To facilitate a more detailed comparison between the model of the RU Lup-like
system and the data of the RU Lup observations, we made use of a custom
developed tool which allowed us to derive the spectral and interferometric
observables at a uv coverage set to correspond to the actual observation of RU
Lup in 2021. The tool has the ability to modify the observables in such a way
that they reflect the impact of varying physical and geometric properties of
the system, such as the inclination of the magnetosphere. These modifications
also include the possible addition of a continuum disk contribution to
compensate for the lack of a K-band emitting dusty disk in the original model,
see Appendix C. The full set of parameters that can be manipulated to produce
the spectrum, visibility amplitudes and differential phases are: the common
inclination and position angle of a system of continuum disk and
magnetosphere, the size of the disk and the flux percentages coming from
either the disk or an unresolved halo component.
Fig. 17 shows the synthetic data versus the observational data from 2021. In
this case, the continuum disk size, disk flux, halo flux, position angle and
inclination were set to the values reported in Gravity Collaboration et al.
(2021b). At those parameters, the peak flux agrees with the observational
results, but the line shape is lacking the broad base seen in the RU Lup data.
The absorption feature in the red wing of the model spectrum translates to a
similar dip in the model visibilities, as the obscuration of the stellar
surface will increase the relative weight of flux coming from a more extended
area. The synthetic visibilities reproduce the observational continuum
visibilities, but overestimate the Br$\gamma$ peak, suggesting that the model
is significantly more compact than the emission region that is reflected in
the observational data. This is also reflected in size profile (Fig. 18),
where we also see that the change across the line is essentially inverted for
the model when compared to the observational profile. The synthetic
differential phase exhibit S-shaped signals at all baselines, with two phase
peaks with opposite signs at red- and blueshifted velocities. The amplitude of
these peaks ranges from 1∘ to 5∘, which is an order of magnitude larger than
the signals we detect in the observational data. This translates into a far
more extended photocenter shift profile than the comparatively compact
observed shifts for RU Lup in 2021, and also shows a clear alignment along one
axis. The model here fails to reproduce the curved shape of the observational
photocenter profile. We also see that the S-shapes produced by the model are
essentially inverted in sign when compared to the observational data,
indicating that the sense of rotation of the model is also essentially
inverted with respect to the observed RU Lup system.
We also made an attempt to reproduce the observational data as closely as
possible by fitting the model parameters, and found that we were not able to
produce synthetic data in agreement with flux, visibilities and differential
phases simultaneously. It is possible to reproduce the visibilities quite
accurately, but only at a significantly higher inclination than what was found
for RU Lup, and in such case the line flux would be severely underestimated.
It is possible to find a much better agreement between synthetic and observed
differential phases, mostly by manipulating the position angle of the system,
but only at the expense of the visibilities. We were also not able to
reproduce the broad base of the spectral line through any combination of
parameters.
It is possible that a non-axisymmetric model, featuring a tilted
magnetosphere, will be better suited to resolve some of these discrepancies.
The photocenter shift profile in particular, and its dependency on the
inclination, is likely to be affected by the formation of preferential funnel
flows per hemisphere, as is predicted in such a case. However, even such a
modification to the model will not be able to reconcile the observed large
scale emission in some objects with the compact scale of typical
magnetospheres, nor would it account for broad base of the emission lines or
explain Br$\gamma$ emission coming from beyond the corotation radius.
With this in mind, based on the observational results and the work done on the
model cube, we reaffirm our original classification of our sample objects:
sources for which magnetospheric accretion is the main driver of Br$\gamma$
emission (TW Hya, DoAr 44), those which show spatially extended emission,
likely predominantly originating in a disk or stellar wind (DG Tau, AS 353),
and those that exhibit features not consistent with the model, but with
characteristic sizes still on the order of the 4-7 R∗ range associated with
typical magnetospheres. For the latter group, a mixed contribution to the
total Br$\gamma$ flux coming from multiple emission components is likely
responsible for the detected size and photocenter profiles.

###  6.3 Extended Br$\gamma$ emission: disk and stellar winds

While hydrogen line emission in T Tauri stars is strongly associated with the
magnetic heating mechanisms in the accretion flows, magnetospheric accretion
alone can only account for Br$\gamma$ radiation produced in a compact region
close to the stellar surface. This implies that spatially extended emission
must arise from one or multiple other components in the system. In the context
of Herbig Ae/Be type YSOs, the hot accretion disk itself is readily proposed
as a possible source of Br$\gamma$ emission (e.g. Kraus et al. (2008),
Tambovtseva et al. (2016)), given the hot temperature environments found in
the inner disk of intermediate mass objects. This idea has been extended to
the case of the T Tauri stars in the early studies published by Eisner et al.
(2010) and Eisner et al. (2014), who chose a Keplerian disk model and an
infall/outflow model when attempting to reproduce their interferometric data
and found no significant difference between the two in terms of the quality of
their fit. However, computations of the radial temperature profile in the
inner disk of T Tauri stars (Bertout et al., 1988) suggest disk temperatures,
even at distances of 0.02 au, do not exceed 4000 K and thus fall short of the
8000 - 10000 K range we would expect from a strong Br$\gamma$ emitter. By
contrast, it has been shown that in magnetically driven winds hydrogen gas can
be heated to temperatures of 104 K via ambipolar diffusion very quickly after
being launched (Garcia et al., 2001b). This makes a wind a suitable candidate
for an extended Br$\gamma$ emitting component in the inner disk system. Such a
wind might come in the form of an actual disk wind, launched from the surface
of the inner disk, or originate from the polar region of the star as a stellar
wind (e.g. Ferreira et al. 2006, Matt & Pudritz 2005). Interactions between
stellar and disk winds, which can occurr simultaneously, can also lead to
material being ejected from the interface region between disk and stellar wind
(Zanni & Ferreira, 2013). Ultimately it is unclear whether it would be
possible to distinguish between different types of outflow purely based on
interferometric observations, as wide opening angles may lead to certain
overlap of the possible scales of wind-driven emission. For the purpose and in
the context of this discussion, we will only refer to a wind as a Br$\gamma$
emission component that is spatially extended beyond the magnetosphere and not
differentiate any further. In any case, given that the rotation of the
accretion disk is likely to transfer to and kinematically dominate the launch
region of a disk wind, such a scenario could potentially reconcile the low
temperatures of the accretion disk with the quasi-rotational photocenter
profiles that we found for some of the our objects with extended emission
regions.
Disk and stellar winds as Br$\gamma$ emission sources have in the past also
been discussed for Herbig Ae/Be stars: Kurosawa et al. (2016) examined
spectro-interferometric data taken with VLTI AMBER, the predecessor to
GRAVITY, and modelled the Br$\gamma$ feature and associated interferometric
observables with both a magnetosphere and a disk wind. Garcia Lopez et al.
(2015) similarly explored AMBER data on the Herbig star VV Ser and found that
a disk wind model was able to best reproduce the interferometric data, but did
not match the line shapes, which were better reproduced by a bipolar outflow
model. Weigelt et al. (2011) examined the scenario of a wind launched from a
rotating disk for a Herbig Be star and produced model images for a Br$\gamma$
disk wind with which they were able to reproduce their interferometric
observables. The images they present show changes in the intensity
distribution across the line that appears very similar to a rotating
magnetosphere observed at low inclination, and would likely lead to a
similarly rotational photocenter profile. Wilson et al. (2022) presented a
spectroscopic study on hydrogen emission lines in a number of T Tauri stars,
including DG Tau, which they compared against radiative transfer models of
accretion with an additional stellar wind. Their model was able to reproduce
the characteristics of the H$\alpha$ line, but showed significant divergence
from the observed infrared emission lines, including the Br$\gamma$ feature.
Most objects in the GRAVITY T Tauri sample show emission lines which are
either comparable in width to those derived from the radiative transfer model
for low to moderate inclinations, or are asymmetric to the point that they are
best described by a superposition of two Gaussians. Of the two Gaussians, one
is typically broader and shifted to a higher blueshifted velocity, whereas the
second one is comparatively narrow and centered. The narrow component is
across the sample consistently similar in width to the magnetospheric model
emission line, extending between $\pm$ 200 km/s, whereas the broad component
appears at quite different widths and positions for the different targets, and
can change significantly even for the same object between two epochs (see Fig.
7).
It is possible that the asymmetry is a manifestation of a multi-component
Br$\gamma$ emission region, consisting of a magnetospheric accretion component
and additional outflow components in the form of disk or stellar winds. Such
an interpretation would help to explain the prevalence of the asymmetries and
broad bases of the line shapes throughout the sample, given that the presence
of outflows is generally expected in accreting objects due to their role in
managing angular momentum in the system. This is especially true for objects
with known associated large scale jets, such as AS 353, DG Tau and RU Lup. We
note that AS 353 and RU Lup feature emission beyond the corotation radius as
well as broad high velocity Gaussian components centered between -150 km/s to
-50 km/s, which is similar to what has been reported by Garcia et al. (2001a)
for the high velocity component of an atomic wind close to the star.
In AS 353, Br$\gamma$ radiation was detected up to a HWHM of more than 13 R∗,
i.e. about twice as far out as the corotation radius of the system. The
photocenter profile is quasi-rotational, which would match the idea of a wind
base in rotation. However, it is not clear if the scale of the photocenter
profile is consistent with such an extended wind, as even the most displaced
photocenter exhibits a shift magnitude of less than 2 R∗.
In RU Lup we equally observe Br$\gamma$ emission from outside the corotation
radius, but at 6.3 R∗ in 2018 and 5 R∗ in 2021 the emission region HWHM is
still consistent with typical magnetospherical scales. If there is an
additional wind present in the system, its contribution to the total
Br$\gamma$ output would likely be comparatively weaker than in the case of AS
353. The strong U-shape in the size profile at both epochs could also be
explained in the context of a wind. In the case of a magneto-centrifugally
driven wind, observed at low inclination, an accelerated flow of hydrogen
ejected at a wide opening angle would reach higher velocities at more extended
spatial scales. The complex photocenter profile of RU Lup is also not well
explained by either a rotating wind or a rotating axisymmetric magnetosphere
alone, as both of these scenarios on their own should produce a profile with a
singular axis alignment. Whether a combination of a rotating magneotsphere and
a disk wind could lead to an RU Lup-like profile is a question that will
ultimately require more advanced modelling to answer.
The line shape observed in DG Tau is peculiar insofar as that it is one of
only two objects in the sample for which the blueshifted excess is not at all
detected, despite the fact that it has one of the largest emission region
HWHMs of 12.9 R∗ and is also the only object which shows Br$\gamma$ radiation
originating beyond the NIR continuum HWHM. The line extends almost to a range
of $\pm$ 250 - 300 km/s and is as such broader than predicted by the
magnetospheric accretion model. It is possible that inclination effects play a
part here, but since DG Tau is observed at a similar inclination as AS 353,
this alone would not account for the discrepancy in line shape between the two
objects. Either way, the case of DG Tau casts some doubt on whether the
Gaussian decomposition into broad and narrow components can be
straightforwardly tied to physical emission components. Within the
uncertainties, the photocenter profile can be interpreted as effectively
rotational and well aligned with the disk profile, which also supports the
idea of a wind launched from a rotating base, especially since in this case we
know the photocenter profile mimics the one reported for the large scale