The GRAVITY young stellar object survey. III. The dusty disk of RY Lup

Context. Studies of the dust distribution, composition, and evolution of protoplanetary disks provide clues for understanding planet formation. However, little is known about the innermost regions of disks where telluric planets are expected to form. Aims. We aim constrain the geometry of the inner disk of the T Tauri star RY Lup by combining spectro-photometric data and interferometric observations in the near-infrared (NIR) collected at the Very Large Telescope Interferometer. We use PIONIER data from the ESO archive and GRAVITY data that were obtained in June 2017 with the four 8m telescopes. Methods. We use a parametric disk model and the 3D radiative transfer code MCFOST to reproduce the spectral energy distribution (SED) and match the interferometric observations. MCFOST produces synthetic SEDs and intensity maps at different wavelengths from which we compute the modeled interferometric visibilities and closure phases through Fourier transform. Results. To match the SED from the blue to the millimetric range, our model requires a stellar luminosity of 2.5 L (cid:12) , higher than any previously determined values. Such a high value is needed to accommodate the circumstellar extinction caused by the highly inclined disk, which has been neglected in previous studies. While using an effective temperature of 4800 K determined through high- resolution spectroscopy, we derive a stellar radius of 2.29 R (cid:12) . These revised fundamental parameters, when combined with the mass estimates available (in the range 1.3–1.5 M (cid:12) ), lead to an age of 0.5–2.0 Ma for RY Lup, in better agreement with the age of the Lupus association than previous determinations. Our disk model (that has a transition disk geometry) nicely reproduces the interferometric GRAVITY data and is in good agreement with the PIONIER ones. We derive an inner rim location at 0.12 au from the central star. This model corresponds to an inclination of the inner disk of 50 ◦ , which is in mild tension with previous determinations of a more inclined outer disk from SPHERE (70 ◦ in NIR) and ALMA (67 ± 5 ◦ ) images, but consistent with the inclination determination from the ALMA CO spectra (55 ± 5 ◦ ). Increasing the inclination of the inner disk to 70 ◦ leads to a higher line-of-sight extinction and therefore requires a higher stellar luminosity of 4.65 L (cid:12) to match the observed ﬂux levels. This luminosity would translate to a stellar radius of 3.13 R (cid:12) , leading to an age of 2–3 Ma, and a stellar mass of about 2 M (cid:12) , in disagreement with the observed dynamical mass estimate of 1.3–1.5 M (cid:12) . Critically, this high-inclination inner disk model also fails to reproduce the visibilities observed with GRAVITY. Conclusions. The inner dust disk, as traced by the GRAVITY data, is located at a radius in agreement with the dust sublimation radius. An ambiguity remains regarding the respective orientations of the inner and outer disk, coplanar and mildly misaligned, respectively. As our datasets are not contemporary and the star is strongly variable, a deeper investigation will require a dedicated multi-technique observing campaign.


Introduction
Because they are the sites of planet formation, the study of protoplanetary disks is one of the primary science drivers for several recent major observing facilities.ALMA, operating at (sub-)millimeter wavelengths (ALMA Partnership et al. 2015), and second-generation adaptive-optics imagers operating in the optical/near-infrared(NIR), reach comparably high angular resolutions, providing clear views of the disks at the 10 au scale for the typical star-forming regions (i.e., at ∼ 140 pc).One of the many outstanding findings to come from such observations is that most of the disks are non-symmetric in nature when angular GTO program with run ID 099.C-0667 resolution is improved: resolved disks show gaps, rings, spirals, vortices, and shadows (Andrews et al. 2018;Long et al. 2018;Benisty et al. 2015;de Boer et al. 2016;Pohl et al. 2017;Benisty et al. 2017;Avenhaus et al. 2018).The origin (and universality) of these features is actively debated, given the clear possibility that they are tracing the dynamical interaction between a forming planet and its parental disk.However, other mechanisms like snow lines (Zhang et al. 2015), non-ideal MHD effects (Béthune et al. 2016), zonal flows (Lorén-Aguilar & Bate 2015), and selfinduced dust-traps (Gonzalez et al. 2015) also have the capacity to produce rings, gaps, and spirals without the need for planets.
While scattered light imaging in the optical and the NIR with, for example, SPHERE (Beuzit et al. 2019) or imaging at (sub-) millimeter wavelengths with ALMA provides de-tailed views of the outer regions of the disks (i.e., from ∼10 to ∼500 au), these instruments are not directly sensitive to the very inner regions.Typically, the inner parts of the disks (< 5-10 au) escape from direct view because of the use of a coronagraph with SPHERE, for example.This is perhaps an unfortunate situation, because the techniques most currently used today to search for extrasolar planets, namely those involving radial velocities and transit surveys, are naturally biased towards planets with small periods and thus small separations.Nevertheless, the results from these techniques clearly indicate that planets are frequent around solar-like stars (Clanton and Gaudi, 2016;Baron et al., 2019;Fernandes et al., 2019).
The standard scenario for planet formation is called core accretion (Pollack et al. 1996;Rice & Armitage 2003).The main caveat with the original version of this scenario is that the timescale to form a Jupiter-like planet (at 5 au from the star) is longer than the typical gas-disk lifetime, with the situation becoming even worse for planets located farther away from the central star.Also, there is not enough mass available in-situ to form giant planets closer in, for example in Earth-like orbits or closer in.Nevertheless, several ways of mitigating the timescale and mass problems have been suggested.On one hand, planets can form further out and migrate inward (Alibert et al. 2005).On the other, the streaming instability in the disk midplane (Youdin & Goodman 2005;Johansen et al. 2014) and/or pebble accretion (Ormel & Klahr 2010;Lambrechts & Johansen 2012) can significantly speed up the planet formation process.Therefore, knowledge of the disk properties in the midplane is critical, and is even more so in the inner regions.
Complementary to SPHERE and ALMA, long-baseline NIR interferometers provide the necessary spatial resolutions of a few milliarcseconds (mas), which is typically a fraction of an astronomical unit at 140 pc.These instruments allow us to study the dust distribution, composition, and evolution of protoplanetary disks, which provides clues for understanding planet formation.Moreover, a sharper view of the dust distribution in the inner astronomical unit is also crucial to better understand the timescales of disk evolution.Previous NIR interferometric observations focused mostly on the brightest end of the pre-main sequence stellar distribution for lack of sensitivity.In practice therefore, Herbig Ae/Be stars constitute the bulk of the observed sample (e.g., Lazareff et al. 2017), with a few rare exceptions for the brightest T Tauri stars observed with the Keck Interferometer (Eisner et al. 2007(Eisner et al. , 2010)), AMBER (Vural et al. 2012), and PIONIER (Anthonioz et al. 2015).With a significant gain in sensitivity, GRAVITY (Gravity Collaboration et al. 2017) provides the opportunity to observe, at high angular resolution, a larger sample of these young, low-mass stars, which appear more compact and fainter in the NIR than the more massive Herbig stars, which is due to their cooler photospheres.
Here, we present the observations of the Solar-like pre-main sequence star RY Lup obtained with GRAVITY in the K-band.RY Lup is located at a distance of 158 pc (Gaia Collaboration et al. 2018).Although a range of spectral-type estimates are available in the literature, RY Lup is now classified as K2 based on X-Shooter spectroscopic data (Alcalá et al. 2017).These latter authors also derived a luminosity of 1.7 L and a stellar mass of 1.47M (from Siess et al. (2000) tracks), and revised the effective temperature to slightly lower values (T eff = 4900 ± 227 K) compared to previous studies (e.g., 5080 K derived from Frasca et al. 2017).Using the evolutionary tracks of Baraffe et al. (2015), Frasca et al. (2017) derived a mass of 1.4 M and an age of 10.2 Ma, while most of the young stellar objects of the Lupus association have an age of 1-3 Ma (see their Fig. 6).
The SED of RY Lup shows a significant IR excess, as well as a modest UV excess (Evans et al. 1982;Gahm et al. 1989).RY Lup also exhibits strong photometric variability.Manset et al. (2009), using simultaneous BV polarimetric and UBV photometric observations, showed that the polarization is high (3.0%) when the star is faint and red (V = 12.0 mag, B-V = 1.3 mag), and low (0.5%) when it is brighter and bluer (V = 11.0 mag, B-V = 1.1 mag).The photometric and polarimetric variabilities also share a common period of 3.75 d, leading the authors to conclude that these variations are produced by an inclined circumstellar disk that is warped close to the star, where it interacts with the stellar magnetosphere, and that co-rotates with the star.
ALMA observations of RY Lup in the 890-µm dust continuum (Ansdell et al. 2016) provided the first direct indication that the disk of RY Lup is seen at high inclination.A ring of dust is resolved around RY Lup, surrounding an inner cavity with a diameter of 0.8" (radius ∼ 60 au).The results of the UV-to-NIR study performed by (Arulanantham et al. 2018) to probe the details of the cavity and inner disk are fully consistent with the presence of a gas gap within the millimeter(mm)-dust cavity, confirming the pre-transition disk nature of the system.The inclination and the position angle of the disk derived from these ALMA data are 67 • and 109 • , respectively (Francis & van der Marel 2020).The accuracy on these angles was estimated to be 5 • by van der Marel et al. (2018).The disk dust mass M dust derived from the ALMA 1.3 mm data is 3.7 10 −4 M and the disk gas mass is 7.1 M Jup (Ansdell et al. 2018), which is equal to 6.8 10 −3 M .Using the ALMA 13 CO and C 18 O spectra, Yen et al. (2018) derived a stellar mass of 1.3 ± 0.1 M .This is compatible with the estimation of 1.47 ± 0.22 M by Alcalá et al. (2017).Yen et al. (2018) also determined the disk orientation through the velocity-aligned stacking method and obtained an inclination of 55 ± 5 • and a position angle of 110 +5 −10 deg.Langlois et al. (2018) presented a scattered light image obtained at H-band with the VLT/SPHERE instrument.The high inclination of the disk is confirmed, and interestingly the surface brightness of the disk reveals no inner gap, as opposed to the ALMA data.This is also a frequent signature of transition disks where the central part of the disk shows a decrease in density but is not completely void of gas and dust, especially small dust best seen in scattered light.This explains why the transition disk nature of the disk was not recognized in early SED studies (e.g., Manset et al. 2009).
We used GRAVITY guaranteed time observations of RY Lup to further constrain the inner geometry of the circumstellar environment of RY Lup.We present our observations in Sect. 2 and our modeling approach in Sect.3. The results of our modeling are detailed in Sect. 4 and discussed in Sect. 5.

GRAVITY
Within the context of the Young Stellar Objects Large Program of GRAVITY (Gravity Collaboration et al. 2017;Eisenhauer et al. 2011), we observed RY Lup with the Unit Telescopes (UTs) of the VLTI on 2017 June 11 for 2 hours.With the science channels (SCs), we recorded seven data sets at high spectral resolution (R ∼ 4000) corresponding to ten exposures, each of 30 seconds integration time.With the fringe tracker (FT), we recorded low-spectral-resolution data (5 channels across the K-band) at a   1) frame rate of about 1 kHz (Lacour et al. 2019).We interlaced our observations with sky exposures and observations of an interferometric calibrator, HD 110878, to retrieve the instrumental transfer function.We used the GRAVITY standard pipeline (Lapeyrere et al. 2014) to reduce and calibrate the observations.In this paper, we focus on the low spectral resolution K-band data provided by the FT for probing the inner dust rim in the K-band continuum.The high cadence of the FT frames allows the atmosphere to be frozen during each exposure of 0.85 ms, reducing the smearing effect on the interferometric fringes that can drastically degrade the visibilities.The FT calibrated data and the (u,v)-plane coverage of our observations are displayed in the top penal of Fig. 1.RY Lup is partially resolved with the UT baselines, with squared visibilities ranging from 0.55 to 0.8, and we detect no clear departure from centro-symmetry since the closure phases are all in agreement with zero.

Complementary data
VLTI/PIONIER RY Lup was observed in the H-band with the PIONIER instrument of the VLTI (Le Bouquin et al. 2011).We used two datasets obtained with two different configurations of the auxiliary telescopes (ATs; Table 1).The data were reduced with the standard PIONIER pipeline.RY Lup is partially resolved in the H-band with visibilities squared between 0.4 and 0.9 (Fig. 1-bottom).The visibilities quickly decline at the base-lines between zero to a few meters, which is generally interpreted as the contribution of a component that is much more extended than the angular resolution of the interferometer.Such a contribution can be interpreted as scattered light, as described in detail in Pinte et al. (2008a) and as applied in Anthonioz et al. (2015) who derive a ratio of 0.41 between the scattered light flux and the stellar flux for RY Lup in their composite model (see their Table 5 and Fig. 2).
Photometry We use nonsimultaneous photometric data from the literature to build the SED (see Table 2) as well as the Spitzer/IRS spectrum (Chen et al. 2016) to cover the mid-IR (MIR) and the silicate features at 10 and 20 microns.In order to minimise the impact of the known stellar variability, we use the X-Shooter spectrum published by Alcalá et al. (2017) to provide data over a wide range, from the B-to the K-bands.

A transition disk model for RY Lup
We used the radiative transfer code MCFOST to model the SED, the NIR, and sub-mm images, and the NIR H-and K-band interferometric data of RY Lupi.MCFOST is a 3D continuum radiative transfer code based on the Monte Carlo method (see Pinte et al. 2006Pinte et al. , 2009, for a complete description and benchmarking).MCFOST includes multiple anisotropic scattering when-Fig.2. Spectral energy distributions corresponding to our models of RY Lup: M50 for an inclination of 50 • with r in = 0.12 au (black solid line) and r in = 0.17 au (red dashed line), and M70 for an inclination of 70 • (black dotted-line) compared with the photometric measurements (orange circles; Table 2), the X-Shooter spectrum (blue crosses), and the Spitzer/IRS spectrum (red plus).The green lines display the SEDs calculated with the disk parameters of M50 (dashed line) and of M70 (dotted line) models but using the stellar parameters (T eff , L * ), and the inclination of 75 • used by Langlois et al. (2018).See text for details.
(W.m − 2) 3.40 1.33 10 −12 2.610 −13  Johnson L (a) 5.03 7.39 10 −13 1.510 −13  Johnson M (a) 11.5 3.87 10 −13 2.310 ever needed with a complete treatment of polarization.It also includes passive dust heating assuming radiative equilibrium, and continuum thermal re-emission.Viscous heating was not considered in our calculations.Very briefly, the code first computes the temperature structure assuming that dust is in radiative equilibrium with the local radiation field.This is done by propagating photon packets, originally emitted by the central star, through a combination of scattering (following Mie theory), absorption, and re-emission events until they exit the computation grid.
We used a parametric disk model for our calculations.The disk geometry includes a surface density Σ (r) = Σ 0 (r/r 0 ) −p and a scale height h (r) = H 0 (r/r 0 ) β prescription, with r 0 = 100 au.During calculations, we ensured that the dust temperature at the inner radius does not exceed the sublimation temperature, T sub , which we set at 1500 K. To achieve this we set the position of the inner disk, iteratively if needed, to ensure that the maximum dust temperature does not exceed T sub .We considered dust grains with sizes distributed according to the power law dn (a) ∝ a −3.5 da and we assumed two different grain size populations, that is, small grains (whose diameters a vary between 0.01 and 1.1 µm) and large grains (a = 10-1000 µm).For simplicity the grains are made of "astronomical silicate" and we used the optical properties published by Draine & Lee (1984).
In our calculations we assumed that the dust is at local thermodynamic equilibrium, and in radiative equilibrium with the stellar radiation field: at a given position, grains of all sizes have the same temperature, T dust .For the photosphere, we adopted a NEXTGEN stellar atmosphere model (Hauschildt et al. 1999) with an effective temperature of T eff = 4800 K and a surface gravity of log g =4.0.
We distributed the two grain populations in two distinct zones in the disk to describe a typical geometry for transition disks.This choice is guided by the ALMA images showing a ring of emission starting at ∼50 au and extending out to ∼125 au.This mm emission is dominated by large grains (zone #2 in our model).On the contrary, the SPHERE NIR image is dominated by scattering and small grains, and shows a disk extending all the way to the center (zone #1 in our model).Furthermore, following the results obtained for HL Tau by Pinte et al. (2016) and for a sample of edge-on disks (Villenave et al. 2020, submitted), we distributed the small dust in a flared disk with a standard aspect ratio in zone #1 (h/r ∼ 0.2) but we assumed significant vertical settling for the larger dust in zone #2.The dust in this zone (i.e., the ALMA ring) is therefore more concentrated toward the midplane with h/r ∼ 0.03.
Disk modeling can be a large multi-dimensional effort where finding the best-fitting, unique solution is CPU-expensive.Ideally, one would like to identify the best and unique solution to best fit all the data available.For RY Lup, that would be the SED, the SPHERE and ALMA images (NIR and mm), the IRS spectrum of the silicate features, and the NIR interferometric visibilities at H-and K-bands for example.However, as demonstrated by Pinte et al. (2008b, see Fig. 9 and Tab. 5) for example, this can only be achieved by running a very significant number of models.Running models of any individual sets of constraints, that is, the SED alone or the images alone, systematically provides solutions of lesser quality, with broader and shallower probability distribution functions (Pinte et al. 2008b, for the case of IM Lupi).It is also relevant to stress here that additional constraints, (1) at a reference radius of 100 au. (2)From north to east.coming a posteriori, can identify different or better solutions that were not indicated by the original modeling effort.This is the case for the models presented in Pinte et al. (2008b) for example, where a posteriori CO observations (Panić et al. 2009) revealed that the gas disk radius is at least twice that found by Pinte et al. (2008b) based on continuum data only.Therefore, we considered here that finding the ultimate model solution to match all the data available simultaneously was beyond the scope of the paper, as the paper focuses primarily on the new results provided by the GRAVITY instrument.We have therefore chosen a simpler approach.First we identified a model that correctly reproduces the SED.The adjustment is made by visual inspection and is not the result of a formal minimization.We focus in particular on correctly matching the optical/NIR part as this is the wavelength range where, in a second step, we calculate images and visibilities to compare with the VLTI interferometric data.The results from the SED fitting are presented in §4.1 below.To model the visibilities, the SED model identified in the first step is then adjusted, and modified as minimally as possible, to match the NIR interferometric data while maintaining a good match to the SED.The results are presented in §4.2.

Spectral energy distributiion, circumstellar extinction, and revised stellar luminosity
An SED was constructed using the photometric data listed in Table 2.We visually tested several models by spanning different ranges of stellar luminosity (L * ), disk dust mass, flaring index (β), scale height (H 0 ), and inclination (i) to broadly match the shape of the SED.The explored ranges are given in the last column of Table 3. Well-defined parameters were kept fixed at values derived from previous studies.Interstellar extinction (A v ) was applied a posteriori to the modeled SED.This is in addition to the circumstellar extinction naturally caused by the disk and fully taken into account during calculations.In this first modeling step, to match the SED only, we fixed the inner radius at 0.17 au.The disk outer radius R out is set fixed to 125 au (Ansdell et al. 2018), and the total dust mass M dust is set to 3.7 10 −4 M (the sum of the dust mass in the two model disk zones), which provides a good match to the SED and the mm-flux.We set the position angle from ALMA observations to 109 The recent scattered light images presented by Langlois et al. (2018) suggest that the disk reduces the apparent stellar luminosity because of circumstellar extinction, by occulting the photosphere.In our calculations, the luminosity was therefore left as a free parameter.We found that in order to match the observed flux density, a luminosity of ∼ 2.5 L is needed.This is a factor of ∼ 1.5 higher than the value of 1.7 L reported by Alcalá et al. (2017).Our revised luminosity translates into a stellar radius of R = 2.29 R , when considering an effective temperature of 4800 K.The circumstellar extinction is calculated self-consistently by the radiative transfer code, which includes the effects of anisotropic and multiple scattering.An additional small amount of interstellar extinction, A v = 0.7 mag, is needed to match the shape of the SED in the optical/NIR (with R v = 3.1).A model that reproduces the overall shape of the SED reasonably well is shown in Figure 2. Its parameters are given in Table 3 in the column entitled M50 (the inclination of this model is 50 • ).As we show in the following section, this inclination and model M50 produce the best agreement with the GRAVITY data.
For comparison, we overplotted the SED computed for our disk model M50, which has the same geometry and interstellar extinction, but takes the inclination, luminosity, and effective temperature used by Langlois et al. (2018) following Manset et al. (2009), that is 50 • , L * = 2.4 L (calculated from a stellar radius of 1.72 R ), and an effective temperature of 5500 K. Clearly, the resulting SED (displayed by the green dashed line in fig.2) is globally too faint to match the photometric data.This is because the model by Langlois et al. (2018) does not include an inner disk.The disk used by these latter authors is a single flat dust ring with a large inner cavity extending out to r = 18 au, as suggested by ALMA (Ansdell et al. 2016).In that case the photosphere is not occulted by the disk, even though the authors are considering a high inclination of 75 • , contrary to a case where the inner disk would be located much closer to the star, as in our M50 model.
However, since the SPHERE and ALMA images of the outer disk suggest an inclination of about 70 • (Langlois et al. 2018;Francis & van der Marel 2020), we also investigated the parameter space needed to reproduce the SED at such a high inclination.The SED of this model, that we call M70, is overplotted in black dotted line in Fig. 2. Its parameters are given in the column M70 of Table 3.The overall agreement with the SED is similar to, though slightly poorer than, the one provided by the M50 model (see the values of reduced χ 2 for the SED fits in Table 3).Of importance, the flux density in the GRAVITY range is correctly reproduced.Due to the higher inclination of the system, the disk blocks more flux from the star compared to the M50 model.Accordingly, the stellar luminosity of the M70 model has to be higher to match the SED data.We find a luminosity of 4.65 L , which translates into a stellar radius of R = 3.13 R for an effective temperature of 4800 K.

Visibility fitting
After identifying a disk model that matches the SED reasonably (Fig. 2), we added one layer of complexity by calculating images and visibilities to verify the capacity of that model to also fit the VLTI interferometric data.We first computed intensity maps in the H-and K-bands with a pixel sampling of 0.05 mas, as illustrated in the top panel of Figure 3 for the K-band.This sampling is about 70 times denser than our K-band spatial resolution and 55 times denser than our H-band spatial resolution.These images were then processed by Fourier transform to produce interferometric visibilities.We used the ASPRO2 tool 1 provided by JMMC to compute the visibilities with the exact same (u,v)coverage as the observations and to take into account the field of view of the telescopes through a convolution with a Gaussian whose half width at half maximum (HWHM) is equal to the point spread function of the telescopes (i.e., 60 mas for the UTs and 250 mas for the ATs).
For model M50, we needed to slightly adjust the inner disk radius to fine tune the adjustment of the visibility curves, while maintaining a proper match to the SED.An inclined disk (with an inclination of ∼50 • ) and an inner dust rim R in at 0.12 au correctly fits the SED and the NIR interferometric observations (Fig. 3).We note that reducing the inner radius to 0.12 au entails an increase in the maximum dust temperature at the inner disk up to 1700 K, slightly beyond to adopted sublimation tempera-1 Available at http://www.jmmc.fr/asproture of 1500 K, but still reasonable.The fit of the high-quality GRAVITY visibilities is excellent.Our model is also in agreement with the PIONIER data, considering the larger error bars of these datasets, even if the difference between the data and the model is larger.Such a difference could come from the variability of the source because the geometry of the inner disk is expected to be very similar in both the H-and K-bands.Contemporary interferometric datasets are needed to investigate this in more detail.We also verified that our M50 model leads to closure phase signals that are consistent with ∼ 0 • , as observed by both GRAVITY and PIONIER (Fig. 3).This is the case because even if the central part of the disk model is not centro-symmetric, its compactness (1 mas or less) is such that it is only partially resolved by GRAVITY as shown by the squared visibilities higher than 0.55, even at the longest baselines.The parameters of our final model M50 are given in Table 3.
For the M70 model, we failed to correctly reproduce the GRAVITY visibilities (Fig 4), even when changing the inner rim position which is about three times further away from the star than for the M50 model due to the higher luminosity.In particular, the modeled visibilities along three baselines exhibiting a position angle of 150-160 • (in green, orange, and yellow in Fig. 1) remain systematically larger than those measured during observations.See the values of reduced χ 2 for the visibility fits in Table 3).

Discussion and concluding remarks
Taking into account the dimming of the photospheric flux caused by the extinction from the inclined disk, we derive a revised luminosity of 2.5 L for an inclination of 50 • (M50).For comparison, Alcalá et al. (2017) report a luminosity of 1.7 L for a similar effective temperature.With the effective temperature and luminosity we estimated for M50, RY Lup registers in the Hertzsprung-Russell diagram immediately outside the range calculated by Baraffe et al. (2015) (see Fig. 6 of Frasca et al. (2017), in-between the 1 and 3 Ma isochrones).These latter calculations stop at 1.4 M , implying a small extrapolation, but such an extrapolation would indicate a mass slightly above 1.4 M and an age of ∼ 2 Ma for RY Lup for M50.Using the Siess et al. (2000) tracks for Lupus, this would correspond to a ∼ 1.6 M and 3-4 Ma pre-main sequence star (see, e.g., Fig. 2 of Alcalá et al. (2017)).A comparison with older pre-main sequence tracks by D' Antona & Mazzitelli (1994) converts into a ∼ 1.3 M and ∼ 1 Ma pre-main sequence object, using their CM convection models In all cases, this brings the age of RY Lup (about 2.0 Ma) more in line with the age of the Lupus association, in opposition to the estimations of for example 10.2 Ma by Frasca et al. (2017) and 12 Ma by Manset et al. (2009).This discrepancy can be explained by the systematic underestimation of the luminosity caused by not taking into account the circumstellar extinction due to the disk, and only assuming a mild interstellar extinction of the photosphere in agreement with the observed photospheric colors.The underluminous nature of highly inclined disk systems in Lupus is discussed in Frasca et al. (2017), although not specifically for RY Lup.Other similar cases of underluminous, highly inclined disk systems in Lupus, whose ages are drastically higher than the age of the Lupus association, include MY Lupi, Sz 112, and Par-Lup 3-4 (Alcalá et al. 2017).
For M50, the evolutionary models mentioned above also predict masses that are in agreement with observations.The estimates range from 1.47 ± 0.22 M to 1.3 ± 0.1 M (Alcalá et al. 2017;Yen et al. 2018, respectively), the latter being a "dynamical mass" estimated from the observed rotation of the gas disk.For the model with an inclination of 70 • (M70), the derived luminosity is higher (4.65 L ) and corresponds to a mass larger than 2 M using the Siess tracks, which is well above the observed estimates.At this stage, we need to recall that our favored model for the system (M50) does not include a detailed modeling of the scattered light image, and in particular the solution for the disk scale height is driven by the SED, and it is therefore possible (if not probable) that the stellar luminosity we derived, although improved with respect to previous estimates, still requires refinement based on a much more thorough analysis of all the data available.Given the intrinsic large variability of the star, such refinement is well beyond the scope of this paper.
Within the range of parameters explored in Table 3, we could not find another disk model that was able to reproduce both the SED and the visibilities properly.Indeed, many parameters are well constrained either by the SED (e.g., the stellar luminosity and the extinction, the disk dust mass, and flaring) and/or by the interferometric data (the inner radius).
It is worth comparing the inclination of the dusty disk we derive to match the SED and interferometric data with the largescale disk images provided by SPHERE and ALMA.Our models with inclinations larger than 50 • produce too much extinction in the SED.More extinction would force the stellar luminosity to higher values than 2.5 L , making the object more massive and in tension with current mass and age estimates, as demonstrated with the M70 model.In all of the models presented above the inner disk is sharp, with a clear density truncation in the radial direction at the inner radius.We explored a more realistic geometry, although parametric, for the inner disk by making the inner wall rounder (larger R in with increasing |z| distance).The jump in density was also made smoother by using a steep exponential taper in the radial direction at the inner radius.Even with this parametrization of the density at the inner rim we failed to correctly fit the visibility curve at all spatial frequencies.Our model considers only a single inclination for the inner and outer parts of the circumstellar disk for now, and we have not yet investigated a potential warp or misalignment to reconcile the observed inclinations of the inner and outer disks.To go further, it would be necessary to improve the details of the vertical structure of the disk.A vertically thinner disk would allow for example to accommodate a larger inclination while maintaining the circumstellar extinction at the current level.Langlois et al. (2018) reports the presence of spiral arms in the H-band scattered light image.If a planet is driving these spirals, then the radial and vertical disk structure can be dynamically altered by the forming planetary companion, leading among other effects to a misalignment between the inner and outer disks that could also reconcile the inclination measurements of the inner and outer disks.Such a misalignment is also invoked to explain the shadows observed in the SPHERE images (Marino et al. 2015).Complementary high-angular-resolution techniques covering different spectral ranges (e.g., GRAVITY, SPHERE, ALMA) allow such inclination differences in protoplanetary disks to be to investigated provided the interferometric (u,v)-coverage is high enough to produce a high-fidelity image of the inner disk, which is not the case for our observations of RY Lup.Within this context, the sensitivity gain brought by the NAOMI adaptive optics systems on the ATs (Woillez et al. 2019) is important for expanding the (u,v)-coverage and also for improving the measurement accuracy to look for variability and/or to better constrain the scattering phenomena that can be investigated through interferometric observations at very short baselines.Additional and complementary data, also including the interferometric instrument MATISSE (Lopez et al. 2018) operating in the MIR range, and probing the intermediate scales between the very innermost regions (< 1 au) probed by NIR interferometry and the few tens of astronomical units probed by ALMA, would allow us to constrain the more complex features of this dusty disk.For this purpose, due to the variability of the source, a multi-technique campaign of simultaneous observations would be required.

Fig. 3 .
Fig. 3. Top.Intensity maps produced by MCFOST for our model M50 of RY Lup in the K band for the whole disk (left) and zoom in at the milliarcsecond scale of the central region (right).Middle.Observed squared visibilities (red circles) with GRAVITY (left) and PIONIER (right) compared with the visibilities computed with our model M50 (black circles).Bottom.Observed closure phases (red circles) with GRAVITY (left) and PIONIER (right) compared with the closure phases computed with our model M50 (black circles).
A disk inclination of 50 • is in agreement with the inclination determined from the ALMA CO spectra analysed by Yen et al. (2018) (55 ± 5 • ).However, the inclinations derived from the aspect ratio of the large-scale outer disk imaged by SPHERE (70 • ; Langlois et al. 2018) and ALMA (67 ± 5 • ; Francis & van der Marel 2020) are marginally discrepant with our model.Finally, the comparison of the GRAV-ITY interferometric data with model M70 clearly shows that this model is unable to reproduce the visibility curve in detail, in particular for the three baselines that directly probe the front of the inner rim due to their orientations.

Fig. 4 .
Fig. 4. Squared visibilities (left) and closure phases (right) as observed with GRAVITY (red circles) and as computed with our model M70 (black circles).

Table 2 .
Photometric data used to build the SED of RY Lup.

Table 3 .
Star and disk parameters of our models of RY Lup , M50, and M70.See main text for the definitions and the corresponding references for the fixed values.