© GRAVITY Collaboration 2021

1 Introduction

The formation and evolution of protoplanetary disks are directly linked to planet formation. The outer disk features of young stellar objects (YSOs) have been thoroughly studied in the past through scattered light imaging (e.g., with SPHERE/VLT, Beuzit et al. 2019) and with ALMA in the (sub-)millimeter range (ALMA Partnership 2015). With both techniques, disk observations have shown rings, gaps, and asymmetric structures up to a few hundred au (e.g., Lodato et al. 2019; Benisty et al. 2017). The inner regions (at ~au scale) of such disks are also of prime interest since key processes like gas accretion flows, winds, outflows, and dust sublimation take place. All these processes affect the dynamics and evolution of the first few au regions where terrestrial planets may form and/or migrate over few million years. Constraints on these processes can be derived indirectly through spectroscopic studies, but at typical distances of a few hundred parsecs only observations with milliarcsecond (mas) resolution, which are required to probe sub-au scales, can discriminate between competing models. Numerous interferometric studies have been conducted in the past in the near- and the mid-infrared (IR), for instance with IOTA (Millan-Gabet et al. 2001), the Palomar Testbed Interferometer (PTI; Eisner et al. 2004), the Keck Interferometer (Monnier et al. 2005; Eisner et al. 2014), and the VLTI (Menu et al. 2015; Lazareff et al. 2017; GRAVITY Collaboration 2019). To date, the new generation of four-telescope instruments including GRAVITY, operating in the K band (GRAVITY Collaboration 2017), and MATISSE, operating in L to N band (Lopez et al. 2014), are pushing further the achievable spectral coverage, the sensitivity, and precision of interferometric measurements. Even though statistical studies of large YSO samples are of high relevance (Lazareff et al. 2017; GRAVITY Collaboration 2019), some of these objects require a more in-depth study enabled by the improved data quality of recent interferometric observations. Some objects, like HD 141569, are unusual in their evolutionary sequence and require dedicated studies. For such systems the exact nature and properties of the very inner regions are still a matter of debate and can be solved in part by modeling the distribution of the warm dust traced in the K band. Furthermore, it is still unclear how to characterize the star–disk mechanisms traced by emission lines like the hydrogen Brackett-γ (Brγ ~ 8000–10 000 K) and CO (T ~ 2000–3000 K) bandheads, and where this emission occurs. The high-quality spectroscopic capabilities of the GRAVITY instrument allows us to study in detail the gas phase and its spatial morphology thanks to the interferometric visibilities and differential phase signals (GRAVITY Collaboration 2020a,b, 2021).

HD 141569 is a Herbig star classified as a B9-A0 spectral type (Augereau & Papaloizou 2004), with an effective temperature of 9750 ± 250 K, an estimated age of 7.2 ± 0.02 Myr, a luminosity between 16.60 ± 1.07 L_⊙ (Vioque et al. 2018) and 27.0 ± 3.6 L_⊙ (Di Folco et al. 2020), a mass of 2.14 ± 0.01M_⊙, and a Gaia distance of 110 ± 1 pc (Arun et al. 2019)². It is a non-flaring disk system with little mid-IR excess classified as a group-II source (Meeus et al. 2001). It is the only known pre-main sequence star characterized by a hybrid disk (Wyatt et al. 2015; Péricaud et al. 2017; Di Folco et al. 2020), an evolutionary disk state between the protoplanetary and debris-disk regimes. Near-IR imaging spatially resolved an optically thin disk consisting of two rings located at about ~ 280 and ~ 455 au from the star Augereau et al. (1999); Weinberger et al. (1999); Biller et al. (2015). A more complex system is shown in the visible, consisting of multiple rings and outer spirals that could be explained through perturbations by two nearby (~ 7.5 arcsec) M dwarfs, or by planetary perturbations (Augereau & Papaloizou 2004; Wyatt 2005; Reche et al. 2009). Fisher et al. (2000) found a warm disk component up to 110 au at 10.8 and 18.2 μm, later confirmed by Marsh et al. (2002). The short-wavelength counterpart of this component was detected by Mawet et al. (2017) through L′ imaging and ranging between 20 and 85 au. Emission at 8.6 μm was detected by Thi et al. (2014), and was interpreted as emission from polycyclic aromatic hydrocarbons (PAHs). NOEMA and ALMA observations in the millimeter range showed continuum emission equally shared between a compact (≲ 50 au) and a smooth extended dust component (~350 au), with large millimeter grains dominating the inner regions and smaller grains in the outer ones (Di Folco et al. 2020). Finally, inner disk features were detected by SPHERE in the Y, J, H, and K bands (Perrot et al. 2016) and by Keck/NIRC2 in the L′ band (Currie et al. 2016) at physical separations of 45, 61, and 88 au. These results point out the high morphological complexity of the outer disk in the HD 141569 system.

Little is known about the central astronomical units of the system. The spectral energy distribution (SED) of HD 141569 alone does not help us in this sense since the IR excess is very small (see Fig. 9 of Thi et al. 2014). The majority of the K-band measurements listed in Table G.1 reflect a featureless SED in the near-IR. Moreover, pure SED fits obtained by different authors (Li & Lunine 2003; Merín et al. 2004; Thi et al. 2014) may suggest at first that the near-IR emission is exclusively photospheric in nature and that the disk contributes only at longer wavelengths. The object was observed at milliarcsecond resolution in the K band with the PTI and the Keck interferometer, but it was spatially unresolved (Eisner et al. 2004, 2009). Monnier et al. (2005) derived a 10 mas upper limit in radius for the spatial extension of the K-band emission. Therefore, trustworthy information on the first 5 au of the system are scarce, and the question arises of whether the inner region of the disk could be already in a debris disk stage, where the SED fits in the near-IR are not accurate enough to detect such a faint excess.

The circumstellar gas has been observed in both atomic and molecular form, which suggests the system has not yet reached the gas-depleted stage characteristic of a debris disk system. Mendigutía et al. (2017) set an upper limit of ~ 0.11 au for the gas region responsible for the spatially unresolved double-peaked Hα emission. A comparable upper limit of ~0.13 au for the Brγ line emitting region is suggested by Eisner et al. (2009). Both lines are observed to be not variable over timescales of days and years (Eisner et al. 2015; Mendigutía et al. 2011b). In addition to hydrogen, CO ro-vibrational emission (v ≥ 1, Δv = 1) was observed by many authors extending from 10 to 275 au (Dent et al. 2005; Goto et al. 2006; Brittain et al. 2007; Flaherty et al. 2016; White et al. 2016; Miley et al. 2018; Di Folco et al. 2020).

We present here the first GRAVITY interferometric observations of this disk, with the goal of revealing the geometry and dynamics of the internal structure of HD 141569, and of gaining insights about the dust and gas properties. Section 2 describes the observations; Sects. 3 and 4present the observational data and the adopted methodology; Sect. 5 describes the results of the possible scenarios along with the corresponding modeling; a discussion is developed in Sect. 6.

2 Observations

HD 141569 was observed with VLTI/GRAVITY (GRAVITY Collaboration 2017) using the four 1.8 m Auxiliary Telescopes (ATs), on March 18 and July 12, 2019, in the intermediate D0-G2-J3-K0 configuration, and on May 23, 2019, in the large A0-G1-J2-J3 configuration. The observations span a spatial frequency range between about 15 and 65 Mλ (see Fig. 1, left panel) with a maximum angular resolution of λ∕2B, about 1.7 mas for the longest baseline (B) of 130 m, which corresponds to about 0.19 au at a distance of 110 pc. The data consist in high spectral resolution (R ~ 4000) observables recorded by the science channel (SC) detector over the whole K band with individual integration times of 30 s and in low spectral resolution (R ~ 20) observables recorded by the fringe tracker (FT) detector (five spectral channels over the K band at 1.908, 2.058, 2.153, 2.256, 2.336 μm) at frame rates of ≈300 and ≈900 Hz (Lacour et al. 2019). Each observation block corresponds to 5 min on the object. In total, three files were acquired in March 2019, one inMay 2019, and eight in July 2019. HD 141569 observations were preceded by the observation of a point-source calibration star, close to our object on the sky, and of similar spectral type and brightness in order to calibrate the atmospheric and instrumental transfer function. We observed the calibrator HD 137006 in March and July 2019, and the calibrator HD 141977 in May 2019. Further details on the executed runs and observing conditions are given in Table A.1.

Fig. 1 HD 141569 FT data; squared visibilities (left panel), closure phases (central panel), and U-V plane coverage (right panel), fromall the observation epochs. Colors refer to the different GRAVITY spectral channels.

3 Data

All the data were reduced and calibrated using the GRAVITY data reduction software (Lapeyrere et al. 2014). For the low-resolution FT data we discarded the first spectral channel, which can typically be affected by the metrology laser operating at 1.908 μm. Figure 1 shows the U-V plane coverage and the FT calibrated squared visibilities and closure phases (right, left, and center panel, respectively). Following GRAVITY Collaboration (2019), we applied a floor value on the error bars of 2% for the squared visibilities and 1° on the closure phases as the error bars computed by the pipeline might be underestimated or correlated. We observe that GRAVITY partially resolved the near-IR emission in HD 141569 with squared visibilities between 0.8 and 1.0. Therefore, the data can be used to estimate the characteristic size of the dust environment (see Sect. 5.1). Moreover, with the inclusion of the May 2019 large configuration data, we observe that the visibility reaches a plateau at almost all spatial frequencies, allowing us to constrain near-IR flux contributions of the star and environment.

Since the closure phases are consistent with 0° at all baselines and for all the epochs, we can confidently consider the emission to be centro-symmetric on the spatial scale of our observations, and we therefore discard the hypothesis of a close companion as the origin of the resolved emission, at least within the 250 mas (~ 28 au) field of viewof GRAVITY with the ATs.

For the high-resolution SC data, we concentrated on the July 2019 dataset only since this is the epoch where we gatheredthe highest number of files. In order to optimally exploit the SC dataset, the eight files from July 2019 have been merged in order to increase the signal-to-noise ratio per spectral channel. Considering that the maximum span in position angle between the two extreme positions of the UV coverage is only ~4°, we do not expect any visibility smearing of the data due to differences in hour angles. The error bars were computed as the standard deviation between the eight files of the corresponding differential quantity (visibility or differential phase) in each spectral channel. For instance, we derived the absolute error on the visibility to about 1.8%. Figure 2 shows the visibility amplitude (left panels) and the differential phase (right panels) for the six baselines in the region of interest of the Brγ line, between 2.15 and 2.18 μm. The top panels show the object spectrum normalized to the continuum, corrected for telluric lines (left plot), and corrected for both telluriclines and photospheric absorption (right plot). The visibilities appear to be spectrally flat with no clear signature at 2.16612 μm. They are measured to vary between 0.9 and 0.96 as a function of the baseline, which is indicative of a compact region well inside the dusty disk. Interestingly, the differential phase signal is more marked at the position of the Brγ line. We observe a clear S-shaped signal through baselines J3-D0 (100.1 m, 220°) and J3-K0 (53.6 m, 152°), and a weaker signal through baselines J3-G2 (62.1 m, 223°) and D0-G2 (38.2 m, 36°). No clear signal is detected within the error bars through the baselines D0-K0 (95.1 m, 72°) and G2-K0 (68.1 m, 92°).

In order to study the Brγ line gas region using the SC high spectral resolution data, we need precise measurements for the line-to-continuum flux ratio. It is essential to perform a proper wavelength calibration of the spectrum and to take out the contribution by the telluric lines. We describe the whole procedure in Appendix C. Errors are from the original data, reduced and calibrated through the GRAVITY data reduction software.

The high-resolution GRAVITY spectrum of HD 141569 (see top panels of Fig. 2) shows a double-peaked Brγ emission line. Since the error associated with the GRAVITY wavelength calibration is ~ 3 Å, the two peak positions can be considered to be symmetric with respect to the Brγ wavelength rest position. Both the double-peaked emission line and the S-shaped feature in the differential phase suggest a scenario where the gas emitting in the Brγ line could be in Keplerian rotation. We explore this hypothesis further in the following sections.

Fig. 2 Science channel data of July 2019. Top left: wavelength-calibrated continuum-normalized spectrum,corrected for telluric lines. Top right: same as top left, but corrected for the Brγ photospheric absorption. Differential phases (left column) and visibilities (right column) along the six GRAVITY baselines. The red lines in the visibility plots show the pure-line visibilities.

4 Methodology

The properties of the spatially resolved continuum emission were first investigated with the help of low spectral resolution FT data (see Sects. 4.1 and 5.1). The squared visibility curve was modeled through chromatic geometrical models accounting for a point-like central star, and simple geometrical shapes (Gaussian geometrically thin rings, Gaussian-convolved infinitesimally thin rings) representing the circumstellar environment. Useful information was obtained such as the star-to-dust flux ratio, the dust spectral index, and the spatial distribution of the dust. To gain further information on the dust emission properties, we investigated through radiative transfer (RT) modeling the impact of such a component on the SED (see Sect. 5.2). We used the RT code MCMax (Min et al. 2009), which solves 2D RT (e.g., Bjorkman & Wood 2001) to calculate the dust density and temperature structure of a given disk setup. In Sects. 4.2 and 5.3 we discuss the high-spectral resolution SC data in the Brγ region used to constrain the spatial scale of the hot gas emitting component through modeling of the visibility curves, and through the analysis of the Brγ spectrum under the assumption of a gas disk in Keplerian rotation. Further information on the gas region size and its dynamical properties were derived through the analysis of the differential phases and the resulting photocenters shifts. Finally, an analytical axisymmetric Keplerian disk model is compared to our observations.

4.1 Dust continuum: low spectral resolution data

Following the work of Lazareff et al. (2017) and GRAVITY Collaboration (2019), we used geometric models that consist of a point-like central star, assumed to be unresolved at all observed baselines, and a circumstellar environment in order to fit the observed visibilities. The complex visibility of the system at spatial frequencies (u, v) and at wavelength λ is therefore described by a linear combination of the two components as $$V(u,v,\lambda )=\frac{F_{\text{s}}(\lambda )V_{\text{s}}(u,v,\lambda )+F_{\text{c}}(\lambda )V_{\text{c}}(u,v,\lambda )}{F_{\text{s}}(\lambda )+F_{\text{c}}(\lambda )},$$(1)

where V_c is the visibility of the circumstellar environment, and F_s and F_c the specific fractional flux contributions of the star and of the circumstellar environment, respectively. The visibility of the star V_s is equal to 1 since we assume it to be unresolved. Since our GRAVITY FT data contains six visibility measurements and four closure phases for each of the four spectral channels and for each file, we can derive the spectral dependence of the circumstellar environment by modeling it as a power law, defined by its spectral index k_c, where k = d log F_λ∕d log λ, and by describing the complex visibility of the system as $$V(u,v,\lambda )=\frac{F_{\text{s}}{(\lambda /{\lambda }_0)}^{k_{\text{s}}}+F_{\text{c}}{(\lambda /{\lambda }_0)}^{k_{\text{c}}}{V}_{\text{c}}(u,v)}{F_{\text{s}}{(\lambda /{\lambda }_0)}^{k_{\text{s}}}+F_{\text{c}}{(\lambda /{\lambda }_0)}^{k_{\text{c}}}},$$(2)

where λ₀ = 2.15 μm is the wavelength of the central spectral channel of the FT, and k_s the spectral index of the star derived assuming that it radiates as a black body at the star effective temperature T_eff = 9750 K, which translates into s spectral index k_s = −3.62 at λ₀ for the central star.

We chose to fit our visibility data with three different geometric models that differ only by the V_c term in Eq. (2). Since the closure phases are basically zero for every baseline, as shown in the central panel of Fig. 1, we do not consider any azimuthal modulation in our models. Therefore, the resulting brightness distributions are centro-symmetric.

The first model consists of a Gaussian disk whose visibility is described, from Berger & Segransan (2007), as $$V_{\text{c}}(u,v)=V_{\text{gauss}}(u,v)=\text{exp}\left(-\frac{{(\pi \Theta r)}^2}{4\text{ln}2}\right),$$(3)

where Θ is the Gaussian full width at half maximum (FWHM) and $r=\sqrt{u^2+v^2}=B/\lambda$, with u and v the spatial frequency coordinates, and B the projected baseline.

The second model consists of a geometrically thin ring whose visibility is described by subtracting an inner smaller uniform disk from a larger one: $$V_{\text{c}}(u,v)=V_{\text{ring}}(u,v)=\frac{F_{\text{outer}}{V}_{\text{outer}}-F_{\text{inner}}{V}_{\text{inner}}}{F_{\text{outer}}-F_{\text{inner}}}.$$(4)

Here the subscript outer refers to the larger disk and the subscript inner refers to the inner hole. After normalization and knowing that F_disk = π D²∕4, where D is the diameter of the disk, we can express Eq. (4) as $$V_{\text{ring}}(u,v)=\frac{D_{\text{outer}}2\frac{J_1(\pi D_{\text{outer}}r)}{\pi r}-D_{\text{inner}}2\frac{J_1(\pi D_{\text{inner}}r)}{\pi r}}{D_{\text{outer}}^2-D_{\text{inner}}^2},$$(5)

where J₁ is the first-order Bessel function. Finally, our last model consists of an infinitesimally thin ring convolved by a Gaussian whose visibility is described as the product between V_gauss(u, v) and $$V_{\text{inf-ring}}(u,v)=J_0(2\pi ar),$$(6)

where J₀ is the zero-order Bessel function and a is the radiusof the infinitesimally thin ring (Berger & Segransan 2007). Since in the Fourier space the convolution of two functions is simply their multiplication, we have $$V_{\text{c}}(u,v)=V_{\text{gauss-ring}}(u,v)=V_{\text{gauss}}(u,v)\cdot V_{\text{inf-ring}}(u,v).$$(7)

The inner hole radius and ring width are defined as r_i = D_inner∕2 and w = (D_outer − D_inner)∕2, respectively, for the geometrically thin ring model while they are defined as r_i = a − Θ∕2 and w = Θ, respectively, for the Gaussian-convolved infinitesimally thin ring. Inclination and position angle of the circumstellar environment are taken into account through the parameter $r=\sqrt{u^2+v^2}$ following Berger & Segransan (2007). The model fitting is based on a Markov chain Monte Carlo (MCMC, Foreman-Mackey et al. 2013) numerical approach and was implemented on the combined dataset of all three epochs in order to maximize the number of experimental points against the number of free parameters. This assumes that the near-IR emission and the disk structures are not variable over a five-month period, which is strengthened by the fact that the star is not variable, either spectroscopically in the optical (Mendigutía et al. 2011a) or photometrically in the mid-IR (Kóspál et al. 2012). Once a global solution was identified, we further checked how well the parameters are indeed constrained by the data. For this purpose, we performed a series of squared visibility fits to the model by fixing the tested parameter to different values and left the other parameters free in the subsequent minimization. In this way we obtain a χ² curve as a function of the tested parameter value.

Finally, the error on the χ² was estimated by treating the quantity T_i as a stochastic variable with $$T_i=N{\left(\frac{y_i-y_{\text{model}}}{{\sigma }_{\text{i}}}\right)}^2,$$(8)

where N is the number of points in the dataset, y_i the individual measurement, y_model the value of the model, and σ_i the error associated to the measurement. The χ² value is given by the mean of T_i, $${\chi }^2=\frac{1}{N}{\sum }_{i=1}^N{T}_i={\sum }_{i=1}^N{\left(\frac{y_i-y_{\text{model}}}{{\sigma }_i}\right)}^2,$$(9)

and the χ² error is given by the error on the mean of T_i, $${\sigma }_{{\chi }^2}=\sqrt{\frac{{\sum }_{i=1}^N{T}_i^2}{N(N-1)}-\frac{{({\chi }^2)}^2}{N-1}}.$$(10)

This applies to the reduced χ² as well.

4.2 Gas: high spectral resolution data

To estimate the gas region size from the SC visibility we extrapolated the pure-line contribution from the total visibility (line+continuum) displayed in Fig. 2. To do this we modeled the total visibility with a three-component model that accounts for the contributions from the star, the circumstellar dust, and the line emitting gas. The total visibility is given therefore by $$\begin{array}{lll}\hfill & \hfill & V_{\text{tot}}(u,v,\lambda )\hfill \\ \hfill & \hfill & =\text{}\frac{\alpha (\lambda )F_{\text{s}}(\lambda )V_{\text{s}}(u,v,\lambda )+F_{\text{c}}(\lambda )V_{\text{c}}(u,v,\lambda )+F_{\text{L}}(\lambda )V_{\text{L}}(u,v,\lambda )}{\alpha (\lambda )F_{\text{s}}(\lambda )+F_{\text{c}}(\lambda )+F_{\text{L}}(\lambda )},\hfill \end{array}$$(11)

where α(λ) is the science star continuum-normalized photospheric absorption (see Appendix C for more details), the subscript c refers to the dust component, and the subscript L refers to the Brγ line gas. From Eq. (11) it can be proven (see Appendix D) that the pure-line visibility is given by $$\begin{array}{lll}\hfill & \hfill & V_{\text{L}}(u,v,\lambda )\hfill \\ \hfill & \hfill & =\frac{V_{\text{tot}}(u,v,\lambda )\{[1+\beta (\lambda )]F_{\text{L/C}}-1-\beta (\lambda )\}-\alpha (\lambda )+1}{[1+\beta (\lambda )]F_{\text{L/C}}-\alpha (\lambda )-\beta (\lambda )},\hfill \end{array}$$(12)

where F_L/C is the line-to-continuum flux ratio $$F_{\text{L/C}}=\frac{\alpha (\lambda )F_{\text{s}}(\lambda )+F_{\text{c}}(\lambda )+F_{\text{L}}(\lambda )}{F_{\text{s}}(\lambda )+F_{\text{c}}(\lambda )}$$(13)

and β(λ) is the disk-to-star flux ratio outside the line $$\beta (\lambda )=\frac{F_{\text{c}}}{F_{\text{s}}}(\lambda )=\frac{F_{\text{c}}}{F_{\text{s}}}{\left(\frac{\lambda }{{\lambda }_0}\right)}^{k_{\text{c}}-k_{\text{s}}}.$$(14)

For clarity, we note that, outside the line emitting region, α(λ) and F_L (λ) tend toward 1 and 0, respectively. Equation (13) corresponds to the line-to-continuum ratio including the photospheric absorption. Finally, we estimated the gas region size by modeling it with an infinitesimal ring model given by Eq. (6). In the same way, we needed to take out the continuum contribution from the total differential phases (Fig. 2). Following Weigelt et al. (2011), the pure-line differential phase is given by $$\text{sin}{\varphi }_{\text{L}}=\text{sin}\varphi \frac{F_{\text{tot}}{V}_{\text{tot}}}{F_{\text{L}}{V}_{\text{L}}}=\text{sin}\varphi \frac{V_{\text{tot}}}{V_{\text{L}}}\frac{F_{\text{L/C}}}{F_{\text{L/C}}-\left(\frac{\alpha +\beta }{1+\beta }\right)},$$(15)

where ϕ_L is the pure-line differential phase, ϕ the total differential phase, F_tot the total flux (star, dust, and gas), and we write the ratio F_tot∕F_L through Eq. (D.5). Following Le Bouquin et al. (2009), we can derive wavelength-dependent photocenter displacements along each baseline from the pure-line differential phases by $$p=\frac{-{\varphi }_{\text{L}}}{2\pi }\cdot \frac{\lambda }{B},$$(16)

where p is the projection on the baseline B of the 2D photocenter vector with origin on the central star. The error bars on the pure-line differential quantities are computed through error propagation in Eqs. (12) and (15).

5 Results

5.1 Disk component inside 2 au

Of the three models discussed in Sect. 4.1, the fit of the squared visibilities to the Gaussian disk model (Eq. (3)) did not converge to a solution (i.e., the marginal posterior distributions for the Gaussian width, inclination, and position angle are flat). Therefore, we discarded this model in the rest of the work. Solutions were found for the geometrically thin ring model (Eq. (5)) and the Gaussian-convolved ring model (Eq. (7)). The two models converge basically toward the same solution, with a reduced ${\chi }_r^2$ = 4.7 for both models.

We performed a wide scan range of the fitted parameters to find convergence toward a global solution. The results of the minimization are presented in Fig. B.1 and Table 1 for the ring model along with the parameter scan range and the 1σ uncertainties. The resulting MCMC posterior distribution is presented in Fig. B.1 (top plots) and allows us to identify an optimal global solution for the six parameters. The fitting process leads to a photospheric near-IR flux contribution of ~ 93.8%, and therefore to a dust ring flux contribution of ~6.2% for both models. Interestingly, the degeneracy typically found between the disk’s flux and the characteristic size in V² is broken here because the constant plateau as a function of spatial frequencies unambiguously determines the level of the disk’s flux contribution. Both models predict a spectral index k_c for the dust ring with a value of − 0.35 ± 0.2. To better illustrate the visibility plateau and the expected modulation due to the modeled thin ring, we show in the bottom plots of Fig. B.1 three visibility curves corresponding to the best model for three selected baseline orientations and for a fixed wavelength value of 2.15 μm. Regarding the geometrical shape of the circumstellar environment, both models lead to a ring inclined from face-on by 58.5°, with a northeast position angle PA ~ 0°. The inner hole radius is estimated to be r_i ≈ 7.4 mas (0.8 au) from both models, while the ring width isestimated to be ~0.24−0.35 mas (0.03−0.04 au, for the Gaussian-convolved ring model and the geometrically thin ring model, respectively).

Importantly, Fig. B.2 shows the ${\chi }_r^2$ curves of each parameter for the geometrically thin ring model, which helps us to evaluate how well each parameter is constrained by our data. The star near-IR flux contribution is very well constrained as we expected from the plateau seen in the squared visibility curve. The ring spectral index is more loosely constrained, since the best value is consistent with values ranging between − 2 and 2. The ring inclination is constrained between ~ 45° and ~ 75°, while the position angle is less well constrained with two possible minima at ~ 0° and 120°, the former considered to be the absolute minimum. The inner hole radius is constrained to be inside the first 2 au of the system, with a global minimum found around 0.8 au (7.4 mas) and a second (almost equally possible) solution, at ~1.7 au (15.4 mas). Taking into account the upper limit of 10 mas for the radius of the K-band emission found by Monnier et al. (2005), we decided to adopt r_i = 0.8 au. According to Fig. B.2, the ring width w tends toward small values, not larger than ~ 0.3 au. This is discussed further in Sect. 6.1.

5.2 Dust properties through radiative transfer modeling

To strengthen the obtained results and to assess the scenario of an inner ring as close as 0.8 au, gaining further information on the emission properties, we investigated through RT modeling the impact of such a component on the SED. We note here that our aim was not to perform a detailed mineralogy study of the system, but rather to understand how the detected inner dust is consistent with both the near-IR flux and overall SED of the system.

We modeled the multiple and complex outer rings with only three rings based on the results from Thi et al. (2014). In their model, the lower limit particle size was set at 0.5 μm and the upper one at 0.5-1 cm for the two outermost rings and the innermost one, respectively. The grain size follows a distribution ∝ a^−3.5, the surface density profile is a modified version from Li & Lunine (2003), and the flaring index is γ = 1. The three rings peak at ~15, 185, and 300au, with the first two rings separated by a 75 au gap.

Our initial disk setup consists of a disk structure similar to Thi et al. (2014), but with updated stellar parameters (see Sect. 1) and a grain population based on DIANA standard dust grains (Woitke et al. 2016) containing 75% amorphous silicates (e.g., Mg_0.7Fe_0.3SiO₃), 25% porosity, and no amorphous carbon. Our modified grain size distribution and surface density is described in next paragraph. The computed SEDs account for both thermal emission and scattered light contributions.

A silicate dust ring

First, we attempted to reproduce the near-IR excess detected by GRAVITY by including dust grains close to the star, but also taking into account the fact that HD 141569 does not show a silicate emission feature at 10 μm (Seok & Li 2017), which is in part connected to the grain size distribution. We tested several models of the inner ring with the following properties: a varying lower-limit dust grain size of 0.6, 1.2, 2.5, 5.0, 10, 20, 40, 80, 158, and 316 μm; an equal upper-limit grain size of 1 cm, with a size distribution ∝a^−3.5; a surface density ∝r⁻¹; and an inner ring radius fixed at 0.8 au with a width of 0.04 au according to our best-fit model.

All models with grains smaller than ~ 20 μm in the inner ~1 au region can be tuned to reproduce the ~6% near-IR excess, but at the same time they still exhibit a clear 10 μm silicate feature, which is not in adequacy with the observations.

Models accounting only for grains larger than 40 μm result in an almost complete quenching of the 10 μm silicate feature. When testing the mass at ~1 au required not to exceed the mid-IR flux, we find that 10⁻¹⁰ M_⊙ (or 3.3 × 10⁻⁵ M_⊕) would be compliant with this condition, but only a ~1% near-IR excess is generated. On the other hand, the dust mass required to reach a ~ 6% near-IR excess is well beyond 8.5 × 10⁻¹⁰ M_⊙ (or 2.8 × 10⁻⁴ M_⊕), but this then produces a too large mid-IR emission inconsistent with the known SED.

Finally, decreasing the percentage of silicates and increasing the carbon percentage in the dust grains up to 25% did not improve the SED fit. Obtaining a more precise fit to the global SED would require further analysis and tuning of the outer ring contribution in the mid-IR, which is beyond the immediate goal of the paper. However, our modeling seems to point out that it is difficult to reconcile the level of near- and mid-IR excess reported with a model of solely silicate dust in the disk ring at ~ 1 au. Three representative cases of our modeling are presented in Fig. E.1.

A ring of quantum heated particles

Another way to produce near-IR emission consistent with the absence of the prominent 10 μm silicate feature and with the presence of the mid-IR PAH bands is to consider quantum heated particles (QHPs, Purcell 1976; Draine & Li 2001). In the context of interferometric observations, this scenario was invoked for HD 100453 where QHPs were detected in the disk gap (Klarmann et al. 2017), and for HD 179218 with the presence of hot QHPs inside the disk cavity (Kluska et al. 2018).

We tested different amounts of QHPs (from 10⁻¹³ to 10⁻¹⁵ M_⊙), and two different QHP particles sizes (10⁵ and 5 × 10⁵ carbon atoms), based on the results obtained by Klarmann et al. (2017).

The highest masses (~ 10⁻¹³ M_⊙) produce an inner rim emission which results in a near-IR excess that is too large and inconsistent with the GRAVITY measurement. The smallest value (~ 10⁻¹⁵ M_⊙) corresponds instead to an optically thin disk with negligible excess at 2 μm. For a ring geometry in agreement with our best fit of the GRAVITY data, a sweet spot is found for a mass of 4.3 × 10⁻¹⁴ M_⊙ (or 1.4 × 10⁻⁸ M_⊕) for which the resulting disk produces a near-IR excess of ~ 7 %. Smaller particle sizes (e.g., 10² carbon atoms) would require a larger mass reservoir of QHPs to reach a near-IR excess of ~ 6%. As a consequence, a higher mass would result in stronger mid-IR PAH features overestimating the HD 141569 SED observed in the Spitzer IRS spectrum (Sloan et al. 2005). The order of magnitude of ~ 10⁵ carbon atoms per particle appears consistent with the SED profile estimated by GRAVITY and IRS/Spitzer.

The resulting disk SED shows a spectral index (d logF_λ/d logλ) of −1.4 at 2.15 μm which is consistent with the minimization curve of the spectral index in Fig. B.2. The details of the QHP model parameters is shown in Table 2. In Fig. F.1 we show the density and temperature structure of the full-disk model. For the silicate dust in the outer second (~ 15 au), third (~ 185 au), and fourth ring (~ 300 au), the equilibrium temperature varies between 125 K at ~5 au to 20 K in the outermost disk regions. The QHPs used to model the innermost optically thin ring are not in thermal equilibrium. Their temperature distribution depends on their size (the smaller, the hotter) and on the strength of the local ultraviolet (UV) radiation field as they absorb UV photons to quickly re-emit in the near-IR, changing their temperature drastically and very fast. The dark red color in Fig. F.1 shows the location of the QHPs, and according to our RT simulations their temperatures are in the range between 1865 and 95 K.

We find that the strongest constraint on the ring width is set by our interferometric measurement. When considering widths from 0.04 to 0.3 au in our RT modeling, the observed tendency remains the same: the population of silicate dust grains does not satisfactorily reproduce the near- and mid-IR excess, contrary to a population of a stochastically heated small grains.

Beside the ring’s width, the important result of this analysis suggests that a tenuous, QHP-dominated, optically thin inner ring may provide a suitable description of the close circumstellar environment of HD 141569 in agreement with existing detailed modeling of the outer disks. On the contrary, a silicate-dominated dusty inner ring, heavier by about four orders of magnitude and composed of large grains, fails to provide a satisfying description of the system, in particular in terms of flux contribution in the mid-IR spectral range. Figure 3 presents the final result of our MCMax modeling with the parameters of Table 2. We note that the strong PAH feature at 7.8 μm in HD 141569 is not seen in this model because the corresponding opacity has not been added to our models of the outer rings, unlike the models of Thi et al. (2014).

Fig. 3 MCMax SED (black line) of the model described in Table 2. The yellow circles are HD 141569 photometric data listed in Table G.1. The blue line represents the star SED modeled as a black body. The cyan line represents the QHPs SED and the red line the silicates SED. The magenta line is the total disk SED accounting for both QHP and silicate emission. The top right plot focuses on the first 10 μm wavelength range of the SED.

5.3 Spatial scale of the Brγ-line emitting region

We exploited the high-spectral resolution data of GRAVITY in the Brγ region to constrain the spatial scale of the hot gas emitting component, following the formalism in Sect. 4.2. Conservatively, the resulting pure-line visibilities plotted in red in Fig. 2 are very close to 1 for all the baselines. Considering the error bars we can say that the gas region is at the limit of spatially unresolved emission. We therefore propose to constrain the size of the gas emitting region by considering the gas emitting at the Brγ line wavelength peaks (i.e., 2.1654 and 2.1672 μm) using an infinitesimally thin ring model and estimating the upper-limit size that would exceed the error bar of ~ 2% on the pure-line visibilities. In this way, we estimated a maximum radius of ~ 0.35 mas (0.0385 au) for the gas emitting region.

While the analysis of the visibility amplitudes only provides us an estimate of the size scale of the gas emitting region, further information on the spatial and kinematic properties of the gas component is found in the differential phase signal. Typically, differential phases provide information on photocenter displacements along the baselines on angular scales that can surpass the nominal resolution of the interferometer. Figure 4 shows the GRAVITY pure-line differential phases. After removal of the continuum contribution (cf. Sect. 4.2), the S-shape becomes clearly visible for the baselines J3-D0 (100 m, 220°) and J3-K0 (54 m, 151°) around the Brγ line, while a weaker trend is seen for the baselines J3-G2 (62 m, 223°) and D0-G2 (38 m, 36°). The strongest signature shows an amplitude in the differential phase exceeding ~ 20° for J3-K0. Since the differential phase signals for the baselines D0-K0 (95 m, 72°) and G2-K0 (68 m, 92°) are consistent with zero at all wavelengths, we fixed the pure-line differential phases to 0°. The typical uncertainties after correcting for the continuum subtraction are ~ 4°. The resulting deprojected photocenter shifts per spectral channel, with reference frame fixed to the star location, are shown in Fig. 5. We clearly observe that all points are aligned along the same direction with an angle of −10° ± 7°. The redshifted points are located along northwest, while the blueshifted points are located toward the southeast. Based on the redshifted maximum extent of the photocenter shifts, we estimated the radius of the gas region from the differential phases to be 0.333 ± 0.039 mas or 0.037 ± 0.004 au. This appears consistent with the less precise upper-limit size set through the analysis of the visibility amplitudes. We recall however that the size estimate derived through the photocenter shift does not correspond to the physical outer radius of the gas region, but to the size where the gas emission is more intense for a given wavelength.

The distribution of the 2D photocenter solution can be interpreted to the first order as being caused by a gas disk in Keplerian rotation orbiting HD 141569. Under this assumption the analysis of the Brγ emission line’s shape provides further clues to the gas kinematics. We can indeed derive an estimate of the gaseous disk’s radius from the separation between the two peaks of the line and the rest position. From Beckwith & Sargent (1993) we use $$R_{\text{g}}=\frac{G{M}_{\star }}{v_{\text{obs}}^2}{\sin }^{2}i,$$(17)

where R_g is the radius, G is the gravitational constant, M_⋆ the mass of the star, v_obs the projected velocity at the line peaks, and i the disk inclination. The 128 ± 42 km s⁻¹ average peaks shift of the Brγ line with respect to the line rest position leads to an outer limit for the gas region of 0.766 ± 0.554 mas in radius, or 0.084 ± 0.061 au. The resulting error accounts for the uncertainty on the stellar mass (0.01 M_⊙), on the distance (1 pc), on the disk inclination (15°), and on the peak position (3 Å), the last being the dominant one.

The three approaches presented above and based on the analysis of the visibility amplitudes, differential phases, and the spectrum all seem consistent with a gas component in Keplerian rotation confined within ~ 0.8 mas (~ 0.09 au, ~ 12.9R_⋆) in radius. We compared our differential phase signals, our strongest measurable quantity, to a simple geometrical model of an axisymmetric disk in Keplerian rotation built with two thin layers that account for the top and bottom sides of the disk, parameterized by an inner radius r_in (varying from 0.008 to 0.03 au), an outer radius r_out (varying from 0.033 to 0.8 au), and a power-law exponent α (varying from 0 to 4.0) for the disk’s intensity radial profile following I(r) ∝ r^−α. The inclination of the disk was fixed to 58.5°, based on the results of the FT data analysis and under the assumption of coplanar dust-gas rings, and its position angle was fixed to −10°, based on the photocenter shift analysis. The hypothesis of an optically thin disk is made, which implies that the disk intensity does not strongly depend on the disk scale height, which is then fixed to H/R = 0.1 throughout the ring, but only on the surface area of the disk: I ∝dS where dS is an area element (for a complete description of the model, see de Valon et al., in prep.). The model delivers a spectral profile that is eventually normalized and fitted to our experimental double-peaked spectrum in Fig. 6. A grid of 2700 models has been explored, and a best-fit model is found for r_in = 0.011 au, r_out = 0.09 au, and α = 0.5. From this best-fit model, we produced 2D velocity maps (Fig. 7) and retrieved the theoretical pure-line differential phase signal (see Fig. 4). Our Keplerian disk model reproduces qualitatively well the observed pure-line differential phases both in the orientation of the sine wave relative to the wavelength and in amplitude reinforcing the Keplerian gas disk scenario. We also observe that the weakest signals indeed correspond to the baselines D0-K0 and G2-K0 for which the peak amplitudes are ~ 3–5°. On a final note, for the best-fit solution model we tested the optically thick hypothesis by accounting for projection effects on the emissivity law, I ∝dS ⋅e_y = dS ⋅ (n ⋅e_y), where e_y is the line-of-sight unit vector and n the normal unit vector of the area element considered (see de Valon et al., in prep.). The resulting model leads to a spectrum that is similar to that derived from the optically thin model (a small difference is seen at the lowest velocities), and to pure-line differential phases and photocenters shifts consistent with those of the optically thin model.

Fig. 4 HD 141569 spectrum (top), total (black circles and line), and pure-line differential phases (colored signs) along the different baselines. D0-K0 and G2-K0 baseline pure-line differential phases are set to zero. The colors refer to the different spectral channels. Dashed magenta lines represent the pure-line differential phases of the analytical Keplerian disk model described in Sect. 5.3.

Fig. 5 Deprojected photocenter shifts. The colors refer to the different spectral channels and velocities, as shown in Fig. 4. The dashed black line, derived through a linear fit of the photocenter shifts, represents the gas region position angle.

Fig. 6 HD 141569 GRAVITY spectrum (black line), and the spectra of the Keplerian disk models described in Sect. 5.3. The blue line represents the model in the optically thin scenario with residuals given by the red dashed line, while the orange line represents the model in the optically thick scenario with residuals given by the green dashed line.

6 Discussion

6.1 A newly detected inner dust ring

Thanks to GRAVITY, an additional inner ring component at ~ 1 au from the star, narrower than ~ 0.3 au, has been discovered, which adds to the multi-ring picture identified for HD 141569. A schematic visualization of the dust and gas distribution obtained from multi-epoch and multi-instrument observations is depicted in Fig. 8. The low level of near-IR excess (~6%) is comparable to the typical 5% accuracy of K-band photometric data (e.g., 2 MASS), which led to considering the near-IR flux of HD 141569 as essentially photospheric in earlier studies. With this result, the presence of dust as close as 1 au from the star, but well beyond the sublimation radius R_sub ~0.25 au is a more robust piece of evidence. The disk inclination and position angle estimated through the GRAVITY FT data analysis are consistent within the error bars with the values found for the outer rings that make up the circumstellar environment of the system (see Table H.1), suggesting an almost coplanar system of circumstellar rings. Considering how many constraints can effectively be set on the position angle (see Fig. B.2), it does not appear that any clear misalignment between the inner and outer disks could be claimed.

The inclination and position angle could be the reason why PTI and the Keck interferometer were not able to spatially resolve the near-IR emission of HD 141569 in earlier measurements, since the alignment of their single baseline was around 42° from north to east, hence not far from the semi-minor axis of the disk. However, the upper limit of 10 mas in radius for the location of the dust proposed by Monnier et al. (2005) is in agreement with our results. They also report a fractional excess of ~ 5%, but compatible within the error bars with pure photospheric flux in the K band.

Our results are substantially different from those obtained by Lazareff et al. (2017) with PIONIER in the H band. The reported flux excess is ~ 50% larger than in our case and their half-flux radius for the dust emission is only 0.03 au, well inside the sublimation radius of the system expected to be at ~0.21-0.26 au (for L_⋆ = 16.6 − 27.0 L_⊙, T_subl ~ 1470 K and a cooling efficiency ϵ ~0.5). Lazareff et al. (2017) modeled the H-band excess emission using an ellipsoid distribution and not a ring (see their Table B.2 and B.3). The choice of the model and the low-quality PIONIER data could potentially explain this unexpected result.

The values derived for the radius and width of this new innermost dusty component allows us to compare the system with other YSOs. The analysis made by GRAVITY Collaboration (2019) for a sample of 27 Herbig Ae/Be stars revealed dusty circumstellar environments with half-flux radii that range between 0.1 and 6 audepending on the stellar luminosity with a median of 0.6 au, and a width-to-radius ratio w ranging from 0.1 to 1 with a median of 0.83, which is interpreted as smooth and wide rings, even though with large error bars.

Our radius estimate for HD 141569 (~0.8 au) is within the range found by these authors, but its peculiarity is also reflected in the fact that its position in the size-luminosity diagram does not coincide with the bulk of the Herbig stars, reinforcing the idea that HD 141569 is a unique system in terms of evolution. From Table B.2, the width-to-size ratio is estimated to 0.05 ± 0.05 for the best-fit model, but it is also noticeable from Fig. B.2 that the ring’s width is difficult to clearly constrain in the region below 0.3 au. This implies that the width-to-size ratio could be seven times larger. Therefore, our estimate of the width-to-size ratio goes on the lower end of the range found in GRAVITY Collaboration (2019), but remains comparable to systems such as HD 114981 (0.10 ± 0.03) or HD 190073 (0.14 ± 0.03).

In addition to the disk’s well-known flux-size degeneracy, which in our case is broken thanks to the constant plateau in the squared visibilities curve as a function of spatial frequencies, Lazareff et al. (2017) found a negative correlation between the ring width and its half-flux radius, meaning that it gets more difficult to detect with the VLTI baseline small-radius ring-like structures than small-radius ellipsoid structures for angular sizes of the K-band emission of ~ 1 mas or smaller. In our case the size is well constrained being beyond the suggested limit, which is one argument favoring the robustness of our modeling as opposed to the model of a Gaussian brightness distribution that did not converge to a solution.

To further test our model findings, first we tried to add the contribution of a fully resolved emission (also known as a halo) to our geometrically thin ring model following Eq. (4) of Lazareff et al. (2017). Our best-fit solution (${\chi }_{\text{red}}^2=4.69\pm 0.32$) led to results similar to our geometrically thin ring model with a halo contribution that actually converges toward zero (${0.103}_{-0.075}^{+0.144}\%$). Interestingly, this is in full agreement with Lazareff et al. (2017) who found null halo flux contribution as well. Second, we also tried a Lorentzian-convolved infinitesimally thin ring (see Table 5 of Lazareff et al. 2017), which leads to practically the same solution (${\chi }_r^2=4.68\pm 0.32$). We conclude from this detailed analysis that the narrow ring-like shape with a width smaller than 0.3 au is a good description of our observations.

Fig. 7 Monochromatic images (from 2.16452 to 2.16852 μm, i.e., from 221 to 332 km s−1) of the Keplerian ring model described in Sect. 5.3. The colored dashed lines refer to the different GRAVITY baselines. Circles represent the 2D photocenter shift for each baseline.

Fig. 8 Visualization of dust and gas distribution in HD 141569 (adapted from Di Folco et al. 2020). Shown in light orange are the optical–IR dust rings detected in scattered light by HTS (e.g., Augereau et al. 1999; Clampin et al. 2003; Konishi et al. 2016). Overplotted in green are the three near-IR dust ringlets detected by VLT/SPHERE (Perrot et al. 2016). The extended millimeter continuum emission detected by ALMA (Miley et al. 2018) and NOEMA (Di Folco et al. 2020) is shown in orange. The large brown ring represents the mid-IR continuum emission detected by VLT-VISIR and modeled by Thi et al. (2014). In dark red is shown the near-IR dust emission detected by GRAVITY and studied in this work. In the same color is the Brγ line emitting gas region detected by GRAVITY, also analyzed in this work. In blue is depicted the Hα line region based on the upper-limit size estimated by Mendigutía et al. (2017). Finally, in purple is shown the CO gas region whose emissions were detected by ALMA (e.g., White et al. 2016; Miley et al. 2018) and NOEMA (Di Folco et al. 2020).

6.2 Nature and origin of the detected ring

We advance in Sect. 5.2 the scenario of a population of stochastically heated particles (e.g., PAH-like very small grains) as the cause for the near-IR excess. A number of arguments can be discussed in this context.

Our best-fit chromatic model shows a spectral index k_c = −0.35 ± 0.21 for the circumstellar emission. Following Lazareff et al. (2017) and GRAVITY Collaboration (2019), we estimate from the parameter k_c a temperatureof the radiating dust under the gray body hypothesis (i.e., wavelength-independent emissivity) and find T_c = 1460 ± 70 K. In the case of a silicate dust ring in thermal equilibrium at ~ 0.8 au, we would expect a cooler temperature³ of T_c ~650–850 K using the stellar parameter of Table 2 and a cooling efficiency ϵ between 0.3 and 1. Therefore, we argue that the near-IR emission is not dominated by emission of dust in thermal equilibrium,which can be explained by the presence of these small particles that are quantum heated by the stellar UV radiation. Even though the spectral index is found to be not very well constrained (see Fig. B.2), the spectral index corresponding to dust in thermal equilibrium at 800 K would be around +3.4 at λ₀, relatively far from the our best-fit model.

Maaskant et al. (2014) interpret the high intensity ratio between the PAH band at 6.2 μm to the band at 11.3 μm as a tracer of predominantly ionized PAH species located in a disk’s gap and exposed to the intense ionizing UV radiation field of the central star. For instance, with a I_6.2/I_11.3 feature peak ratio of ~3−4 (Seok & Li 2017) the PAH sources IRS 48 and HD 179218 present emission from such predominantly ionized PAHs located in part inside the gap or disk cavity (Maaskant et al. 2014; Klarmann et al. 2017; Kluska et al. 2018; Taha et al. 2018). Interestingly, the I_6.2/I_11.3 peak ratio of HD 141569 derived from Seok & Li (2017) is high as well, estimated to ~ 5−6. This may indicate the presence of PAH species close to the star, with a predominantly ionized state due to the direct irradiation by UV stellar flux, bringing further support to our QHP-dominated inner ring model.

Comparing the IRS spectrum to our model, we find that the our model accounts for 28, 25, and 4% of the observed PAH peak emission for the features at 6, 8, and 11 μm, respectively. The remaining emission would come from PAHs located in the outer rings. Comparing the outer ring PAH mass reservoir estimated by Thi et al. (2014) to that of our innermost ring model, we find that our estimate (4.3 × 10⁻¹⁴ M_⊙) is smaller than their ~15 au and ~ 300 au rings by three orders of magnitude (2.0 × 10⁻¹¹ M_⊙ and 2.1 × 10⁻¹¹ M_⊙, respectively), and by four orders of magnitude for their ~ 185 au ring (1.2 × 10⁻¹⁰ M_⊙) and the entire outer disk environment (1.6 × 10⁻¹⁰ M_⊙). We recall however that in their model Thi et al. (2014) do not account for any dust located at ~ 1 au and suggest a 5 au dust-free inner gap, which could result in overestimated values on their side.

Regarding the origin of the ring structure, the case of HD 141569 is particularly interesting under the aspect of our proposed QHP-dominated inner component: one could question the presence of QHPs in an inner narrow ring since this kind of particle is expected to be coupled to the gas component and to date have been mostly invoked in more extended emission (e.g., Klarmann et al. 2017; Kluska et al. 2018).

Several authors have detected CO emission beyond ~ 10 au (see Introduction), whereas little information is available on the presence of CO within the first 10 au. It is likely however that this inner region is not gas depleted. The [O_I ] λ6300 emission detected by Acke et al. (2005), and suggested by these authors as a dissociation product of OH in the circumstellar disk, would originate between ~0.05 and 0.8 au under the assumption of a gas disk in Keplerian motion (Brittain et al. 2007). Hydrogen recombination and sodium lines were also detected (van der Plas et al. 2015). Both Mendigutía et al. (2017) and our results (see Sect. 6.3) find excited atomic hydrogen in the dust-free cavity, which implies the existence of a replenishment mechanism from the outer regions. Moreover, Brittain et al. (2007) set an upper limit on the column density of CO inside 6 au of N(CO) < 10¹⁵ cm⁻², which translates into a gas mass < 5.9 × 10⁻¹³ M_⊙. Considering our QHPs mass (4.3 × 10⁻¹⁴ M_⊙), a gas-to-dust ratio of 100, and a [CO/H₂] ratio of 10⁻⁴, this would translate in a CO mass of 4.3 × 10⁻¹⁶ M_⊙, which is well below the upper limit and therefore not detected yet. These arguments suggest that a gaseous component may exist and coincide with the proposed QHP component to which it would be coupled.

Furthermore, even though QHPs have been invoked in more extended emission in past works, we cannot exclude the possibility of more compact or narrow components. For example, Khalafinejad et al. (2016) modeled its inner circumstellar region through an optically thin spherical halo extending from 0.1 to 1.7 au in order to explain the near-IR flux of HD 100453 and to fit simultaneously its Q-band flux. The choice of the spherical halo was based on the fact that their data poorly constrained the structure of the inner disk (and so it is the halo extension estimate) and the optically thin hypothesis was set in order to not affect the Q-band flux modeling. This component was also suggested by Klarmann et al. (2017). Their QHP model for HD 100453 underestimates the observed flux in the 1-5 μm wavelength range by up to 30%, and slightly overestimates the long-baseline visibility data, indicating that the missing flux is emitted on short spatial scales. Closer results to those we obtained for HD 141569 were found by Maaskant et al. (2013). These authors suggest a compact optically thin spherical halo for HD 169142 (0.1–0.2 au), HD 135344 B, and Oph IRS 48 (0.1–0.3 au) to reproduce the observed near-IR flux. Several scenarios have been proposed to explain the different structures that protoplanetary disks exhibit, such as gaps, spirals, or rings. Fragmentation of wide rings into narrow ones by secular gravitational instability (e.g., Tominaga et al. 2020), self-induced pileup of particles by aerodynamical feedback (e.g., Gonzalez et al. 2017), and dust traps at local maxima in the gas density due to a reversal of the pressure gradient by dynamical clearing from a companion (e.g., Pinilla et al. 2012) could explain a structured nature of the disk. Our observations leave this matter as an open question since GRAVITY informs us solely on the spatial properties of the detected K-band continuum emission.

6.3 HD 141569 Brγ-line emitting gas region

Our analysis of the kinematic and spatial distribution (via the differential phase) of the hot hydrogen gas is in line with a scenario of a Keplerian disk inside the dust-free cavity. The distribution of the photocenter shifts shown in Fig. 5 agrees well with the behavior expected from a Keplerian disk (Mendigutía et al. 2015). The position angle of the photocenter shifts distribution (−10° ± 7° north to east) is also found to be in overall agreement with the position angle of the inner ring responsible for the near-IR excess, and of the outer rings. Moreover, the photocenters of the Brγ line emitting gas region are located as the photocenters of the outer CO regions, blueshifted ones along the southeast and redshifted ones along the northwest (White et al. 2016). The profiles of pure-line differential phase signals depart a bit from the perfect S-shaped signal expected for a pure Keplerian disk. We believe that it is also limited ultimately by our spectral calibration. In this sense, more accurate measurement of the differential phases in HD 141569 using GRAVITY with the 8 m Unit Telescopes could certainly improve the accuracy of this analysis. In order to evaluate the quality of our spectral data, we chose to compare the GRAVITY profile measured with the ATs to other high-quality spectra obtained with the ISAAC spectrograph at the VLT (Garcia Lopez et al. 2006), the NIRSPEC echelle spectrograph at the Keck Observatory (Brittain et al. 2007), and with SINFONI/VLT from archival data. The comparison is shown in Fig. 9 and we observe that the GRAVITY spectrum exhibits a mild asymmetry between the blue and red peaks. This would suggest that our spectrum could be still affected by some calibration effects, either telluric or instrumental. We then further explored how far our resulting differential phases might be impacted by the slight spectrum asymmetry and tested the derivation of the pure-line differential phases using the SINFONI spectrum, which has a very similar spectral resolution, instead of the GRAVITY spectrum. We found that the 2D distribution of the photocenter shifts remains unchanged within the error bars reported in Fig. 5.

The star is known to be a fast rotator (222.0 ± 7.0 km s⁻¹, Folsom et al. 2012), which results in a small co-rotation radius around 2.38 ± 0.53 R_⊙ (0.011 ± 0.002 au), assuming R_⋆ = 1.5 ± 0.5 R_⊙ (Fairlamb et al. 2015). We cannot exclude that part of the Brγ line emission comes from magnetospheric accretion flows, but the small co-rotation radius compared to the size of the Brγ line emitting region estimated from the SC data analysis (~0.09 au) would not favor this scenario, as opposed to what has been recently found for TW Hya (GRAVITY Collaboration 2020b). Interestingly, the scenario of magnetospheric accretion was also tested by Mendigutía et al. (2017) to explain the Hα double-peaked emission line, but they were not able to reproduce the observed profile with any set of input parameters. Comparing the extent of the Brγ emission to the continuum emission (R_Brγ/R_cont ≈ 0.1), we find that the case of HD 141569 is in contrast with the findings of Kraus et al. (2008). These authors found for a small sample (5 objects) of Herbig Ae/Be stars that those showing a P Cygni Hα line profile and a high mass-accretion rate (> 10⁻⁷ M_⊙ yr⁻¹) seem to show compact Brγ-emitting regions (R_Brγ/R_cont < 0.2), from which the emission stems from magnetospheric accretion or recombination line from ionized hydrogen, while stars showing a double-peaked or single-peaked Hα line profile show a more extended Brγ-emitting region (0.6 ≤  R_{Br γ}/R_cont ≤ 1.4), which would trace a stellar or disk wind. Our system shows mixed features, a Brγ and Hα double-peaked emission line that originates from compact disks in Keplerian rotation where magnetospheric accretion is not the most likely main emission mechanism. Therefore, recombination line emission from ionized hydrogen in an inner gaseous accretion disk as hinted by GRAVITY is a more supported scenario. If that is the case, considering the age of the system and the reported accretion rates between 10⁻⁷ and 10⁻¹¹ M_⊙ yr⁻¹ (Merín et al. 2004; Garcia Lopez et al. 2006; Mendigutía et al. 2011a; Thi et al. 2014; Fairlamb et al. 2015), the inner gaseous disk requires some sort of replenishment mechanism to explain its presence and for it to survive. Replenishment flows, planet-boosted or not (Mendigutía et al. 2017), connecting the inner and outer disk as already observed in other Herbig stars like HD 142527 (Casassus et al. 2013) could be investigated in the future.

Fig. 9 HD 141569 continuum-normalized spectrum taken at different epochs and with different instruments. In blue is depicted thedata taken on July 12, 2019, by GRAVITY; in red data taken in June 2019 by SINFONI; in green data taken in 2002 by the KECK NIRSPEC (Brittain et al. 2007); and in cyan data taken in 2004 by the VLT/ISAAC (Garcia Lopez et al. 2006).

6.4 Hybrid or debris disk?

As mentioned in the introduction, HD 141569 is the only known pre-main sequence star characterized by a hybrid disk. The main characteristic of a hybrid disk is the weak fractional excess of IR emission (8.4 × 10⁻³ for HD 141569, Sylvester et al. 1996), that stems from the optically thin second-generation grains dust component, coupled with the presence of a significant gaseous component, believed to be primordial. Unfortunately the other known best hybrid disk system candidates, which are 49 Cet, HD 21997 (Moór et al. 2011), HD 131835 (Moór et al. 2015), HD 121617, HD 131488 (Moór et al. 2017), and HD 32297 (Moór et al. 2019), have not been studied in their hot dust content. As HD 141569 does, they show a featurelessSED in the near-IR suggesting systems depleted of material inside the first 5 au. However, as we saw in this work, interferometric observations could reveal as-yet-undetected dust, so this scenario should not be excluded for the other mentioned objects. The potential of long-baseline interferometry to detect small levels of hot circumstellar dust emission in supposedly dust-free systems has been indeed exploited for older main sequence stars (Ertel et al. 2014).

In this context it is interesting to compare our results for HD 141569 with Vega, the most iconic debris disk system with a similar spectral type (A0V), but that is significantly older than our source (400–700 Myr). Near-IR excess from Vega was detected and constrained to ~ 1.29 ± 0.19% by interferometric observations with CHARA/FLUOR in theK band (Absil et al. 2006). These authors suggest by SED modeling that the excess comes from hot small grains starting at ~ 0.2–0.3 au with a total dust mass of 8 × 10⁻⁸ M_⊕. Knowing the age of Vega, it is clear that its circumstellar dust is of second generation, a characteristic of debris disks. In the scenario of a silicate-dominated inner ring for HD 141569, the total mass required to induce a larger near-IR excess is, as expected, significantly larger (>10⁻⁴ M_⊕). Considering the younger age of HD 141569 (~7.2 Myr), a more massive inner disk is compatible with a system at an earlier stage of disk evolution. Since the timescale for disk dissipation is known to be about 5−10 Myr (Wyatt 2008), it is plausible that part of the dust in the inner region of HD 141569 is of first generation and a remnant of the primordial circumstellar environment. This would confirm that HD 141569 is closer to a system in the final stage of the protoplanetary disk phase than a debris disk system.

The fact that gas, both H and CO, is detected in HD 141569 is an factor in favor of a system in the (late) protoplanetary disk stage rather than in the debris disk stage. The system, similarly to other hybrid disks, does not follow the correlation between the CO flux density and the millimeter continuum emission followed by T Tauri, Herbig Ae, and debris disks, but lies systematically above the correlation line (Péricaud et al. 2017). The authors suggest that the dust and gas evolution are decoupled, with the dust evolving faster than the gas, leading to an unusual high gas-to-dust ratio (between 135 and 2370 for HD 141569, Di Folco et al. 2020). Other than the primordial origin scenario, a secondary origin of CO was proposed (Kral et al. 2019) in which the gas is self-shielded and shielded by accumulated neutral carbon produced through photodissociation of molecular gas released by planetesimals. This model is able to explain the estimated CO masses in all the hybrid disk candidates (Kral et al. 2019; Moór et al. 2019) except for HD 141569, for which it was not tested. Since CO molecule photodissociation occurs at UV wavelengths, we note that C⁰ could shield these molecules, and also QHPs (Woitke et al. 2016), which are known to absorb UV photons and cool down very quickly by emitting photons in the near-IR. The new detected innermost ring, which we propose in this work to be dominated by a small amount of QHPs, could contribute to the shielding process.

7 Summary

We presented the first GRAVITY interferometric observations of HD 141569. Here we summarize the main conclusions of our work:

• The system was resolved by GRAVITY with squared visibilities down to V² ~ 0.8. If before these observations the near-IR flux contribution of a dust disk was considered absent because no feature was seen in the SED of the object, now, thanks to interferometry, the presence of dust in the first au of the system is a more robust piece of evidence and the flux excess is clearly detected and constrained to be ~6.2% of the total flux.

• Data modeling suggests that the dust is located in a thin ring (≲0.3 au in width) at a radius of ~1 au from the star. The ring shares, within the errors of Fig. B.2, the same inclination (~ 58°) and position angle (~ 0°) as the outer rings observed in the past.

• MCMax SED modeling suggests that this innermost ring could be made of a small amount (1.4 × 10⁻⁸ M_⊕) of QHPs. Large silicate grain models, with and without carbon, can reproduce the 6% flux excess at 2.15 μm, but at the same time they show a significant emission in the mid-IR not consistent with the SED of the system.

• The SC data analysis confirms the significant amount of Brγ line emitting gas already observed in the past. The gas region is spectrally resolved, but spatially unresolved.

• The pure-line differential phases constrain the gas to be in a Keplerian disk-like structure, as hinted by the double-peaked line shape, confined within ~0.09 au (~12.9R_⋆) and oriented in the same was as the outer rings (PA_NE ~−10°).

These results confirm the complexity of the HD 141569 circumstellar environment also at milliarcsecond scale, making the system a unique astronomical laboratory to investigate the missing steps of disk evolution and planet formation theories.

Acknowledgements

We thank the anonymous referee for reviewing this work in detail and for providing useful comments that have improved the quality of the paper. This work is based on observations made with ESO Telescopes at the La Silla Paranal Observatoryunder program IDs 0102.C-0408(D), 0103.C-0347(C), and 0103.C-0347(A). 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. received fundings from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 743 029). L.K. acknowledges the helpful discussions that took place during the ISSI meeting “Zooming In On Rocky Planet Formation” (team 482). We thank Michiel Min for allowing us to use his MCMax code and for his help on understanding the code itself. We thank Anne-Marie Lagrange, Mickaël Bonnefoy and Florian Peissker for providing the SINFONI spectrum of HD 141569. We acknowledge the Gemini Observatoryfor the use of the IR spectrum model of the atmospheric transmission above Cerro Pachon. This research has made use of the model atmosphere grid NeMo, provided by the Department of Astronomy of the University of Vienna, Austria (http://www.univie.ac.at/nemo/). NeMo was funded by the Austrian Science Fonds.

Appendix A GRAVITY observations

Fig. A.1 HD 141569 FT data, squared visibilities, closure phases, and U-V plane coverage, from the three different observation epochs. The colors refer to the different GRAVITY spectral channels.

Appendix B MCMC posterior distribution functions and visibility modeling

Fig. B.1 Results of the FT data fit. Top: MCMC marginal posterior distributions of the fitted parameters. The blue lines identify the median of the distributions for the geometrically thin ring model (left) and for the Gaussian-convolved infinitesimally thin ring model (right). Center: Comparison between the model squared visibilities (blue dots) and the observational data (white–black dots). The red dots represent the absolute residuals. Bottom: Vis2 best fit to the geometrically thin ring model (blue continuous line) at λ=2.153 μm for three different position angles identified in the insets. The visibility data are shown for all spectral channels by the red circles.

Fig. B.2 Reduced chi-squared ${\chi }_r^2$ from the geometrically thin ring model squared visibility fit as a function of the different model parameters. Each point represents a squared visibility fit with the respective parameter fixed to that value while the other parameters are set free.

Appendix C Spectrum wavelength calibration, telluric line correction, and photospheric absorption correction

The 2.15-2.18 μm part of the HD 141569 spectrum is shown in Fig. 2 after continuum-normalization, wavelength calibration, telluric line correction (left panel), and photospheric absorption correction (right panel). The continuum-normalization was done in both the science object spectrum and in our inteferometric calibrator HD 137006 spectrum (α = 230.932°, δ = −01.022°, A5II/III SpT) by fitting the slope of each raw spectrum and dividing the raw spectrum by the resulting fit. The first step of the wavelength calibration was done by comparing the observed wavelength positions of the telluric lines in the HD 141569 spectrum with respect to the positions of the telluric lines present in the IR spectrum of the atmospheric transmission above Cerro Pachon, generated using the ATRAN modeling software (Lord 1992) accounting for a 4.3 mm water vapor column and a 1.5 mm airmass column, available on the Gemini Observatory website. The atmospheric transmission spectrum was convolved by a Gaussian with FHWM of 6 Å to have the same resolution as GRAVITY. The same correction, which results in a 5 Å blueshift, was also applied on the SC visibilities and differential phases, and on the HD 137006 spectrum, our calibrator star. The calibrator was also used as a telluric spectroscopic standard. It shows a photospheric absorption feature at 2.16612 μm that was taken out before the telluric correction on the science object through a NextGen spectrum model (Allard et al. 1997; Hauschildt et al. 1999), available on the GAIA archive, with the following parameters: 7600 K effective temperature, 6.0 surface gravity logarithm, and a − 2.0 solar metallicity. Moreover, the spectrum shows an absorption feature at around 2.16712 μm that is not observed in the spectrum of the other calibrators observed that night and that significantly affects the shape and the intensity of the HD 141569 Brγ red peak. For this reason we chose to calculate the atmospheric transmission function of the night by taking the average of three calibrator spectra, HD 137006, HD 149789, and HD 157029, shown in Fig. C1. The HD 141569 spectrum was then corrected for its radial velocity (−6.4 km/s) and its proper motion with respect to the local standard of rest, resulting in a correction of − 8.7 km/s. Finally, the HD 141569 spectrum was corrected for atmospheric absorption through a spectrum model from the Vienna New Model Grid of Stellar Atmospheres (Heiter et al. 2002)⁴. The model accounts for a star effective temperature of 9800 K, a surface gravity logarithm of 4.4 (Fairlamb et al. 2015), solar metallicity, and microturbulence of 2.0 km/s (Folsom et al. 2012). Since the star is a fast rotator having a projected linear velocity of v sin i = 222.0 ± 7.0 km/s (Folsom et al. 2012), we included the rotation broadening effect on the spectrum model using SPECTRUM (Gray & Corbally 1994). We decreased the intensity of the absorption model by 3 % to better fit our spectrum.

Fig. C.1 Atmospheric transmission functions derived from calibrator spectra: HD 137006 (top), HD 149789 (center), and HD 157029 (bottom). The last plot is the average transmission function of the three spectra.

Appendix D Pure-line visibility in presence of photospheric absorption

The total visibility three-component model used to fit the SC visibilities accounts for the contributions from the star, the circumstellar dust, and the line emitting gas. The total visibility as a function of the wavelength is given by $$\begin{array}{l}{V}_{\text{Tot}}(u,v,\lambda )=\\ \frac{\alpha (\lambda )F_{\text{s}}(\lambda )V_{\text{s}}(u,v,\lambda )+F_{\text{c}}(\lambda )V_{\text{c}}(u,v,\lambda )+F_{\text{L}}(\lambda )V_{\text{L}}(u,v,\lambda )}{\alpha (\lambda )F_{\text{s}}(\lambda )+F_{\text{c}}(\lambda )+F_{\text{L}}(\lambda )}.\end{array}$$(D.1)

The parameter α(λ) is the star continuum-normalized photospheric absorption, which implies that α(λ) = 1 outside the Brγ line and α(λ) < 1 inside the line. This parameter is estimated from the Vienna New Model Grid of Stellar Atmospheres fitting the stellar parameters of HD 141569, as described in Section C. F_s (λ), F_c (λ), and F_L (λ) are the wavelength-dependent fluxes of, respectively, the stellar continuum (i.e., outside the line and temperature-dependent), the dust ring continuum, and the Brγ line emitting gas. F_L (λ) varies across the line and vanishes to zero outside the line. V_s(λ), V_c (λ), and V_L (λ) are the intrinsic visibility functions of each of the three components taken individually. From now on, the star is considered unresolved; therefore, V_s is equal to 1, and we drop the explicit parameter dependencies for convenience. Outside the Brγ line (i.e., in the continuum region) the total visibility is given by $$V_{\text{Tot}}^{\text{Cont}}=\frac{F_{\text{s}}+F_{\text{c}}{V}_{\text{c}}}{F_{\text{s}}+F_{\text{c}}}.$$(D.2)

Using Eq. D.2 to replace V_c, we can rewrite Eq. D.1 as $$V_{\text{Tot}}=\frac{\alpha F_{\text{s}}+V_{\text{Tot}}^{\text{Cont}}(F_{\text{s}}+F_{\text{c}})-F_{\text{s}}+F_{\text{L}}{V}_{\text{L}}}{\alpha F_{\text{s}}+F_{\text{c}}+F_{\text{L}}},$$(D.3)

and using the definition of the line-to-continuum flux ratio (Eq. 13), Eq. D.3 becomes $$V_{\text{Tot}}=\frac{F_{\text{s}}(\alpha -1)}{F_{\text{Tot}}}+\frac{V_{\text{Tot}}^{\text{Cont}}}{F_{\text{L/C}}}+\frac{F_{\text{L}}{V}_{\text{L}}}{F_{\text{Tot}}},$$(D.4)

where F_Tot = α F_s + F_c + F_L. We note once again that F_L/C is the raw line-to-continuum ratio including the photospheric absorption. This quantity corresponds to the top left spectrum in Fig. 2. We also note that outside the line (i.e., α=1 and F_L =0), Eq. D.4 simplifies to Eq. D.2. Now, making use of the parameter β of Eq. 14, we can write $$\frac{F_{\text{L}}}{F_{\text{Tot}}}=1-\frac{\alpha +\beta }{1+\beta }\frac{1}{F_{\text{L/C}}},$$(D.5)

and finally Eq. D.4 becomes $$V_{\text{Tot}}=\frac{\alpha -1}{1+\beta }\frac{1}{F_{\text{L/C}}}+\frac{V_{\text{Tot}}^{\text{Cont}}}{F_{\text{L/C}}}+V_{\text{L}}(1-\frac{\alpha +\beta }{1+\beta }\frac{1}{F_{\text{L/C}}}).$$(D.6)

Solving Eq. D.6 for V_L, and noting that from our data $V_{\text{Tot}}=V_{\text{Tot}}^{\text{Cont}}$ since the SC visibilities are spectrally flat for all wavelengths and baselines, we obtain Eq. 12: $$V_{\text{L}}=\frac{V_{\text{Tot}}[(1+\beta )F_{\text{L/C}}-1-\beta ]-\alpha +1}{(1+\beta )F_{\text{L/C}}-\alpha -\beta }.$$(D.7)

This equation tells us that the pure-line visbility can be estimated from the total visibility in the line (which in our case is comparable to the total visibility in the continuum) if the photospheric absorption profile can be estimated, the continuum disk-to-star flux ratio β is known and the continuum-normalized spectrum of the line is accessible. We note, in the case where $V_{\text{Tot}}=V_{\text{Tot}}^{\text{Cont}}$, that the absence of photospheric absorption (i.e., α=1) leads simply to V_L = $V_{\text{Tot}}^{\text{Cont}}$=V_Tot. Finally, we see that the pure-line visibility V_L is lower than 1 only when F_L/C is greater than 1, which is equivalent to detecting the line above the continuum.

Appendix E MCMax silicate model

Here are shown three representative cases of our SED modeling through a silicate dust ring, as described inSection 5.2. The first model accounts for particles with a lower-limit size of 1.26 μm, a dust mass of 10⁻⁵ M_⊕, and 0 % carbon with a resulting temperaturedistribution in the range 570 − 1070 K. The second model accounts for particles with a lower-limit size of 40 μm, a mass of 1.17× 10⁻⁴ M_⊕, and 0 % carbon with a resulting temperaturedistribution in the range 620 − 840 K. The last model accounts for particles with a lower-limit size of 40 μm, a mass of 3.33× 10⁻⁴ M_⊕, and 25 % carbon with a resulting temperaturedistribution in the range 440 − 770 K. We note that in all three cases the near-IR excess is smaller than that obtained through the GRAVITY data analysis, while the mid-IR emission exceeds the photometry data. Decreasing the mid-IR excess would at the same time decrease further the near-IR value, suggesting that a pure silicate innermost ring model is not compatible with the nature of the HD 141569 system.

Fig. E.1 MCMax models that account for the emission of the star and only the first innermost silicate ring located at 0.8 au and 0.04 au in width. The models differ for the particle size lower limit, dust mass, and carbon percentage. The top plot model accounts for a dust mass of 10−5 M⊕ and no carbon; the center plot model accounts for a mass of 1.17 × 10−4 M⊕ and no carbon; the bottom plot model accounts for a mass of 3.33 × 10−4 M⊕ and 25 % carbon.

Appendix F MCMax model density and temperaturestructure for the QHP model

Fig. F.1 Density and temperature structure of the QHP model described in Section 5.2 and Table 2. The x-axis shows the radius in au, the y-axis the height divided by the radius. The colors of the top plots indicate the mass density, the colors of the bottom plot the temperature. QHPs are located at 0.8 au in a very narrow (0.04 au) ring. The optical depth at the observation wavelength is indicated by the radial (dotted line) τ = 1 surface. Themagenta dashed line represents the sublimation radius of the system. The dark red color in the temperature plot refers only to the location of the QHPs and not their temperature, since QHPs do not have an equilibrium temperaturevalue.

Appendix G Photometric data

Appendix H Inclination and position angle of outer disks from the literature

References

All Tables

All Figures

	Fig. 1 HD 141569 FT data; squared visibilities (left panel), closure phases (central panel), and U-V plane coverage (right panel), fromall the observation epochs. Colors refer to the different GRAVITY spectral channels.
In the text

	Fig. 2 Science channel data of July 2019. Top left: wavelength-calibrated continuum-normalized spectrum,corrected for telluric lines. Top right: same as top left, but corrected for the Brγ photospheric absorption. Differential phases (left column) and visibilities (right column) along the six GRAVITY baselines. The red lines in the visibility plots show the pure-line visibilities.
In the text

Fig. 3 MCMax SED (black line) of the model described in Table 2. The yellow circles are HD 141569 photometric data listed in Table G.1. The blue line represents the star SED modeled as a black body. The cyan line represents the QHPs SED and the red line the silicates SED. The magenta line is the total disk SED accounting for both QHP and silicate emission. The top right plot focuses on the first 10 μm wavelength range of the SED.

In the text

	Fig. 4 HD 141569 spectrum (top), total (black circles and line), and pure-line differential phases (colored signs) along the different baselines. D0-K0 and G2-K0 baseline pure-line differential phases are set to zero. The colors refer to the different spectral channels. Dashed magenta lines represent the pure-line differential phases of the analytical Keplerian disk model described in Sect. 5.3.
In the text

	Fig. 5 Deprojected photocenter shifts. The colors refer to the different spectral channels and velocities, as shown in Fig. 4. The dashed black line, derived through a linear fit of the photocenter shifts, represents the gas region position angle.
In the text

	Fig. 6 HD 141569 GRAVITY spectrum (black line), and the spectra of the Keplerian disk models described in Sect. 5.3. The blue line represents the model in the optically thin scenario with residuals given by the red dashed line, while the orange line represents the model in the optically thick scenario with residuals given by the green dashed line.
In the text

	Fig. 7 Monochromatic images (from 2.16452 to 2.16852 μm, i.e., from 221 to 332 km s−1) of the Keplerian ring model described in Sect. 5.3. The colored dashed lines refer to the different GRAVITY baselines. Circles represent the 2D photocenter shift for each baseline.
In the text

Fig. 8 Visualization of dust and gas distribution in HD 141569 (adapted from Di Folco et al. 2020). Shown in light orange are the optical–IR dust rings detected in scattered light by HTS (e.g., Augereau et al. 1999; Clampin et al. 2003; Konishi et al. 2016). Overplotted in green are the three near-IR dust ringlets detected by VLT/SPHERE (Perrot et al. 2016). The extended millimeter continuum emission detected by ALMA (Miley et al. 2018) and NOEMA (Di Folco et al. 2020) is shown in orange. The large brown ring represents the mid-IR continuum emission detected by VLT-VISIR and modeled by Thi et al. (2014). In dark red is shown the near-IR dust emission detected by GRAVITY and studied in this work. In the same color is the Brγ line emitting gas region detected by GRAVITY, also analyzed in this work. In blue is depicted the Hα line region based on the upper-limit size estimated by Mendigutía et al. (2017). Finally, in purple is shown the CO gas region whose emissions were detected by ALMA (e.g., White et al. 2016; Miley et al. 2018) and NOEMA (Di Folco et al. 2020).

In the text

	Fig. 9 HD 141569 continuum-normalized spectrum taken at different epochs and with different instruments. In blue is depicted thedata taken on July 12, 2019, by GRAVITY; in red data taken in June 2019 by SINFONI; in green data taken in 2002 by the KECK NIRSPEC (Brittain et al. 2007); and in cyan data taken in 2004 by the VLT/ISAAC (Garcia Lopez et al. 2006).
In the text

	Fig. A.1 HD 141569 FT data, squared visibilities, closure phases, and U-V plane coverage, from the three different observation epochs. The colors refer to the different GRAVITY spectral channels.
In the text

Fig. B.1 Results of the FT data fit. Top: MCMC marginal posterior distributions of the fitted parameters. The blue lines identify the median of the distributions for the geometrically thin ring model (left) and for the Gaussian-convolved infinitesimally thin ring model (right). Center: Comparison between the model squared visibilities (blue dots) and the observational data (white–black dots). The red dots represent the absolute residuals. Bottom: Vis2 best fit to the geometrically thin ring model (blue continuous line) at λ=2.153 μm for three different position angles identified in the insets. The visibility data are shown for all spectral channels by the red circles.

In the text

	Fig. B.2 Reduced chi-squared ${\chi }_r^2$ from the geometrically thin ring model squared visibility fit as a function of the different model parameters. Each point represents a squared visibility fit with the respective parameter fixed to that value while the other parameters are set free.
In the text

	Fig. C.1 Atmospheric transmission functions derived from calibrator spectra: HD 137006 (top), HD 149789 (center), and HD 157029 (bottom). The last plot is the average transmission function of the three spectra.
In the text

Fig. E.1 MCMax models that account for the emission of the star and only the first innermost silicate ring located at 0.8 au and 0.04 au in width. The models differ for the particle size lower limit, dust mass, and carbon percentage. The top plot model accounts for a dust mass of 10−5 M⊕ and no carbon; the center plot model accounts for a mass of 1.17 × 10−4 M⊕ and no carbon; the bottom plot model accounts for a mass of 3.33 × 10−4 M⊕ and 25 % carbon.

In the text

Fig. F.1 Density and temperature structure of the QHP model described in Section 5.2 and Table 2. The x-axis shows the radius in au, the y-axis the height divided by the radius. The colors of the top plots indicate the mass density, the colors of the bottom plot the temperature. QHPs are located at 0.8 au in a very narrow (0.04 au) ring. The optical depth at the observation wavelength is indicated by the radial (dotted line) τ = 1 surface. Themagenta dashed line represents the sublimation radius of the system. The dark red color in the temperature plot refers only to the location of the QHPs and not their temperature, since QHPs do not have an equilibrium temperaturevalue.

In the text