© The Authors 2024

1 Introduction

Top-tier optical/infrared interferometers, such as the Center for High Angular Resolution Astronomy (CHARA) Array or the Very Large Telescope Interferometer (VLTI), have benefitted from the great technical advancements in long baseline interferometry over the past 20 yr to the point where once inaccessible regions close to the surface of a in young stellar object (YSO) can be explored in increasing detail. The current generation of instruments now routinely allows astronomers to probe even the innermost disk regions of relatively faint sources, such as T Tauri stars, at angular scales in the sub-milliarcsond regime. Their improved capabilities present the astrophysical community with opportunities to put long held beliefs about the physical processes that govern the star-disk interface to the test. The past years have seen a growing number of studies focussing on spatially resolved observations of the innermost circumstellar disk YSOs (GRAVITY Collaboration 2023a, 2021a,b, 2020; Setterholm et al. 2018). Theoretical frameworks describing the dynamics of ejection and accretion on such small scales have long been discussed in the context of spectroscopic surveys (e.g. Muzerolle et al. 1998, 2001), but only now does long baseline interferometry provide the means to also directly trace the spatial signatures of those mechanisms.

The investigation of the magnetospheric accretion paradigm, according to which the flow of matter close to magnetically active YSOs should be heavily dominated by the stellar magnetic field, has been an important starting point in this context. For T Tauri stars with their kG order magnetic field strengths in particular, the notion is that matter would be funnelled from the inner edge of a magnetically truncated disk along the magnetic field lines onto the stellar surface at high latitudes, see the detailed treatments in, for example, Bouvier et al. (2007) and Hartmann et al. (2016). The temperatures and densities within these funnel flows would then give rise to higher order hydrogen line emission, as in the form of the Brγ line, which would otherwise be absent in the relatively cool environment of classical T Tauri stars.

Initial observations of T Tauri YSOs such as TW Hya (GRAVITY Collaboration 2020) and DoAr44 (Bouvier et al. 2020) with VLTI GRAVITY confirm the assumption that Brγ emission in these systems could thus be spatially constrained to a condensed region close to the stellar surface. In both cases, the authors relied on the co-rotation criterion to determine whether the spatial scale of the emission region is consistent with the scenario of stable magnetospheric accretion. The criterion dictates that the magnetosphere needs to be truncated within the co-rotation radius r_co for stable accretion columns to form, i.e. r_mag ≤ r_co (Romanova & Owocki 2016). In GRAVITY Collaboration (2023b) we follow up on these earlier works with a comparative study of Brγ line emission in seven T Tauri stars with strong Brγ signals. There the interpretation of the derived interferometric sizes and positional offsets of the Brγ emission region is supplemented by the use of a simple axisymmetric magnetospheric accretion model. This model is used to produce synthetic interferometric observables in order to investigate and identify the characteristic behaviour of the magnetospheric accretion paradigm as seen through the lens of interferometry. As one of the key results, the study finds that, while the weakest accretors in the sample, TW Hya and DoAr 44, are in good agreement with the accretion model and co-rotation criterion, the remaining targets show emission on spatial scales beyond typical magnetospheric radii.

Given the apparent limitations of the axisymmetric magnetospheric accretion scenario to explain many of the spectrointerferometric signatures obtained with GRAVITY, it is obvious that a more complex model is needed to approximate the observations. Firstly, the assumption of axisymmetry in the model could be questioned. It is well established that the magnetic dipole in many magnetically active YSOs is tilted with respect to the stellar rotational axis (Donati et al. 2007; Johnstone et al. 2014; McGinnis et al. 2020), meaning real observed systems typically exhibit non-axisymmetric magnetospheric geometries. However, while a non-zero dipole tilt may very well have some degree of impact on the interferometric signatures, it is not clear how this in itself would resolve the discrepancy between the extended emission region size and the comparatively small co-rotation radius. It is more likely that an additional large scale Brγ emission component, such as a disk wind being launched from the magnetosphere-disk interface, is needed to address this issue.

Such winds, among other outflows, have been thoroughly discussed in the context of young stars for decades due to their connection to large scale jet structures that are observed observed around multiple YSO systems, such as RU Lup or AS 353 (Takami et al. 2003; Whelan et al. 2021), as well as their presumed role in managing excess angular momentum in star formation (Hartmann & Stauffer 1989; Shu et al. 1994; Ferreira et al. 2000; Matt & Pudritz 2005). Winds constitute a strong contributor to mass loss and thus disk dispersal in YSOs as these objects evolve towards the main sequence (Alexander et al. 2014; Tabone et al. 2022). As such, accretion and ejection are fundamentally related and mass ejection rates are typically proportional to the accretion rate. Both are known to decline with the age of the YSO as the system matures out of its premain sequence (PMS) phase (see e.g. Watson et al. 2016). In GRAVITY Collaboration (2023b), we show that the younger systems from among the sample, with higher accretion rates, deviate from the pure magnetospheric accretion scenario more strongly. The idea of a wind component potentially acting as an additional emitter of Brγ radiation in younger systems, then diminishing in relative strength with increasing age, could serve as a possible explanation of the distinction we find among the sample objects.

Previous theoretical works by, for example, Garcia et al. (2001) and Lima et al. (2010) suggest that a magnetically driven wind, launched through magnetohydrodynamical (MHD) processes, could reach temperatures on the order of 10⁴ K close to the disk surface. Such an MHD wind could thus provide a sufficiently hot and dense environment to produce a significant amount of Brγ radiation on spatial scales close to the magnetosphere but larger than the co-rotation radius (Romanova & Owocki 2015). MHD wind models have in the past been tested predominantly against spectroscopic observational data (Wilson et al. 2022; Weber et al. 2020), whereas interferometric studies that incorporate radiative transfer simulations are so far scarce. Works that feature models of wind or accretion based hydrogen emission are largely restricted to Herbig Ae/Be stars (e.g. Weigelt et al. 2011; Kurosawa et al. 2016), although more recent contributions have begun to focus on connecting simulations and interferometry in the context of classical T Tauri stars (Tessore et al. 2023; GRAVITY Collaboration 2023b).

In this paper, we build not only on our own previous work, but expand upon the body of interferometric studies of inflow/outflow mechanisms in general (e.g. Eisner et al. 2010). For the first time, we combine sophisticated modelling efforts of the near-infrared (NIR) signatures of disk winds and magneto-spheric accretion in T Tauri stars with the unparalleled spatial resolution of long-baseline interferometry. We re-investigate the comparison of synthetic spectroscopic and interferometric observables to VLTI GRAVITY observational data for the exemplary case of the young star RU Lup in the light of this comprehensive modelling approach. The goal of the analysis is to determine whether these types of complex wind-magnetospheric accretion hybrid models are suited to address the discrepancies between the simple accretion scenario and the RU Lup data.

Section 2 of this work recapitulates the most important observational results on RU Lup. Section 3 details the use of our models. This includes an introduction to the MCFOST radiative transfer code and the technical choices made in preparation of the simulations, as well as a description the model parameters and of the treatment of the synthetic interferometric data. Section 4 contains the results of our simulations in the form of six exemplary model configurations and their comparison to the RU Lup GRAVITY data. Section 5 gives an overview of the effects of parameter variations in the context of the synthetic interferometric data in more general terms. Section 6 presents a thorough discussion of the most important results.

2 Observational results on RU Lup

RU Lup is a classical T Tauri star, located in the Lupus star forming region at a distance of 157.5 pc (Gaia Collaboration 2020). The system is observed at an almost pole-on configuration with an inclination of only 18.8° according to an interferometric study of the large-scale disk structure in the 1.25 mm continuum with ALMA (Huang et al. 2018). RU Lup has a known large scale outflow component in the form of a jet whose blue-shifted component has been determined to be at a position angle of 229° north to east (Whelan et al. 2021). It is generally considered to be a strong accretor, with mass accretion rates in the literature ranging up to 30 × 10⁻⁸ M_⊙ yr⁻¹ (Siwak et al. 2016). Stock et al. (2022) proposed that the accretion rate of RU Lup, as estimated from an empirical relationship between line luminosity and accretion luminosity for certain tracing spectral lines, can vary by as much as a factor of two over a 15 day period, making accretion in RU Lup potentially also highly variable (Sousa et al., in prep.).

RU Lup was observed with GRAVITY in April 2018 and May 2021, combining the four 8.2 m Unit Telescopes (UTs) of the VLTI to record spectrally dispersed interferometric fringes with 6 baselines. The observations made use of GRAVITY’s single field mode, delivering simultaneously low spectral resolution (R ~ 22) fringe tracker (FT) and high spectral resolution (R ~ 4000) science channel (SC) data of RU Lup across the NIR K-band between 1.98 and 2.4 μ m. The effective baseline lengths of the observation range from 46.48 m (UT3-UT2) up to 128.5 m (UT4-UT1), translating into maximum angular resolutions in the K-band between around 4.6 and 1.7 mas, since the angular resolution scales with the projected baseline length B as $\frac{\lambda }{2B}$. In both 2018 and 2021, the total observation time ran for a little over an hour, preceded by a single calibrator measurement in each case.

In the following we summarise the most important results from the GRAVITY continuum and Brγ line studies on RU Lup (GRAVITY Collaboration 2021b, 2023b). The continuum visibility data for RU Lup was fitted with a multi-component geometrical model (Berger & Segransan 2007). It consists of a point source to account for the star itself, a ring-like region which represents the interferometrically partially resolved region around the inner dusty rim, and an extended component which we refer to as the ‘halo’. The halo component is thought to represent the large scale scattered light emission from the disk. It takes into account K-band continuum flux contributions which appear in the spectrum, but are spatially overresolved at the projected baseline lengths of the telescope configuration. The relevance of the halo and the continuum disk will be further detailed in Sec. 3.4. The ring-like morphology of the compact environmental component was chosen due to the association of the K-band continuum emission with the hot dusty wall at the sublimation radius. The continuum analysis yielded a halo contribution of 12 ± 3% and a ring contribution of 30 ± 10% to the total continuum emission. The fit of the ring region, for which a Gaussian brightness profile was assumed, yielded a size in the form of a half width at half maximum (HWHM), or half flux radius, of 0.21 ± 0.06 au in 2018 and 2021. The inclination of the inner disk continuum was determined to be ${16}_{-8}^{+6}$ degrees for the 2018 data and ${20}_{-8}^{+6}$ degrees for 2021 data. The position angles of the ring were constrained at 99 ± 31 degrees and 101 ± 31 degrees for 2018 and 2021, respectively.

Figure 1 shows the observational results obtained for the Brγ emission region. The Brγ line profile shows a small variation in its equivalent width (−16.82 Å to −13.47 Å) from 2018 to 2021, features an asymmetric shape with blueshifted excess emission, and no visible inverse P Cygni characteristic. The visibility data was fitted using a geometrical Gaussian disk model. As such, we presumed that the brightness distribution of the emission region could be described by a disk-like morphology with a centrally peaked brightness profile, which then falls off as a Gaussian function with increasing distance from the centre. In this manner we obtained spectrally dispersed characteristic sizes, again in the form of a half flux radius, for up to 13 velocity channels across the emission line. Similarly, we derived the positional offset of the barycentre of the brightness distribution (i.e. the photocentre) for the same spectral channels from the continuum-subtracted differential phase data. Both the visibility and differential phase data were corrected for continuum contributions, meaning the resulting quantities describe the pure line Brγ emission region. Hereby we found that the Brγ half flux radii (2018: 6.32 R_*, 2021: 5.01 R_* in the centre channel) were more extended than the co-rotation radius (R_co = 3.31 R_* for P_* = 3.7 d, R_* = 2.6 R_⊙, and M = 0.6 M_⊙, as assumed in GRAVITY Collaboration 2021b and GRAVITY Collaboration 2023b), by a factor of more than 1.5. In addition, the emission region sizes were increasing towards the high velocity channels, which was inconsistent with the centrally peaked size profile predicted by the simple axisymmetric accretion model. The distribution of emission region photocentres across the line showed no clear alignment with the known jet axis or the model profile. The magnitude of the displacements was significantly smaller, by a factor of about 4 in the most extreme channel, than that obtained from the synthetic model data.

These clear deviations from the simple magnetospheric accretion scenario make RU Lup a suitable test case to investigate the effect of an additional Brγ emitting disk wind on the synthetic observables. The low inclination configuration, confirmed by GRAVITY, promises to be particularly advantageous compared to some of the other objects from the GRAVITY T Tauri sample for which we found similar results, given that it more clearly separates the accretion and ejection emission regions when the system is viewed near pole-on. In this work we focus on the comparison of the new radiative transfer simulations to the 2021 dataset, since the fundamental trends in the observable profiles are retained between 2018 and 2021. In both cases, we detect an asymmetric line, a spatially extended emission region with increasing sizes at high velocities and a compact multi-axial distribution of emission region photocentres across the Brγ line. For a more detailed description of the reduction process, as well as a more thorough discussion on the results for both continuum and line emission data, we refer to GRAVITY Collaboration (2021b, 2023b).

Fig. 1 Previous results for CTTS RU Lup, derived from VLTI GRAVITY observations in 2018 and 2021. Shown here are (from top to bottom) the line-to-continuum flux ratio of the Brγ spectral line, the characteristic size of the emission region across the velocity channels of the line, and the corresponding positions of the emission region photocentres.

3 Synthetic data modelling

In the following we introduce the methodology used to compute synthetic interferometric observables from the analytical mag-netospheric accretion and disk wind models we used. We first introduce the radiative transfer code MCFOST, before we give a succinct description of the models and their relevant parameters. This section concludes then with a description of the model observable computation and treatment.

3.1 Radiative transfer code

We made use of MCFOST (Pinte et al. 2006, 2009) to produce intensity maps, such as shown in Fig. 2, in the wavelength range of the Brγ line. MCFOST is a dust and atomic line Monte Carlo and ray-tracing based radiative transfer code. It is capable of computing both line and continuum fluxes in multi-level atomic systems under non-local thermodynamical equlibrium (non-LTE) conditions (Tessore et al. 2021, 2023). The atomic hydrogen model used here includes the ground state and 15 excited bound states, leading to a total of 101 discrete bound-bound transitions and 15 bound-free transitions, of which the latter contribute to the total system continuum emission while the former produce discrete spectral line emission. While this hydrogen model file contains a large number of other prominent hydrogen lines, such as Ha, we only consider line emission produced by the Brγ electronic 7–4 transition at a vacuum wavelength of 2.16612 μm.

The images were computed across 101 wavelengths, centred on the line position, with a channel width of 0.723 Å. We chose to compute the atomic populations and radiative transfer on a spherical grid with a logarithmic distribution of points in the radial direction to maximise the density of cells in the central regions of the image, where the relatively complex brightness distribution of the star-magnetosphere system is located. As the axisymmetric magnetospheric model and the wind model are also symmetric with respect to rotation about the stellar axis, they were computed on a 2D spherical grid in order to save on computation time. By contrast, the non-axisymmetric magnetospheric accretion models required the use of a 3D spherical grid due to the azimuth-dependence of the hemispheric accretion flows under a tilted magnetic dipole. On the 2D grids we set the number of radial and polar angle grid points to 150 each, while on the 3D grids an additional 64 points in azimuth are included. We confirm, based on comparative simulations of an axisymmetric model on both a 2D and 3D grid, that resolution effects due to the smaller number of azimuth points do not affect the results. The extent of the grid constitutes a cut-off for the contributions coming from the outer regions of the wind, so the grid needs to be sufficiently large to include all significant flux, but should not be larger than necessary for the sake of computational time efficiency. To ensure this, the outer edge of the grid was tied to a multiple of the outer radius of the model. The exact number depends on the model used and was chosen by manually increasing the grid size until further increases did not significantly change the resulting line profile any longer. The effective range of grid sizes for the models discussed in Sec. 4 lies between 36 and 47 R_*. The final images were produced for multiple inclinations and, where applicable, azimuth angles. The field of view of the image was also anchored to the outer radius of the model. It was again set to be sufficiently wide to capture the entire flux at any possible inclination, but as small enough to minimise the resulting file sizes. The pixel size was, due to similar considerations as for grid size and field of view, set to a constant 0.03 au per pixel, meaning the number of pixels per image varies with the radial size of the chosen model.

Combining two different model components, such as the accretion region and the disk wind, necessitates a careful definition of both model regions via their geometric parameters as to avoid any kind of density overlap between them. Adding another model to the grid will overwrite existing cell information, meaning that any given cell can only be filled with material from either the accretion column or the wind. Beyond this point, all interactions between the different components come purely in the form of photon propagation along the line of sight through the different regions, by which way they may interact with each other.

Fig. 2 Example of a magnetospheric accretion model image, produced with the MCFOST radiative transfer code. Shown here is the continuum-subtracted Brγ emission region of a purely axisymmetric rotating magnetosphere with no additional emission components. Channel line maps such as this can be used to compute artificial interferometric observables which we are then able to compare to real observational data. The image here specifically shows the line map for the −69 km s−1 velocity channel.

Fig. 3 Schematic depiction of the magnetospheric accretion and disk wind models. Shown here are the geometrical parameters defining each component. These are: the rotational axis Ω, the magnetic dipole axis μ, the magnetic obliquity βma, the wind focal point (S) displacement d, the inner and outer radius Rmi and Rmo, and width δRmag, of the magnetosphere, the inner and outer radius Rwi and Rwo, and width δRwind, of the wind, the dark disk radius Rdd beyond which the disk becomes completely opaque, and the cutoff zcrit at which the isothermal wind instantly reaches its 104 K temperature. The lower hemisphere parts of wind and magnetospheric funnel are omitted in this depiction for the sake of clarity.

3.2 Magnetospheric accretion model

Under the magnetospheric accretion paradigm, the inner disk is magnetically truncated at a certain distance from the central star, at which point gas is funneled along the magnetic field lines and transported onto the stellar surface near the poles. In order to compute the emission produced in such a system with the MCFOST code, we employ the model presented in Hartmann et al. (1994); Kurosawa et al. (2006, 2011), and Lima et al. (2010) to define the hydrogen mass density ρ, the temperature profile T(R), and the poloidal velocity υ_p along the field lines. We refer to those works for a more quantitative description of the model.

In the axisymmetric scenario, the system is completely characterised by any 2D slice of the 3D distributions and can be parameterised on a 2D spherical coordinate grid centred on the stellar position. The fundamental assumption of magnetospheric accretion is that the stream of matter within the magnetosphere follows the geometry of the magnetic dipole field lines. The gas is assumed to have no velocity at the truncation radius and then accretes onto the star under the gravitational pull, which defines then the absolute velocity parallel to the field line rooted at the anchor point R_m. Since the magnetospheric funnel is not infinitely thin, but rather extends over a range of anchor points, we define an inner and an outer magnetospheric radius R_mi and R_mo, respectively, and a width δR so that R_mo = R_mi + δR. The velocity along any field line with an anchor point R_m between these boundaries is fully defined by the free fall motion, depending on the stellar mass and radius. In addition, we assume that the magnetosphere rotates as a rigid body, so that the velocity component in the plane perpendicular to the field lines is determined by the stellar rotational period, leading to line broadening effects at non-zero inclinations. The velocity field and funnel geometry also define a shock region on the stellar surface, where a certain amount of accreting material reaches the star per unit time and releases its energy. The density along the accretion funnels is normalised to ensure consistency between this local mass flux onto the star and the global mass accretion rate ${\dot{M}}_{\mathrm{a}\mathrm{c}\mathrm{c}}$. The shock region itself is another source of continuum emission, acting as an additional black body with a temperature determined by the energy released by the infalling material.

With regards to the temperature profile in the funnel flows, we follow Kurosawa et al. (2006) and Lima et al. (2010) by adapting the Hartmann et al. (1994) temperature profile based on a volumetric heating rate $\propto \frac{1}{r^{-3}}$. The temperature along the funnel is then computed based on an assumed energy balance between the radiative cooling rates presented in Hartmann et al. (1982) and the volumetric heating rate. The profile is normalised to a maximum temperature along the funnel T_mag, which we set as a parameter in our model.

For the non-axisymmetric scenario depicted in Fig. 3, we introduce an additional free parameter: the dipole tilt angle β_ma relative to the direction of the stellar rotational axis. We refer to this angle as the ‘magnetic obliquity’ of the system. The dipole tilt also leads to a toroidal component of the magnetic field, which interacts with the velocity field in the plane perpendicular to the field lines by adding an additional toroidal velocity component. For a more thorough description of the density profile and magnetic field components in the non-axisymmetric case, we refer to Tessore et al. (2023) and Mahdavi & Kenyon (1998). In summary, our magnetospheric accretion model is fully characterised by the stellar parameters R_*, M_*, P_*, the mass accretion rate ${\dot{M}}_{\mathrm{a}\mathrm{c}\mathrm{c}}$, the inner anchor point R_mi and width δR_mag of the accretion columns, and the obliquity β_ma.

3.3 Disk wind model

In order to ensure a high degree of flexibility when adjusting the outflow component in the system, we chose to adapt the kinematic disk wind model described by Knigge et al. (1995). Their model follows the basic principle of the magneto-centrifugal disk wind as proposed by Blandford & Payne (1982), which features a mass outflow arising from a disk in Keplerian rotation along a range of open field lines anchored at the disk midplane. It is fundamentally a parametric description of a disk wind, designed to allow for straightforward manipulation of the parameter space to be explored, whilst still approximately reproducing the attributes of a proper, self-consistently treated, MHD wind model. We refrain from reproducing the exact equations here, as they have been introduced in great detail not only in the original publication, but also other works such as Kurosawa et al. (2006) and Kurosawa et al. (2016).

The geometry of the model is biconical, featuring gas streams which follow a set of magnetic field lines lying on conically shaped surfaces, see Fig. 3. This structure is well described by three geometric parameters, which come in the form of the focal point distance d, as well as an inner and outer wind boundary radius R_wi and R_wo. The focal point (S in Fig. 3) defines the origin of the magnetic field lines for one hemisphere, which is vertically displaced from the midplane along the disk axis by the focal point distance d. Adjusting this parameter allows us to effectively set the angle of the outflow and the degree of collimation between the wind stream lines. The inner and outer radii of the wind then give the radial distance from the star at which the closest and furthest field lines intersect the midplane, thus defining the effective wind region in our system as the set of conical surfaces between these two radii.

We further extend these parameters by an additional quantity z_crit, which describes the critical height z above the disk at which the wind reaches its final constant temperature of 10⁴ K. It is depicted in Fig. 3 as the distance between the lower edge of the purple shaded wind region and midplane. Since our implementation treats the wind as isothermal, the critical height acts as a vertical lower end cutoff, effectively defining a starting height for the Brγ emitting component. This approach, which essentially approximates a temperature profile along the vertical axis with little heating at first and then a very sharp rise in temperature close to z_crit, was chosen to simplify the question of the temperature structure in the gas streams. Below z_crit, the gas is assumed to be cold and effectively transparent at NIR K-band wavelengths. The isothermal temperature of 10⁴ K was chosen due to practical considerations, as the Brγ wind emission is significantly reduced below this threshold at typical wind densities, and in particular drops off rapidly below 9000 K. The range at which the wind temperature could be reasonably treated as a free parameter of the model is thus narrow to the point that we only consider this temperature.

The velocity profile of the wind can, as previously for the magnetosphere, be separated into a poloidal component, which defines the velocity along the open field line, and a perpendicular component in the disk plane. The latter is largely driven by the fact that any field line, at its point of emergence from the disk surface, is effectively in Keplerian rotation about the stellar axis. Above the midplane, the rotational velocity of the gas stream deviates from a strict Keplerian profile and decreases with height and distance from the rotational axis. In order to compute the local rotational component of the velocity, we follow Knigge et al. (1995) by assuming that the specific angular momentum with respect to the z-axis is conserved for any one stream line.

The poloidal velocity in the direction of the gas stream is defined by a radial velocity law, controlled by the exponential parameter β_wind. This quantity effectively controls how quickly the wind reaches its terminal velocity as a function of radial distance from the centre, as well as of horizontal distance of the stream line anchor point at the midplane. The terminal velocity itself is defined as a factor f of the local escape velocity at the point of emergence. In our models, we again follow Knigge et al. (1995) and keep f = 2 constant.

The gas density profile of the wind is computed from the local mass loss rate per unit surface area and the geometrical configuration of the wind. The local mass loss rate $\dot{m}$ itself is proportional to the temperature profile of the disk: $$\dot{m}(R)\propto T(R{)}^{4\alpha }.$$(1)

The exact nature of the temperature profile is not relevant in this context, as the continuum disk is not part of the model and the gaseous part is always below the cutoff height z_crit. The α parameter on its own, however, allows us to adjust the local mass loss, by which we can influence the radial distribution of the wind emission to some degree. The exact proportionality between horizontal distance and local mass loss is determined from normalising the profile so that integrating the local mass loss rates of all the wind stream lines between the inner and outer radius adds up to the global mass loss rate ${\dot{M}}_{\text{loss }}$. The density distribution is then fully described by the wind geometry, the velocity field, the global mass loss rate and the alpha parameter.

In addition to the cutoff height z_crit, we lastly introduce a second new parameter that is not part of the original model description in Knigge et al. (1995). The dark disk radius R_dd defines the distance from the central star at which the midplane of the cell grid becomes completely opaque to all photons, effectively creating a thin black layer that blocks all radiation from the other side of the disk. This ‘dark disk’ extends from R_dd to the edge of the cell grid, while any cell within the radius is completely transparent to all radiation. This again is very much a simplified approach to approximate the real influence of optically thick dust at certain distances in the midplane, but is sufficient to allow us to modulate the proportion of radiation from the back side of the disk that is going to contribute to our profiles.

We finally note that the disk wind model is axisymmetric and thus effectively two-dimensional. While it is technically possible to combine a 2D wind with a 3D non-axsisymmetric magnetospheric accretion model, we point out that such an approach would ignore the interactions between wind and accretion that arise in numerical treatments of such combined systems. In such a scenario, we would expect an azimuth dependence in the density of the wind and thus also an azimuth dependence of z_crit, which we do not consider in our analytical approach. In summary, we define the disk wind model component by setting the stellar parameters R_*, M_*, P_*, the global mass loss rate ${\dot{M}}_{\text{loss }}$, the exponent regulating the local mass loss per unit area α, the exponent regulating the radial velocity law β_wind, the inner radius R_wi and width of the wind δR_wind, the dark disk radius R_dd, the focal point displacement d, and the temperature cutoff z_crit.

3.4 Synthetic interferometric observables

The channel maps computed with MCFOST were used to extract synthetic spectral and interferometric data. The model spectral data can be computed straightforwardly by integrating the individual pixels of the brightness distribution over the entire image at each wavelength. Doing so yields a spectrum at the previously defined spectral sampling of the model, which contains the Brγ line emission of the model components and a level of continuum emission defined by the stellar component. If the parameters set the gas in either accretion or wind component to be particularly hot and dense, this may add an additional smaller gaseous continuum contribution.

As the goal is to compare the synthetic to the observational data set obtained with GRAVITY for RU Lup in 2021, the synthetic spectrum needs to be degraded to the correct spectral resolution. This is achieved via convolution with a Gaussian kernel and interpolation on the observational wavelength grid. After normalising the spectrum to the continuum level of the image, this total line-to-total continuum flux ratio $F_{L/C}^{I\mathrm{m}}$ in each channel is then modified to account for contributions coming from the dusty disk and the overresolved halo component to the total continuum which are present in the GRAVITY data but which we do not include in the model image: $$F_{L/C}^{\mathrm{O}\mathrm{b}\mathrm{s}}=\frac{F_{L/C}^{\mathrm{I}\mathrm{m}}+F_{\text{IREx }}}{1+F_{\text{IREx }}},$$(2)

where $F_{\text{IREx }}=\frac{F_{\text{Disk }}+F_{\text{Halo }}}{F_{\ast }}$ is the infrared excess caused by the interferometrically partially resolved component F_Disk and the overresolved halo component F_Halo. For the RU Lup 2021 data, we set the fraction of disk flux to the total continuum to 30% and the fraction of halo flux to 12% in accordance with the results reported in Sec. 2.

The synthetic interferometric observables can be obtained from the Fourier transform of the image. A set of baselines mimicking the configuration of the GRAVITY UTs at the time of the observation were chosen to evaluate the Fourier integral at the corresponding points in the uv plane. In this manner we computed the visibilities and phases in each velocity channel, which were then also spectrally degraded akin to the treatment of the spectrum. We derived the sizes of the Brγ emission region per channel from the synthetic visibilities in the same manner as from the observational data by removing the influence of the continuum first: $$V_{\text{Line }}=\frac{F_{L/C}^{I\mathrm{m}}{V}_{\mathrm{t}\mathrm{o}\mathrm{t}}^{I\mathrm{m}}-V_{\mathrm{C}\mathrm{o}\mathrm{n}\mathrm{t}}^{I\mathrm{m}}}{F_{L/C}^{I\mathrm{m}}}.$$(3)

Here V_Line is the pure line visibility, i.e. the part of the visibility from which the continuum contribution has been removed and which can be directly attributed to the line emission region. As Eq. (3) indicates, V_Line does not require knowledge of the observational line-to-continuum ratio (Eq. (2)) and the associated continuum contributions. The modification based on disk and halo contribution is only strictly necessary for the comparison of synthetic and observational spectral data.

It is possible to relate the pure line visibilities to a characteristic angular size by presuming a certain morphology of the emission region brightness distribution, as was described in Sec. 2. In this case, we follow the treatment of the Brγ GRAVITY data by employing a simple Gaussian disk geometrical model to derive the characteristic size as the half width at half maximum, or half flux radius, of the radial brightness profile. Position angle and inclination of these centrosymmetric brightness distributions can be derived from such a model fit by taking into account the difference in size at different baseline angles. However, while this approach remains valid for disk-like structures, Tessore et al. (2023) and GRAVITY Collaboration (2023b) show that for the magnetospheric region close to the star the result of such a fit is difficult to interpret and does not immediately correspond to the physical inclination angle of the larger inner disk region. To remain consistent with our previous work and to ensure the best possible comparability, we remain with the same approach as before and fix the inclination of the system at the value obtained from UT data continuum fits, i = 20°. We compute the observables for a 90° position angle, as the uncertainty on the PA measurement is relatively high due to the low inclination of the system and the effect on the derived sizes is minimal. The stellar rotational axis in the model is then aligned with the north axis of the image and the disk would be aligned with the east axis.

The synthetic differential phases were computed by subtracting the continuum phase from the total line phase of the image. As the continuum is flat across the line, the continuum phase was taken as the average phase in the first and last image channel. The continuum contribution was then removed to obtain the pure line differential phases, which is related to the angular offset on the sky plane via (Le Bouquin et al. 2009) $$p_i=-\frac{{\mathrm{\Phi }}_i\lambda }{2\pi B_i},$$(4)

where p_i is the magnitude of the offset along the i-th baseline as determined from the projection of the photocentre shift vector onto the baseline vector, Φ_i is the differential phase measured with that baseline in the λ channel, and B_i is the length of the baseline. The so computed angular offsets are relative in nature, they describe the positional shift between the barycentre of the Brγ emission region relative to the barycentre of the underlying continuum brightness distribution. If the continuum region is centrosymmetric with regards to its emission, then the photocentre coincides with the position of the star and the differential phase yields the photocentre offset from the stellar position. While this is generally the case in the models, the observational data can potentially feature an additional continuum offset, caused by, for example, a local density fluctuation in the dust of the inner rim. The magnitude and direction of the total photocentre shift vector can be determined by deprojecting the individual shifts measured along the six baselines of GRAVITY. For a more detailed description of this process, its technical steps and the derivation of the pure line quantity expressions, we refer to GRAVITY Collaboration (2023b).

4 Comparison to GRAVITY RU Lup data

In this section, we present a number of different models in an attempt to recreate the observational profiles obtained from the GRAVITY data for RU Lup from 2021. As relevant observables we consider here the continuum-normalised line profile, the characteristic emission region size obtained by fitting a geometric Gaussian disk model to the visibilities, and the photocentre shift of the emission region as reconstructed from the differential phases.

Exploring the possible model variations is challenging given the large potential parameter space and the computational demands of MCFOST. Even a basic large scale model grid search over the entire parameter range is currently not possible within the technical capacities available and would require a substantial additional optimisation effort that is beyond the scope of this work. Instead, we utilised a mixed approach. We first considered the effect of isolated individual parameter variations from a common reference model to identify those quantities with the strongest response. This analysis indicates that there is a clear distinction between parameters which influence primarily the line flux, but have minor impact on the interferometric quantities, and those that strongly affect also the characteristic size and photocentre shift. We refer to Appendix A for a more thorough discussion of those results. Second, we employed a semi-manual fitting routine which involves running partial model grids for which we varied only a subset of the total model parameters simultaneously within a manually defined bracket of possible values. The best fitting models from those partial grids were then determined by χ² minimisation. This process was then iterated on by computing a new partial grid with different parameter variations to check if the χ² could be further improved.

Given that the instrumental error on the data points computed by the GRAVITY pipeline lie on the order of less than 2%, even the reduced χ² is typically very large. We relied on these quantities exclusively as relative indicators in order to rank different parameter combinations and models against each other and did not derive any statement on the global statistical significance of the models from them. With the simplified nature of our models and the limitations on the parameter exploration in mind, we did not necessarily expect to find a model that would be a full quantitative match of the observational data. Instead, the idea was to determine whether the main trends derived from the observation could be reproduced in principle on an acceptable level and to this end we prioritised achieving a good visual fit over a more quantitative approach.

In this section, we describe in more detail two variants each of a pure magnetospheric accretion model, a pure disk wind model, and of a hybrid combination of both models. Table 1 summarises the parameters of the different models.

4.1 Stellar parameters

The photospheric contribution to the continuum emission was modelled as a blackbody which radiates at an effective temperature T_* over a surface area defined by the stellar radius R_*. In addition, the stellar mass M_* and stellar rotational period P_* affect the velocities in the magnetospheric funnel flows and wind outflows. While in GRAVITY Collaboration (2023b) we assume those parameters to be relatively well defined for the test case of RU Lup, a more thorough review of the literature shows a significant spread of reported figures for this object. Consequently, we chose to adapt a broader range of possible parameter values, treating them effectively as a semi-free parameter in our model exploration. The values adopted for the six models presented in this section are listed in Table 1.

4.2 Scenario 1: Pure magnetospheric accretion model

As a starting point, we revisit the case of the simple axisymmetric magnetospheric accretion model from GRAVITY Collaboration (2023b). In this work we explore a significantly expanded parameter space for this scenario, resulting in the model configuration depicted in the left column of Fig. 4 and referred to as MA I in Table 1. The normalised line flux is single peaked, with a centre channel line-to-continuum flux ratio of 1.76, compared to the observational peak flux ratio of 1.83. The line width between model and observation is comparable at the 50% flux mark, with a model full width at half maximum (FWHM) of 168 km s⁻¹ versus the observational profile at 187 km s⁻¹. At the base of the line, where the observational asymmetry is most pronounced, we note a larger disparity, with a width at 10% flux (W10%) of 320 km s⁻¹ in the model versus a W10% of 420 km s⁻¹ in the GRAVITY data. Visually this difference is clearly localised in the blue wing of the line, where the model profile lacks the excess amount of blueshifted Brγ emission that we detect in the observation. By contrast, the red wing of the line is well reproduced by the model.

In the centre channel, the fit of the Gaussian disk geometric model yields a characteristic size of 0.055 au. As such, the half flux radius sits at half the outer magnetospheric radius of 9 R_*, which at R_* = 2.5 R_⊙ translates to 0.11 au. This result is again comparable to the observational centre channel half flux radius of 0.061 au, although still beyond the estimated 1σ fit uncertainty of 0.001 au. The increase in characteristic size in the wings, however, which we see in the observational data, is not recovered by this accretion model. From the observation we derive a maximum size of 0.1 au at the extreme blue end of the wavelength range, while the model profile is centrally peaked and drops off to 0.045 au in the corresponding channel.

In the photocentre profile, the offset magnitude ranges up to ± 0.022 au, which is by about a factor of two larger than the largest photocentre shift we obtain observationally. The distribution of photocentres across the line has a similar ‘crescent’-like shape, where the offsets at the most extreme velocities become smaller in magnitude again, but compared to the GRAVITY data the high velocity channels are still mostly aligned along a similar axis as the low velocity channels. The alignment of the overall model photocentre distribution appears almost perpendicular to the observational one, but the large uncertainty of the PA estimate obtained with GRAVITY (± 30°) leaves some ambiguity about the relative orientation.

Our second magnetospheric accretion model, referred to as MA II in Table 1 and depicted in the right column of Fig. 4 for an azimuth angle of −90°, addresses the question concerning the effect of a tilted magnetic dipole on the synthetic observables. Most of the fundamental parameters are unchanged with respect to MA I, but the model variant now features an obliquity of 30°. The change results in an overall drop in line flux, with the peak flux reduced to 1.71 and the FWHM to 152 km s⁻¹. At an azimuth angle of −90 degrees, this predominantly affects the red wing of the line profile. For the same angle, the centrally peaked shape of the size profile is flattened to a relatively constant size of 0.051 au across the line. A similar effect is seen for the photocentre profile, where the crescent arch is also compressed, although the maximum shift magnitude remains close to 0.022 au.

The introduction of the dipole tilt breaks the axisymmetric nature of MA I. The image and resulting synthetic observables now depend on a selected azimuth angle, which effectively traces the rotation of the magnetosphere. At −90° in azimuth, the observer’s line of sight is aligned with the two hemispheric funnel flows at non-zero inclinations. In Fig. 5 we concentrate on two specific examples at ± 45°, which in conjunction with the − 90° plot from Fig. 4 cover the relevant potential differences, while in Appendix B we present the full set of images for five different azimuth angles between ± 90°. From these figures it is clear that even qualitatively the nature of the impact of a dipole tilt depends on the azimuth. The localisation of the variation in flux ratio between axisymmetric and non-axisymmetric model changes between the red and blue wing as we move from negative to positive azimuth angles. The flattening of the size profile only appears at negative azimuths, while the positive orientations retain the centrally peaked shape. At +45° azimuth, the photocentre distribution is more strongly aligned with the disk axis, which signifies a shift of about 30° towards east relative to −90° azimuth and also relative to the axisymmetric case.

While the dipole tilt can clearly affect the overall shape of the observable profiles, its impact on the sizes and photocentre shift magnitudes at this level of inclination is even in the best case on the order of 10% or less. In itself, this is not sufficient to fully resolve the limitations of the pure magnetospheric model in matching the GRAVITY data, as the discrepancies between model and data described for the non-axsisymmetric scenario largely remain the same.

Fig. 4 Synthetic observables of the axisymmetric magnetospheric accretion model MA I (left) and the non-axisymmetric accretion model MA II (right) when viewed at an azimuth angle of −90°. From top to bottom: centre channel image, line-to-continuum flux ratio, characteristic size, and photocentre shift per velocity channel of the Brγ emission region. The ellipse depicted in the image signifies the half flux radius of the geometric Gaussian disk model used to derive the characteristic size. In the depiction of the photocenter shifts, the stellar rotational axis, not the north axis, is aligned with the vertical coordinate axis. The observational photocentres were rotated by −9° with respect to Fig. 1 to put them in the same frame of reference.

4.3 Scenario 2: Pure disk wind model

We further present two distinct pure disk wind models in order to explore the question whether the RU Lup GRAVITY data might be better described by a disk wind rather than the magnetospheric accretion scenario. Both models, referred to as DW I and DW II in Table 1 and depicted in Fig. 6, are mainly differentiated by the wind velocity, as represented by a difference in the β_wind parameter of 0.62 (DW I) to 0.05 (DW II). In practical terms, the differing wind velocities result in the line profile of DW I being single peaked, whereas the Brγ line in DW II appears strongly double peaked.

The DW I model manages to replicate the central channel line flux ratio almost exactly at 1.83, although the line is more narrow with an FWHM of only 152 km s⁻¹ and a W10% of 231 km s⁻¹. The width at the base of the line in particular is then underestimated by a factor of almost two when compared to the W10% of 420 km s⁻¹ of the observational data.

The size profile follows the centrally peaked shape of the Brγ line itself, although the central half flux radius of the model overestimates the observational data by again a factor of two (0.121 au compared to 0.061 au, respectively). The large drop in size in the bluest channel in this case is a result of the low line flux at high velocities, which can threaten the integrity of the observable computation as the extraction of the pure line quantities requires division by a (F_L/C − 1) term.

The jumps in the photocentre positions at the most extreme velocities are likely equally caused by this effect, although we note that even at the centre channels, the shift magnitude is on the order of 0.025 au and thus almost three times as large as the most extreme observational offsets. At high velocities this discrepancy can grow to a factor of five, leaving the DW I profile much more extended than the observational profile.

The main purpose of the DW II model is to illustrate that the disk wind model can in principle reproduce the rising sizes at high velocities. The parameters of DW II were specifically selected to emulate the characteristic size profile of RU Lup as much as possible, which in this case comes at the expense of the line profile fit. The increase in size from a centre channel half flux radius of 0.07 au to a maximum size of 0.13 in the blue wing requires the line profile to be double peaked, which does not agree with the observational data. Additionally, while the centre channel sizes are close to the observational profile, those obtained in the wings of the profile are then larger by about 20%. The entire photocentre distribution is rotated by about 45° with respect to DW I, where the photocentres were mostly aligned parallel to the disk axis. Given that the shift magnitudes are comparable between both disk wind models, we find a similarly large disparity between model and observation for DW II as for DW I.

Fig. 5 Synthetic observables of the non-axisymmetric accretion model MA II when viewed at azimuth angles of −45° (left) and +45° (right). From top to bottom: centre channel image, line-to-continuum flux ratio, characteristic size, and photocentre shift per velocity channel of the Brγ emission region. The ellipse depicted in the image signifies the half flux radius of the geometric Gaussian disk model used to derive the characteristic size. In the depiction of the photocenter shifts, the stellar rotational axis, not the north axis, is aligned with the vertical coordinate axis. The observational photocentres were rotated by −9° with respect to Fig. 1 to put them in the same frame of reference.

Fig. 6 Synthetic observables of the slow disk wind model DW I (left) and the fast disk wind model DW II (right). From top to bottom: centre channel image, line-to-continuum flux ratio, characteristic size, and photocentre shift per velocity channel of the Brγ emission region. The ellipse depicted in the image signifies the half flux radius of the geometric Gaussian disk model used to derive the characteristic size. In the depiction of the photocenter shifts, the stellar rotational axis, not the north axis, is aligned with the vertical coordinate axis. The observational photocentres were rotated by −9° with respect to Fig. 1 to put them in the same frame of reference.

4.4 Scenario 3: Hybrid models

Finally, we introduce a pair of hybrid scenarios in which we combine the disk wind and magnetospheric accretion region into a single common model brightness distribution. While both models, which in Table 1 are referred to as HM I and HM II and which are depicted in Figs. 7 and 8, differ in a number of parameters, we first and foremost distinguish between them primarily based on the co-rotation criterion. For the model HM I, which is characterised by a cooler and more compact magnetospheric configuration than HM II, the size, rotational period and stellar parameters were selected to allow the magnetosphere to be truncated within the co-rotation radius. For HM II, we disregard this requirement and attempt to achieve the best fit to the data regardless of whether the co-rotation criterion is respected.

For HM I we find a centre channel flux ratio of 1.85 along with a FWHM of 180 km s⁻¹ and a W10% of 340 km s⁻¹. Flux ratio and FWHM deviate from the observational results by less than 5%, although we see that W10% is still underestimated. However, visually the disparity in the line base can no longer be clearly attributed to the lack of blue excess emission, but is also caused by the presence of an inverse P Cygni feature in the red wing which we do not detected observationally. On the contrary, the line flux in the blue wing is reproduced well up to velocities of about −150 km s⁻¹ before it diverges significantly from the RU Lup data.

The centre channel characteristic size of 0.047 underestimates the observational result by less than 25%, while the model still reproduces the increase in size towards the edges of the line. In the red wing the observational result is well matched at 0.09 au, although in the blue wing we do note a drop in size below −150 km s⁻¹, which does not agree with the data.

The largest photocentre shifts at high velocities show a displacement magnitude of 0.045 au, which is still four times larger than the largest shifts observed for RU Lup. However, only the highest velocities show such a strong photocentre offset, whereas above −150 km s⁻¹ the largest shift is already only about half as large (0.023 au). There is also a clearly visible change in the alignment of the profile at the same velocities. While the high velocity channels are aligned at an angle of about 45° relative to the disk axis, the remaining channels show photocentre offsets which are almost perfectly oriented along the stellar rotational axis.

The left column of Fig. 8 depicts the decomposition of the three model observables of HM I into their the constituent components. In the blue wing the wind becomes the dominant contributor to the line flux below −150 km s⁻¹, whereas in the red wing the influence is more mixed, as the inverse P Cygni feature of the magnetosphere still appears clearly in the combined line profile. The decomposition of the interferometric quantities shows that the shape in the wings largely follows the wind, but the combined region characteristic size is adjusted downwards by the magnetospheric component. Equally, we see that the alignment of the photocentres at low velocities results from a superposition of differently aligned shift vectors between wind and magnetosphere. The high velocity channels follow the 45° alignment of the disk wind profile, although the magnetospheric influence again reduced the magnitude of the shift compared to the pure wind component.

For HM II, the centre channel flux ratio and FWHM (1.77 and 209 km s⁻¹, respectively) fit the observation similarly well as they did for HM I, but the W10% is now increased to almost 400 km s⁻¹, which compares more favourably to the observational W10% of 420 km s⁻¹. The line shape lacks the P Cygni characteristic of HM I, which improves the fit in the red wing to give an almost perfect match to the data points at velocities above +150 km s⁻¹. The flux ratio in the blue wing remains close to the observational ratio down to almost −200 km s⁻¹ in this model.

The deviation in centre channel size is reduced to about 7% and the sharp drop off in half flux radius at −150 km s⁻¹ is much attenuated to provide an overall visually clearly improved fit to the data. In the photocentre shift profile of HM II, the largest offsets outside the extreme blue spectral channel are reduced in magnitude. In the red, the maximum offset is now 0.023 au, compared to the more than 0.06 au in the same channel of HM I. By contrast, the low velocity channels appear to feature photocentre shifts that are larger by up 20% in magnitude with respect to their corresponding offsets in HM I and the axis of alignment does no longer coincide with the stellar rotational axis but is rather at a −45° degree angle relative to the disk axis.

From the decomposition of HM II shown in the right column of Fig. 8 we see that the magnetospheric component in itself lacks the P Cygni characteristic and that it produces a line that is about 25% broader when compared to HM I, while the relative contribution of the wind component is similar. The characteristic sizes of both wind and magnetospheric component are closer to the combined profile than before and the dominance of the magnetosphere in the centre channels is even greater as now up to three channel total characteristic sizes are almost fully defined by the magnetospheric component. This is again equally true for the low velocity channels of the photocentre shift profile, where the decomposition shows that the alignment of the total profile is also close to identical to the alignment of the magnetospheric photocentre offsets. A plot showing the mass density, temperature and velocity fields of HM II is included in Appendix C.

5 Discussion

In Sec. 4 we present results for a total of six different model combinations, featuring two variants each of a pure magnetospheric accretion model, a pure parametric disk wind model, and a hybrid combination of both. The synthetic spectro-interferometric observables derived from these model images indicate that neither the magnetospheric accretion model nor the disk wind model are individually capable of reproducing the observational interferometric data obtained for RU Lup in 2021 with VLTI GRAVITY. While it is to some degree possible to approximate either the GRAVITY spectrum or the characteristic sizes of the emission region with specific components and parameter settings, a combination of wind and magnetosphere is ultimately required to reconcile model and observation. Both of the hybrid scenarios provide significantly better fits to the observational data and do not only reproduce the trends in principle, but come close to matching the normalised line flux and characteristic sizes.

5.1 Non-axisymmetric magnetospheric accretion models

The introduction of the dipole tilt has a significant effect on the geometry of the Brγ brightness distribution, although the effective impact on the synthetic observables is tied strongly to its azimuth dependency. There is a qualitative change in the derived sizes across the line as we move from negative to positive azimuth angles (Fig. 5, see also Appendix B for a larger range of azimuths). This behaviour is caused by the interplay between the system inclination and the position of the two hemispheric accretion flows relative to the observer’s line of sight. At −90° azimuth, the upper hemisphere accretion column is inclined away from the observer, while at +90° azimuth it is inclined towards the observer. The latter case seems to favour the centrally peaked size profile, while the former shows that the sizes remain more constant at different velocities. However, the reason why the characteristic sizes are slightly larger at negative azimuth is not easily deduced, as there is also the baseline configuration to consider. While for the axisymmetric models the orientation of the baselines has a negligible impact, the azimuth dependency can change how their configuration probes the two hemispheric columns.

On the other hand, due to the nature of the photocentre shift computation, the effect of the baseline orientation on the derived offset position remains negligible even in the non-axisymmetric case. The differently oriented photocentre profiles at different azimuth angles can in this case then be directly attributed to the variation in the brightness distribution. There is once again a principal distinction between positive and negative azimuth angles, as the former lead to photocentre shift distributions that are more preferentially aligned with the disk axis, whereas the latter retain essentially the same angle of about 30° relative to the disk axis that was observed in the axisymmetric model.

While there is no obvious mechanism at play here that would directly explain why these models are associated with these specific alignments, the detection of different photocentre alignments at different azimuths is in itself noteworthy. We observe a similar change in the orientation of the photocentre profiles between the 2018 and 2021 datasets. In addition, the observational size profile appears relatively more flat in 2018, whereas the 2021 data shows a stronger dip at low velocities (Fig. 1). The difference in centre channel size between both epochs is on the order of 0.02 au at comparable uv planes, which is significant compared relative to the uncertainty of 0.001 au on those points. It is also similar to the level of size variation we derive at different azimuth angles. We do note that this azimuth-dependent change in half flux radius primarily affects in the wings of the Brγ feature rather than the centre.

Still, these azimuth-dependent effects may indicate that the two GRAVITY observations probe the RU Lup system at different points during its rotation. So while the switch to the non-axisymmetric dipole configuration does not resolve the fundamental limitations of the axisymmetric scenario with regards to fitting the RU Lup data, the inclusion of a tilted dipole offers one way to address the apparent time dependency of the observational results.

Fig. 7 Synthetic observables of the cool, compact hybrid model HM I (left) and the more extended, hotter hybrid model HM II (right). From top to bottom: centre channel image, line-to-continuum flux ratio, characteristic size, and photocentre shift per velocity channel of the Brγ emission region. The ellipse depicted in the image signifies the half flux radius of the geometric Gaussian disk model used to derive the characteristic size. In the depiction of the photocenter shifts, the stellar rotational axis, not the north axis, is aligned with the vertical coordinate axis. The observational photocentres were rotated by −9° with respect to Fig. 1 to put them in the same frame of reference.

Fig. 8 Decomposition of the synthetic observables into their constituent model component profiles for the first (left, HM I) and second (right, HM II) hybrid model. In this depiction the stellar rotational axis, not the north axis, is aligned with the vertical coordinate axis. The observational photocentres were rotated by −9° with respect to Fig. 1 to put them in the same frame of reference.

5.2 Disk wind models

The shortcomings of the pure wind models DW I and DW II are obvious, as they either fail to reproduce the behaviour of the characteristic sizes in the wings of the line or they produce a type of double peaked line profile inconsistent with the RU Lup data. The single peaked profile of DW 1 essentially replicates the fundamental characteristics of the magnetospheric accretion model. At the same time, it provides a worse fit to the data in terms of line width and especially photocentre shifts, thus offering no advantage over the magnetospheric accretion model. The double peaked configuration is more relevant, as it is the only type of low inclination model capable of matching trend of increasing sizes at high velocities. Indeed, DW II is identical in terms of model parameters to the disk wind component of the hybrid model HM II, where we successfully emulate the observational sizes in the wings by combining it with a hot magnetospheric central region. The line-to-continuum ratio depicted in Fig. 6 appears much larger than in the decomposition in Fig. 8 due to the lack of the magnetospheric continuum contribution, as will be discussed in Sec. 5.3.

We remind that we consider the half flux radii obtained with the Gaussian disk model to be characteristic sizes first and foremost, and do not preoccupy ourselves with their exact relationship to the true physical size of the emission region for the purpose of this work. Still, it is worth pointing out for the sake of prudence that the disk wind model images essentially depict a ring-like region, which we fit with a disk-like geometric model. This may create a greater discrepancy between physical size of the region and characteristic size, especially in the wings of the line. Figure 9 shows the HWHM of the Gaussian disk as an ellipse around the photocentre, compared to the underlying emission region. Especially at higher velocities, the relationship between both becomes more abstract, which may be a result of the mismatch in model morphologies. While this is of no concern to the comparison with the GRAVTIY data due to the consistent approaches between model and observational data treatment, a geometric ring model may better serve to derive values closer to the spatial extent of the true wind region.

5.3 Physicality of the hybrid model parameters

It is immediately obvious from Fig. 7 that the hybrid approach is superior to any of the other models detailed in this work in terms of matching the observational data. It can both reproduce the general trends that we see in the observational line and characteristic size profiles and can, in the case of HM II, even fit the observational data reasonably well, given the limitations of the semi-manual fitting approach described in Sec. 2. However, as the models are parametric in nature, there is a danger of enforcing a physically implausible parameter configuration in order to achieve a certain result. In this section we discuss the selected parameters in the context of the available literature on RU Lup.

While an examination of the physicality is difficult for parameters without explicit connection to quantifiable observables, there are those for which observational evidence or implicit constraints are more readily available. Chief among them are the stellar parameters, which are only partially shared between the hybrid models. For HM I and II the stellar mass was set to 1.2 M_* and 0.8 M_*, respectively, the stellar radii to 2.4 R_* and 2.5 R_*, and the stellar temperature to 4050 K. These masses are well within the broad range of possible values proposed by Alcalá et al. (2017). They estimate stellar mass and mass accretion rates from four different evolutionary models for pre-main sequence stars, thus yielding four different mass estimates between 0.43 M_* and 1.21 M_*. They also report a stellar radius of 2.39 ± 0.55 R_* and a photospheric temperature of 4060 ± 187 K.

The mass accretion rate for HM I was set to 10 × 10⁻⁸ M_⊙ yr⁻¹ and for HM II to 23 × 10⁻⁸ M_⊙ yr⁻¹. This agrees with our own estimates of the accretion rate for RU Lup based on the GRAVITY Brγ data. In GRAVITY Collaboration (2023b) we use the empirical accretion luminosity to line luminosity relationship given by Alcalá et al. (2014) to determine $M_{\mathrm{a}\mathrm{c}\mathrm{c}}={13.38}_{-6.66}^{+15.79}\times {10}^{-8}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$ for the 2021 data and ${17.31}_{-8.54}^{+20.78}\times {10}^{-8}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$ for the 2018 data of RU Lup. The highest estimate available in the literature is found by Siwak et al. (2016) at 30 ×10⁻⁸ M_⊙ yr⁻¹.

Since the magnetospheric accretion model does not treat the gas temperature in the funnel flows self-consistently, one needs to be cautious to choose a maximum accretion flow temperature consistent with observations. Hard limits to the chosen combination of accretion rate and temperature can be easily enough determined by looking at extreme cases. At either very low mass accretion rates or temperatures (≤6000 K, ≤ 10⁻⁹ M_⊙ yr⁻¹) barely any Brγ emission is produced and the image is almost completely continuum dominated. Increasing temperature and accretion rate simultaneously will lead to veiling effects as the gas starts emitting in the continuum and becomes increasingly optically thick. This behaviour was already observed in the magnetospheric accretion study of Muzerolle et al. (2001). They found that, for all the investigated hydrogen transitions, a critical threshold was reached at temperatures ≥10000 K and accretion rates of $ge; 10⁻⁶ M_⊙ yr⁻¹, at which point the emission line is turned into an absorption feature. At lower densities however, the gas never becomes sufficiently optically thick for this effect to occur. In their work they derive optimal ranges of mass accretion rates and temperatures for T Tauri stars from the continuum luminosity and the Paβ/Brγ line ratios. For our hybrid models, we paired the mass accretion rate of HM I (log(M_acc) = −7) with a maximum magnetospheric temperature of 7100 K and the accretion rate of HM II (log(M_acc) = −6.6) with a temperature of 8600 K. According to Fig. 16 in Muzerolle et al. (2001), the excluded temperature ranges for log(M_acc = −7) are given as below 6500 K and above 10 000 K and for log(M_acc = −6.6) below 6000 K and above 9500 K. Our chosen accretion and temperature parameters for both hybrid models lie well within this range of possible, or not excluded, combinations.

Fang et al. (2018) show for a population of T Tauri YSOs, RU Lup among them, that ${\dot{M}}_{\text{wind }}/{\dot{M}}_{\mathrm{a}\mathrm{c}\mathrm{c}}$ ratios can be on the order of slightly above 1, depending on the gas temperature, with a median ratio of 0.1 derived from a number of atomic emission lines for their sample. This is similar to the findings of Watson et al. (2016) and is consistent with the global mass loss rate we set for the hybrid models, where we use ${\dot{M}}_{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}={\dot{M}}_{\mathrm{a}\mathrm{c}\mathrm{c}}$ for HM I and ${\dot{M}}_{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}\approx 0.5{\dot{M}}_{\mathrm{a}\mathrm{c}\mathrm{c}}$ for HM II. However, these values are derived for the low velocity component of the outflows, with velocities on the order of a few tens of km s⁻¹. For the high velocity component, which is more in line with the velocities we require for our wind components, Fang et al. (2018) presents maximum ratios on the order of 0.1 or less. This would suggest that our wind component is significantly stronger than what has been observed in these studies. We also note that the high ${\dot{M}}_{\text{wind }}/{\dot{M}}_{\mathrm{a}\mathrm{c}\mathrm{c}}$ ratios found for the low velocity component are based on more extended wind regions than the ones used in our models.

It may be useful to discuss the dark disk radius, which in both hybrid models is set to 13 R_*, in the context of gas and dust in the innermost parts of the disk. Historically, the simplified view of inner disk opacities assumes that the close environment of a young star is divided into an inner optically thin purely gaseous disk and an outer optically thick dusty region, separated by the sublimation radius (Natta et al. 2001). Although some studies propose that the opacities of the gaseous cavity inwards of the sublimation radius behave in a more complex manner than this image suggests (see e.g. Muzerolle et al. 2004), the gaseous inner disk is only likely to become optically thick if the accretion rate is sufficiently large. If we, in first order, assume that in our case the midplane within the inner cavity is close to transparent to Brγ radiation, then the change in opacity would be predominantly connected to the presence of dust. In GRAVITY Collaboration (2021b), the K band continuum emission is traced to a region with a half flux radius of 0.21 au using a geometric Gaussian ring model. The general assumption is that this region is a proxy for the hot inner dust wall, which is thought to dominate the disk continuum emission at NIR K-band wavelengths. In the same publication, we compare the K-band size to a sublimation radius derived from a simple model of dust absorption and emission (Monnier & Millan-Gabet 2002) and this way find a smaller value of 0.05–0.1 au. At a dark disk radius of 13 R_* (0.15 au), the disk in the hybrid models becomes transparent just in between these two regions.

We set the rotation period in HM I to 9.05 days and in HM II to 7 days, which is consistent with known general distributions of rotational periods in low mass YSOs. These can range up to about 10 days, as was shown in a number of studies based on the data from the Kepler space telescope (e.g. Rebull et al. 2020). However, for the specific case of RU Lup, the rotational period has been estimated to be on the lower end of this dispersion. Stempels et al. (2007) investigated the periodicity of radial velocity variations for this object and derived a stellar period of 3.71 days, consistent with earlier photometric periodicity studies, for example by Hoffmeister (1965) and Drissen et al. (1989). Frasca et al. (2017) use the VLT X-Shooter spectrograph to determine the v sin(i) of RU Lup to be 8.5 km s⁻¹, which translates to a stellar period of 5.1 days when using the inclination fitted from the GRAVITY NIR continuum data (${20}_{-8}^{+6}$ degrees, see GRAVITY Collaboration 2023b). Even using the highest deviation of 26 degrees within the stated 1σ range will only yield a rotational period of 6.5 days in this case. It is not clear why our best fit models point to stellar periods that are significantly slower than known observational results, but the question of the model stellar period is ultimately also tied to the size of the magnetosphere and thus to the accretion regime of RU Lup.

Fig. 9 Images of the disk wind model DW I at different velocities across the Brγ line and at 20° (top row) and 45°(bottom row) inclination. The green ellipses depict the HWHM of the Gaussian disk model used to derive the characteristic sizes of the emission region. The Gaussian model is fixed at 20° inclination in both cases to remain consistent with the treatment of the original GRAVITY data. It is centred on the photocentre of each image, as recovered from the differential phases. The blue points signify the HWHMs obtained for the individual six baselines of the GRAVITY UT configuration that was used to observe RU Lup in 2021.

5.4 The accretion regime of RU Lup

The fundamental distinction between the two hybrid models HM I and HM II lies with how they treat the co-rotation criterion. HM I represents a parameter combination that fully respects the criterion: the sizes, stellar parameters, and rotational period were selected so that R_mo ≈ R_co. The magnetosphere is then truncated on the order of the co-rotation radius R_co, which is consistent with the idea of stable magnetospheric accretion. In numerical MHD simulations, a further distinction arises between the regimes of stable (0.7 R_co ≤ R_mo ≤ R_co) and unstable (R_mo < 0.7 R_co) accretion (Romanova & Owocki 2015). If the magneto-spheric truncation is more compact, then the disk can penetrate into the magnetosphere and the gas is accreted from close to the star along a set of spontaneously forming accretion tongues. Additionally, for R_mo > R_co, angular momentum considerations dictate that the accreting matter would be at least partially ejected at the disk-magnetosphere interface and thus lead to less steady, more time-variable accretion (Blinova et al. 2016). This is referred to as the ‘propeller regime’. In Romanova & Owocki (2015), a ‘fastness parameter’ $w_s={\left(\frac{R_{\mathrm{m}\mathrm{i}}}{R_{\mathrm{c}\mathrm{o}}}\right)}^{(3/2)}$ is used to distinguish between different accretion regimes and also between strong and weak propellers. They define a system in a weak propeller regime as one where w_s ≈ 1, while a w_s at multiples of 1 indicates a strong propeller with ejection exceeding the remaining accretion.

For HM II, we decided to relax the co-rotation criterion and set the truncation radius and the stellar period independently to obtain the best possible fit to the GRAVITY data. The chosen truncation radius of 7+2 R_* compares to a co-rotation radius of 5.72 R_* at the selected P_rot = 7 days. Our configuration of HM II translates into a fastness parameter of 1.35, which would indicate that – if we straightforwardly apply these concepts to our analytical model – it still exists in the regime of a weak propeller. HM II could then still be consistent with a system that is principally in accretion, even if some material is ejected at the magnetospheric boundary. Such a scenario of non-steady accretion in a weak propeller regime was discussed in Zanni & Ferreira (2013), see also Pantolmos et al. (2020).

However, this would contradict the idea that RU Lup is in a regime of unstable accretion, as proposed by Siwak et al. (2016). Their examination of RU Lup photometric data from MOST observations taken in 2012 and 2013 show highly irregular light curves for this object. Based on this lack of a clear periodicity in the brightness variability, they argue that the chaotic behaviour is the result of a randomised distribution of temporary accretion hot spots on the stellar surface caused by the unstable accretion tongues. This is in agreement with the theoretical predictions of unstable accretion light curves from numerical simulations (Romanova et al. 2008), and indeed this regime of unstable accretion is associated predominantly with higher accretion rates such as found for RU Lup. It is not clear whether the remaining accretion in a system in the propeller regime might produce a similar effect.

Another question as to whether the analytical stable accretion model yields plausible physical characteristics of the magnetic field. While we do not set a magnetic field strength explicitly, the chosen stellar parameters, mass accretion rate, and truncation radius can be implicitly related to a magnetic dipole field strength. Using the formulation in Hartmann et al. (2016), we can determine the field strength in kG via $$B_{\mathrm{k}\mathrm{G}}={\left(\frac{R_{\mathrm{t}\mathrm{r}}{M}_{0.5}^{(1/7)}{\dot{M}}_{-8}^{(2/7)}}{18\gamma R_2^{(12/7)}}\right)}^{(7/4)},$$(5)

where M_0.5 is the stellar mass in units of 0.5 $M_{\odot },{\dot{M}}_{-8}$ the mass accretion rate in units of 10⁻⁸ M_⊙ yr⁻¹, R₂ the stellar radius in units of 2 R_⊙, and γ is an uncertain factor ≤1. Assuming the most favourable case with γ = 1, and using the inner magnetospheric radius as truncation radius (i.e. 6 R_* for HM I and 7 R_* for HM II), the field strength for HM I comes out to 1.5 kG, whereas for HM II we compute a 2.6 kG dipole field. Although there are no published measurements of the RU Lup field in the literature, a default assumption for the magnetic fields of T Tauris stars is a range of about 1–2 kG, so the implicit field strength of the HM II configuration would be atypically large. We do note, however, that this formula is approximate and based on simplifications of the mechanics involved. If indeed the accretion process in RU Lup is complex beyond the case of simple stable magnetospheric accretion, as the previous considerations imply, then these magnetic field values should not be adopted without question.

The question of the RU Lup accretion regime is ultimately not easily answered. The stable analytical accretion model is also only an approximation of a more complex truth, and the fact that some of the parameters we derive from it are not immediately consistent with other observations is not unexpected. Another example for a potential additional effect that is not considered at all outside certain numerical simulations is the idea of a failed disk wind, meaning a wind that is launched close the disk-magnetosphere interface, but which fails to escape the gravitational field of the star. Instead, the wind material falls back onto the magnetosphere and is accreted onto the star, effectively enlarging the magnetospheric emission region and masking a more compact actual magnetosphere in the process (Takasao et al. 2022).

A different potential explanation for the size of the best fitting magnetosphere may lie with the treatment of the observational data rather than with the models themselves. Contrary to the treatment of the continuum data in GRAVITY Collaboration (2021b), the visibility modelling of the Brγ line emission data did not consider the effects of a large scale halo component that might account for a significant part of the line flux. Such a halo component, which may for example reflect the influence of Brγ photons scattered on the disk surface, was omitted from the treatment as the primary intent behind GRAVITY Collaboration (2023b) was to compare general trends among a number of objects for a simple geometric approach. Including a halo would effectively reduce the derived characteristic size for the resolved emission region compared to the results presented in Fig. 1. Whether this effect could be sufficiently strong to reduce the observational characteristic sizes to a point where we may be able to fit them similarly well as we did with HM II, without exceeding the co-rotation radius, is at this point entirely speculative and would depend on a number of other factors. Still, it may be a point worth considering for future investigations into similar questions.

5.5 The photocentre shift problem and the possible role of other emission components

One of the primary issues of the model analysis in GRAVITY Collaboration (2023b) is that the synthetic photocentre positions obtained from the images are offset by a factor four to five when compared to the GRAVITY photocentre shifts. Neither the consideration of a dipole tilt, nor the introduction of an additional wind component, do much to resolve this issue. On the contrary, the photocentre profile associated with the disk wind shows even substantially larger photocentre shifts than the magnetosphere. As a consequence, the hybrid photocentre profiles become even more elongated, as shown via the decomposition in Fig. 8.

The divergence between the observational and synthetic profiles is potentially rooted in some aspect of the observed system that is not reflected in the models, such as the presence of an additional Brγ emission component. Such a component could, for example, come in the form of a hot MHD stellar wind, launched along the open magnetic field lines from the polar regions of the star (Ferreira et al. 2006). However, it is not clear whether such a stellar wind would even address the issue of the photocentres in principle, let alone be consistent overall with the GRAVITY data. In a low inclination system such as RU Lup, a collimated outflow, launched from the stellar pole, almost directly along the line of sight would be expected to cause significant blueshifted absorption, whereas GRAVITY on the contrary picked up blueshifted excess emission in the Brγ spectrum. Given that the redshifted stellar wind component is also likely to be at least partially blocked by the star itself at this inclination, there should be a substantial difference in both characteristic size and photocentre shift profiles between the red and the blue arm. It is equally not clear whether such a collimated wind could simultaneously reproduce the increase in size in the wings and the more compact distribution of photocentres.

Any component that may help to resolve the question of the overestimated photocentre shifts fundamentally needs to lead to a decrease in offset magnitudes without an accompanying reduction in characteristic size. This would be achieved by a system that remains relatively centrosymmetric across the line at scales even up to the truncation radius, but then features asymmetries in the extreme close environment of the star. A disk-like structure that would extend to the innermost regions of the system, viewed at low inclination, may provide the necessary degree of centrosymmetry across the line, but lacks the asymmetric elements. In this context we could consider again the possibility that RU Lup may be in the regime of unstable accretion. If the gaseous matter penetrating the magnetosphere behaved essentially as such a disk, then the formation of accretion tongues very close to the star could introduce the necessary asymmetry that would lead to offsets in the photocentre positions. Answering this question would require some dedicated modelling, as without knowledge of the relative Brγ flux between the tongues and the remaining emission region, it is difficult to comment on whether this would be a viable solution. This is also true for the introduction of a halo component, which is another possible way to address this question. The scattered light halo may in principle have a similar compressing effect on the photocentre shifts, assuming that it itself is largely centrosymmetric. Whether this effect could be sufficiently large to explain the degree of divergence we detect between model and observation is again a question of the flux ratios between halo and resolved region. Still, as of now, these ideas may hold the most promise to resolve the issue of the photocentre shifts.

6 Summary

We computed the spectro-interferometric characteristics of six models of the Brγ line emission region in the innermost environment of low mass young stars. These models featured different configurations of emission components in the form of disk wind outflows and magnetospheric accretion. From the synthetic model data, which included the normalised line-to-continuum flux ratio, the interferometric visibilities, and the differential phases, we determined characteristic sizes and the position of the emission region photocentres for ten spectral channels across the Brγ line. We then compared these synthetic quantities to observational results derived from VLTI GRAVITY data for the classical T Tauri star RU Lup.

Our goal was to explain a number of features in the observational data that were not well reproduced by previous modelling efforts presented in GRAVITY Collaboration (2023b). Specifically, we tried to determine whether the addition of a larger scale outflow component onto an accreting magnetosphere would deliver results consistent with an asymmetric line profile with blueshifted excess emission, an increase in characteristic radii towards the high velocity wings of the line profile, and a compact distribution of photocentres across the line. To summarise the most important results of this paper:

• Both axisymmetric and non-axisymmetric magnetospheric accretion models produce either centrally peaked or flat characteristic size profiles across the line. No parameter configuration leads to a size increase at high velocities and a minimum size at the line centre, as observed for RU Lup. While it is possible to fit the central component of the Brγ line with a hot magnetosphere, the blueshifted excess emission that we detect observationally cannot be emulated.

• The variability in equivalent width, peak flux ratio, and photocentre shifts that we find between the 2018 and 2021 epochs of the observation is on the order of the variation that we detect when comparing the line profile of the nonaxisymmetric magnetospheric accretion model at different azimuth angles. This could indicate that between 2018 and 2021 we are in observing a tilted dipole at different phases during the stellar rotational period.

• Only disk wind models with sufficiently fast winds to produce double peaked spectral features at the spectral resolution of GRAVITY are associated with an increase in emission region sizes at high velocities. We do not manage to reproduce both the centrally peaked line profile and the increase in the size profile simultaneously with a disk wind alone.

• Hybrid models, combining a strong magnetospheric component that dominates the central channels of the line and a fast double peaked wind component, can reproduce the trends in both line profile and sizes at the same time as parameters that appear mostly in agreement with previous observational findings on the physical characteristics of RU Lup. The magnitude of the model photocentre shifts overestimates the observational photocentre profile across all model variations by factor of four or more, including the hybrid models.

• Our best fitting hybrid model manages to well reproduce the 2021 Brγ spectral line and characteristic sizes across the line. Still, in that case, the truncation radius exceeds the co-rotation radius by more than 50%, meaning that part of the gas is found beyond the centrifugal barrier, where it is expected to be expelled and not accreted.

In conclusion, we find that the success of the hybrid approach to modelling the Brγ region strongly supports the idea of a multi-component emission environment for certain class II low mass YSOs with large mass accretion rates such as RU Lup. Discrepancies between the expected magnetospheric accretion characteristics and our best fitting parameters likely indicate that the accretion-ejection process in RU Lup is more complex than what is captured by our analytical models. Another explanation may lie in the potential effects of a possible unresolved Brγ halo component, which was not included in the treatment of the GRAVITY data in GRAVITY Collaboration (2023b). Neglecting the halo when deriving the observational results would potentially lead to an overestimation of the size and underestimation of the photocentre shift magnitudes, depending on the flux ratio between resolved and unresolved region. If that is indeed the case, it may be possible to attain a similarly good fit to the remodelled data with a smaller magnetosphere that respects the co-rotation criterion.

Acknowledgements

This work was supported by the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU co-funded by CEA and CNES. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 742095; SPIDI: StarPlanets-Inner Disk- Interactions, http://www.spidi-eu.org, and grant agreement no. 740651 NewWorlds). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the Jean-Marie Mariotti Center Aspro service¹. V.G. was supported for this research through a stipend from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne, and from the Bonn-Cologne Graduate School of Physics and Astronomy (BCGS). A.C.G. has been supported by PRIN-INAF-MAIN-STREAM 2017 “Protoplanetary disks seen through the eyes of new-generation instruments” and by PRIN-INAF 2019 “Spectroscopically tracing the disk dispersal evolution (STRADE)” A.A. and P.G. acknowledge supported by Fundação para a Ciência e a Tecnolo- gia, with grants reference UIDB/00099/2020, SFRH/BSAB/142940/2018 and PTDC/FIS-AST/7002/2020. R.G.L. acknowledges support by Science Foundation Ireland under Grant No. 18/SIRG/5597.

Appendix A Parameter dependencies of the models

While the parameter descriptions laid out in sections 3.2 and 3.3, and their associated equations given in Hartmann et al. (1982) and Knigge et al. (1995), may outline the general impact of a parameter change on the densities and velocities in the system, they cannot straightforwardly connect the effect of such a variation with the resulting emission line profiles and possibly even less so with the interferometric quantities. Due to the complex interplay between atomic populations, radiative transitions, and the spatial distribution of velocities and densities, it is more instructive to directly analyse the effective change in line strength, emission region size, and photocentre distribution when increasing or decreasing a single parameter by a certain margin. In this manner, both the total sensitivity of the model as well as the partial sensitivities of these individual quantities to changes in a specific parameter may be constrained and compared. The limitations to this approach lie predominantly in the need to compute a new model for each parameter change and the time investment required to do so. Given the relatively large amount of parameters that we consider, in principle, fit for variation up to a point, we limit this analysis to one increase and one decrease relative to the parameter value of an arbitrarily chosen reference model. The step size of the variation is not set consistently but rather chosen based on relevant upper and lower limits from the literature for parameters which have known observational constraints, such as the stellar radius or mass. Other parameters, for which such constraints are not available, are varied to a degree that results in a similar change in the line-to-continuum flux ratio. While a purely theoretical study of the models would certainly look at a broader range of parameter changes, we remind that the purpose of this analysis was primarily to determine which quantities were the most relevant to the interpretation of specific observational results.

We differentiate between individual parameters based on their influence on the results. To this end, we present multiple ways to quantify the impact. First, we consider the relative change ΔO as a percentage of the respective reference model parameter: $$\mathrm{\Delta }O=\left|\frac{O_v-O_r}{O_r}\right|.$$(A.1)

Here, O is one of the observables we compute from the model images, i.e. line-to-continuum flux ratio, characteristic size of the emission region or photocentre offset of the emission region. The subscript v designates that the observable belongs to the variation model, while r designates the reference model. We choose to employ a relative measure of the change primarily due to the very different scales between these three observables. If we were to consider the absolute change, the large flux ratio values would dominate the average over the small sizes and photocentre shifts. In particular, we consider here the mean relative change between the parameter increase and the parameter decrease. We further compute the unweighted average of this metric over all considered wavelengths across the Brγ feature. Then, we subsequently set this value in relation to the relative parameter step ΔP, which is computed analogously according to Eq. A.1 for a variation model parameter P_υ and a reference model parameter P_r. Thus, we arrive at a measure for the relative mean observable change, adjusted by the parameter step size. We refer to this quantity as the model observable O response ξ_O to a parameter P change ΔP.: $${\xi }_O=\frac{\mathrm{\Delta }O}{\mathrm{\Delta }P}.$$(A.2)

A.1 Magnetospheric accretion

In Table A.1 we present the relative effects of a parameter change on the magnetospheric accretion model observables, as well as the response metrics computed based on Eq. A.2. When ranked by their respective ξ_tot values, the maximum temperature in the magnetospheric funnel flows T_mag, the stellar Radius R_*, and the inner magnetospheric truncation radius R_mi come out top in terms of their overall impact on the resulting model observables. Figure A.1 shows the deviations of the line profile, the size, and photocentre shifts from the reference model for the across all channels of the Brγ feature for these quantities. Here, the effect of a magnetospheric temperature variation on all three observables stands out as particularly noteworthy when taking into account the relatively small variation step width. The relative changes in photocentre shifts and sizes are on the order of factor 2 to 4 higher when compared even to the other two high impact parameters and on the order of factor 10 when compared to the mass accretion rate of the system, for which the variation step width of 22% is similar to T_mag at 16%. Figure A.1 shows that an increase in magnetospheric temperature from the reference value of 8600 K leads to a decrease in continuum-normalised flux and simultaneously to an increase in the emission regions characteristic size. The former suggests that, at such high temperatures, a significant amount of continuum emission is produced in the accretion columns due to bound-free transitions as the hydrogen becomes increasingly ionised. The larger characteristic size stems from a flatter distribution of flux across the entire magnetosphere, which leads to a spatially more extended full width at half maximum of the Gaussian disk model used to fit the characteristic size. For the same reason, the photocentre appears slightly less extended, as the more evenly distributed emission throughout the very hot magnetosphere symmetrises the emission region relative to the reference model.

Similar reasoning can be applied when breaking down the effects of other parameter changes, although in the case of the stellar radius one should take care to separate the physical influence of the radius from model scaling effects. Since we define the size of our truncation radius and magnetospheric widths in units of stellar radii, we are effectively changing the absolute physical size of the magnetosphere by varying the stellar radius and thus also the position of the field line anchor points as they appear in those equations. This in turn explains why the impact a stellar radius variation seems to beat an adjustment of the truncation radius, both in terms of ξ_tot and also in the partial metric for the characteristic size alone. The characteristic size of the system is a product of both R_mi and, albeit to a significantly lesser degree, δR_m, which are both defined in units of R_*, whereas setting the truncation radius directly will only affect R_mi.

By contrast, it appears that the stellar mass M_* and photo-spheric temperature T_* only marginally influence the outcome of the computations. This is particularly true for the photospheric temperature, which, even when accounting for the very small parameter variation of about 2%, produces a negligible change in observables compared to the other quantities. The low impact of T_* on the interferometric quantities is not surprising, given that the characteristic size of the emission region is based on the continuum extracted pure line visibilities and the pure line differential phase. The fact that some variation in photocentre shift with changing photospheric temperature can be observed at all is likely caused by the shifting weights between the stellar continuum and the magnetospheric continuum. Indeed, the former will produce a small non-zero continuum phase at higher inclinations due to the presence of the ring-like shock region close to the pole, while the magnetospheric continuum is centrosymmetric for an axisymmetric configuration, such as used here.

Finally, the stellar period does not affect the flux ratio or characteristic size to the same degree as R_*, but has a comparable effect on the photocentre shift, due to increased line broadening at higher rotational velocities of the star. Since the magnetospheric funnels rotate at the same velocity, observing the rotating gas stream at a non-zero inclination will lead to a broader line profile. The change in flux ratio caused by this is not necessarily large, as demonstrated by the low ξ_Flux of only 7%. The redistribution of flux to higher velocities, however, translates to a relatively large change in magnetospheric symmetries at those velocities and thus the more drastic variations in the photocentres per velocity component.

A.2 Disk wind

In Table A.2 we present an analogous ranking of the response metric ξ_tot for the disk wind model. It demonstrates that again three parameters perform significantly stronger compared to the rest in terms of their impact on the three observables. For the wind model, we find that the critical height z_crit, the velocity law parameter β_wind, and the global mass loss rate ${\dot{M}}_{loss}$ have far above average impact on the resulting line profiles, characteristic sizes, and photocentre shifts. The response metrics for these three parameters are broadly similar across the three observables. For each of them, the step width-weighted parameter variation causes a particularly large relative change in the photocentre profile, for which the ξ_PC is about a factor two larger than the corresponding flux ratio metric, whereas the respective ξ_Size values are similar to, but slightly smaller than, the flux ratio metric. In Figure A.2 we show the corresponding changes in flux ratio, size, and photocentre offset.

It is obvious that a higher cutoff z_crit of the wind leads to a weaker line-to-continuum ratio. This is expected, as the density along a stream line falls with increasing distance from the midplane. As a consequence, by choosing a higher cutoff, we effectively remove the denser parts of the wind from the model, thus lowering the Brγ emission contribution of the wind. The size, on the contrary, increases slightly when the cutoff is higher, since the effective inner diameter of the brightness distribution is enlarged due to the collimation angle of the field lines, see Fig. 3. The maximum shift magnitude among the photocentre offsets is even more extended for the same reason.

By contrast, for both mass loss rate and the velocity law exponent, parameter increases lead to stronger flux ratios, but smaller emission region characteristic sizes and more compact photocentre shift profiles. For the wind velocity exponent β_wind, this behaviour is straightforward to understand. A higher value for beta;_wind decreases the acceleration of the wind. In turn, in requires the wind to travel further to reach the same velocity, which leads to the size increase due to the angle of the stream lines. The velocity also influences the density inversely, as the density of gas particles in a slow moving outflow is higher, and so is subsequently then the gas emission. The change in characteristic size and photocentre offsets following such a parameter variation is most pronounced at high velocities. As there is now more or less, for decreasing or increasing β_wind respectively, gas moving at those high velocities, the region appears larger or smaller at those wavelengths. The parallel decrease in flux ratio that comes with a larger region is largely irrelevant to the characteristic size, as long as the hydrogen is dense enough to produce Brγ emission, since the characteristic size is obtained from the continuum-corrected visibilities.

An analogous argument can be made for the changes induced by a variation of the global mass loss rate. The local mass loss rate $\dot{m}(\mathrm{R})$ at a certain distance R from the central star is $\propto {\dot{M}}_{loss}$, so a high global mass loss rate leads to an increase in local mass loss at every distance and thus increased wind density at every distance. However, in this case, the α_wind parameter, as well as the chosen collimation angle of the wind stream lines, lead to a stronger increase in density closer to the inner edge of the disk. Under such circumstances, a fit of a Gaussian disk model would yield a smaller characteristic size.

The stellar parameters appear to again have some of the weakest effects on the three observables. The rotation period in particular does not affect the model at all, since the angular velocity component of the wind is fully determined by Keplerian velocity at the stream line anchor point. As a notable difference compared to the magnetospheric model, the stellar radius seems to rank below both inner radius and width of the wind region, even though both are defined in units of R_* and thus naively the same reasoning as laid out there should apply. The equations governing the density and velocity profiles of the wind model, however, show that dependencies on length scales and sizes mostly come in the form of a ratios between different different projected radial distances, which are all defined in stellar radius units. It follows then that these quantities are relatively less affected by a change in R_*. It is also the case here that the most impactful distance scale, the cutoff height z_crit, is defined in absolute astronomical units and does not vary with the stellar radius.

Finally, we note that the dark disk radius parameter R_dd, which we introduced to regulate the amount of detected flux from the back side of the disk, is of below average influence according to its computed ξ_tot. While this is true when considering the overall profile changes induced by a variation in R_dd, it does not properly illustrate the effect on profile symmetry. The dark disk radius R_dd is the only parameter that allows us to directly adjust the relative strength between blue and red wings of the observable profiles. This unique quality of the parameter is not well reflected in the response metric, given that we consider average metrics over the entire relevant wavelength range.

Fig. A.1 Observable dependencies on parameter changes for the top three magnetospheric accretion parameters as ranked in Table A.1. The left column depicts the changes in continuum normalised line flux. The centre column shows a characteristic interferometric size, obtained as the half width at half maximum (HWHM) of a geometric Gaussian disk model. The right column illustrates the spatial distribution of line emission photocentres at different velocity channels.

Fig. A.2 Observable dependencies on parameter changes for the top three disk wind parameters as ranked in Table A.2. The left column depicts the changes in continuum normalised line flux. The centre column shows a characteristic interferometric size, obtained as the half width at half maximum (HWHM) of a geometric Gaussian disk model. The right column illustrates the spatial distribution of line emission photocentres at different velocity channels.

Appendix B Azimuth dependency of the non-axisymmetric MA model

Fig. B.1 Line-to-continuum flux ratio, characteristic size, and photocentre shift of the MA II model at five different azimuth angles from −90° to +90°. The inclination of the image is fixed at 20° to conform to the observational value obtained from VLTI GRAVITY data for the inner disk of RU Lup in 2021. The images depict the emission region at 0 km/s.

Appendix C Grid profiles for HM II

Fig. C.1 2D model profiles for hybrid model II. Shown are, from top to bottom: Mass density, temperature, and the three spherical components of the velocity field.

References

All Tables

All Figures

	Fig. 1 Previous results for CTTS RU Lup, derived from VLTI GRAVITY observations in 2018 and 2021. Shown here are (from top to bottom) the line-to-continuum flux ratio of the Brγ spectral line, the characteristic size of the emission region across the velocity channels of the line, and the corresponding positions of the emission region photocentres.
In the text

Fig. 2 Example of a magnetospheric accretion model image, produced with the MCFOST radiative transfer code. Shown here is the continuum-subtracted Brγ emission region of a purely axisymmetric rotating magnetosphere with no additional emission components. Channel line maps such as this can be used to compute artificial interferometric observables which we are then able to compare to real observational data. The image here specifically shows the line map for the −69 km s−1 velocity channel.

In the text

Fig. 3 Schematic depiction of the magnetospheric accretion and disk wind models. Shown here are the geometrical parameters defining each component. These are: the rotational axis Ω, the magnetic dipole axis μ, the magnetic obliquity βma, the wind focal point (S) displacement d, the inner and outer radius Rmi and Rmo, and width δRmag, of the magnetosphere, the inner and outer radius Rwi and Rwo, and width δRwind, of the wind, the dark disk radius Rdd beyond which the disk becomes completely opaque, and the cutoff zcrit at which the isothermal wind instantly reaches its 104 K temperature. The lower hemisphere parts of wind and magnetospheric funnel are omitted in this depiction for the sake of clarity.

In the text

Fig. 4 Synthetic observables of the axisymmetric magnetospheric accretion model MA I (left) and the non-axisymmetric accretion model MA II (right) when viewed at an azimuth angle of −90°. From top to bottom: centre channel image, line-to-continuum flux ratio, characteristic size, and photocentre shift per velocity channel of the Brγ emission region. The ellipse depicted in the image signifies the half flux radius of the geometric Gaussian disk model used to derive the characteristic size. In the depiction of the photocenter shifts, the stellar rotational axis, not the north axis, is aligned with the vertical coordinate axis. The observational photocentres were rotated by −9° with respect to Fig. 1 to put them in the same frame of reference.

In the text

Fig. 5 Synthetic observables of the non-axisymmetric accretion model MA II when viewed at azimuth angles of −45° (left) and +45° (right). From top to bottom: centre channel image, line-to-continuum flux ratio, characteristic size, and photocentre shift per velocity channel of the Brγ emission region. The ellipse depicted in the image signifies the half flux radius of the geometric Gaussian disk model used to derive the characteristic size. In the depiction of the photocenter shifts, the stellar rotational axis, not the north axis, is aligned with the vertical coordinate axis. The observational photocentres were rotated by −9° with respect to Fig. 1 to put them in the same frame of reference.

In the text

Fig. 6 Synthetic observables of the slow disk wind model DW I (left) and the fast disk wind model DW II (right). From top to bottom: centre channel image, line-to-continuum flux ratio, characteristic size, and photocentre shift per velocity channel of the Brγ emission region. The ellipse depicted in the image signifies the half flux radius of the geometric Gaussian disk model used to derive the characteristic size. In the depiction of the photocenter shifts, the stellar rotational axis, not the north axis, is aligned with the vertical coordinate axis. The observational photocentres were rotated by −9° with respect to Fig. 1 to put them in the same frame of reference.

In the text

Fig. 7 Synthetic observables of the cool, compact hybrid model HM I (left) and the more extended, hotter hybrid model HM II (right). From top to bottom: centre channel image, line-to-continuum flux ratio, characteristic size, and photocentre shift per velocity channel of the Brγ emission region. The ellipse depicted in the image signifies the half flux radius of the geometric Gaussian disk model used to derive the characteristic size. In the depiction of the photocenter shifts, the stellar rotational axis, not the north axis, is aligned with the vertical coordinate axis. The observational photocentres were rotated by −9° with respect to Fig. 1 to put them in the same frame of reference.

In the text

	Fig. 8 Decomposition of the synthetic observables into their constituent model component profiles for the first (left, HM I) and second (right, HM II) hybrid model. In this depiction the stellar rotational axis, not the north axis, is aligned with the vertical coordinate axis. The observational photocentres were rotated by −9° with respect to Fig. 1 to put them in the same frame of reference.
In the text

Fig. 9 Images of the disk wind model DW I at different velocities across the Brγ line and at 20° (top row) and 45°(bottom row) inclination. The green ellipses depict the HWHM of the Gaussian disk model used to derive the characteristic sizes of the emission region. The Gaussian model is fixed at 20° inclination in both cases to remain consistent with the treatment of the original GRAVITY data. It is centred on the photocentre of each image, as recovered from the differential phases. The blue points signify the HWHMs obtained for the individual six baselines of the GRAVITY UT configuration that was used to observe RU Lup in 2021.

In the text

Fig. A.1 Observable dependencies on parameter changes for the top three magnetospheric accretion parameters as ranked in Table A.1. The left column depicts the changes in continuum normalised line flux. The centre column shows a characteristic interferometric size, obtained as the half width at half maximum (HWHM) of a geometric Gaussian disk model. The right column illustrates the spatial distribution of line emission photocentres at different velocity channels.

In the text

Fig. A.2 Observable dependencies on parameter changes for the top three disk wind parameters as ranked in Table A.2. The left column depicts the changes in continuum normalised line flux. The centre column shows a characteristic interferometric size, obtained as the half width at half maximum (HWHM) of a geometric Gaussian disk model. The right column illustrates the spatial distribution of line emission photocentres at different velocity channels.

In the text

	Fig. B.1 Line-to-continuum flux ratio, characteristic size, and photocentre shift of the MA II model at five different azimuth angles from −90° to +90°. The inclination of the image is fixed at 20° to conform to the observational value obtained from VLTI GRAVITY data for the inner disk of RU Lup in 2021. The images depict the emission region at 0 km/s.
In the text

	Fig. C.1 2D model profiles for hybrid model II. Shown are, from top to bottom: Mass density, temperature, and the three spherical components of the velocity field.
In the text