The GRAVITY young stellar object survey in 51 at sub-au scales

Context. 51Oph is a Herbig Ae/Be star that exhibits strong near-infrared CO ro-vibrational emission at 2.3 µ m, most likely originating in the innermost regions of a circumstellar disc. Aims. We aim to obtain the physical and geometrical properties of the system by spatially resolving the circumstellar environment of the inner gaseous disc. Methods. We used the second-generation Very Large Telescope Interferometer instrument GRAVITY to spatially resolve the continuum and the CO overtone emission. We obtained data over 12 baselines with the auxiliary telescopes and derive visibilities, and the differential and closure phases as a function of wavelength. We used a simple local thermal equilibrium ring model of the CO emission to reproduce the spectrum and CO line displacements. Results. Our interferometric data show that the star is marginally resolved at our spatial resolution, with a radius of ∼ 10 . 58 ± 2 . 65 R (cid:12) . The K -band continuum emission from the disc is inclined by 63 ◦ ± 1 ◦ , with a position angle of 116 ◦ ± 1 ◦ , and 4 ± 0.8mas (0.5 ± 0.1au) across. The visibilities increase within the CO line emission, indicating that the CO is emitted within the dust-sublimation radius. By modelling the CO bandhead spectrum, we derive that the CO is emitted from a hot ( T =1900–2800K) and dense ( N CO = (0.9–9) × 10 21 cm − 2 ) gas. The analysis of the CO line displacement with respect to the continuum allows us to infer that the CO is emitted from a region 0.10 ± 0.02au across, well within the dust-sublimation radius. The inclination and position angle of the CO line emitting region is consistent with that of the dusty disc. Conclusions. Our spatially resolved interferometric observations conﬁrm the CO ro-vibrational emission within the dust-free region of the inner disc. Conventional disc models exclude the presence of CO in the dust-depleted regions of Herbig AeBe stars. Ad hoc models of the innermost disc regions, that can compute the properties of the dust-free inner disc, are therefore required.


Introduction
Circumstellar discs are crucial for understanding how stars and planets form. In spite of the enormous wealth of data provided by a new generation of instruments, such as Atacama Large Millimeter/submillimeter Array (ALMA) and VLT Spectro-Polarimetric High-contrast Exoplanet REsearch (VLT/ SPHERE), understanding disc properties and evolution is still challenging. This is mostly due to the lack of knowledge about the physical properties and structure of the innermost disc regions. Even if these new instruments would allow us to spatially resolve the disc emission in the dust continuum and in a number of molecular lines down to scales of a few au, it is still demanding to resolve the structure of the innermost disc Corresponding author; e-mail: maria.koutoulaki@eso.org regions, and this has only been possible for a few objects (Kraus et al. 2008;Eisner et al. 2007;Lazareff et al. 2017;GRAVITY Collaboration 2019). Still, this is a region of great importance, as we expect that both the accretion of disc matter onto the central star and the ejection of magnetically controlled winds occur in a region smaller than about one au, even for young stars in nearby star-forming regions (120-140 au). Spatially resolving the inner disc requires optical interferometry. The accessible tracers are limited to the continuum emission of dust close to the sublimation radius and to a few emission lines, the most prominent of which is the HI Brγ line. It has been found so far that with very few exceptions, the Brγ line mainly originates in a wind (Weigelt et al. 2011;Caratti o Garatti et al. 2015;Garcia Lopez et al. 2015). In a few objects, however, the overtone ro-vibrational emission of CO at ∼2.3 µm is detected and has A&A 645, A50 (2021) Notes. (a) Auxiliary Telescope, (b) high spectral resolution in K-band, (c) uniform-disc diameter derived from the software package SearchCal from Jean-Marie Mariotti Center (JMMC).
been proved to trace gas in Keplerian rotation very close to the star, very likely in the disc itself (Carr 1989;Kraus et al. 2000;Bik & Thi 2004;Tatulli et al. 2008; GRAVITY Collaboration 2020). By model-fitting the line profile of the CO bandhead emission obtained with high spectral resolution observations, several authors showed that the emission comes from a relatively warm (T ∼ 2000-5000 K) and dense gas (N CO ∼ 10 20 −10 22 cm −2 ) in Keplerian rotation (Carr 1989;Najita et al. 1996;Ilee et al. 2014;Koutoulaki et al. 2019). However, this information has so far not been used to further constrain the physical, chemical, and thermal processes occurring in the inner disc. This is mostly because we lack of a solid knowledge of the stellar properties and disc geometry, which prevents an accurate estimation of the location of the CO emission from the inferred Keplerian velocity. Further progress requires directly measuring the location of the CO emission, using near-IR interferometry.
In this paper, we present a near-IR interferometric study of the star 51 Oph (HD 158643) using the ESO VLTI instrument GRAVITY. 51 Oph is a fast-rotating star (v rot sin i = 267 ± 5 km s −1 ; Dunkin et al. 1997) of about 4 M , located at a distance D = 123 ± 4 pc (Gaia Collaboration 2018). The spectral type is B9.5-A0 (Dunkin et al. 1997;Gray & Corbally 1998) and the luminosity is ∼250 L (van den Ancker et al. 2001). The spectrum is rich in molecular and atomic emission lines in the mid-and far-IR (Thi et al. 2013a). 51 Oph shows bright 2.3 µm CO overtone emission (Thi et al. 2005;Berthoud et al. 2007;Tatulli et al. 2008). Tatulli et al. (2008) spatially resolved the CO overtone emission in this star, using the first generation of near-IR interferometric instrument, VLTI/AMBER. They find that the CO emission originates within 0.15 au of the star and rotates in a Keplerian manner. However, the number of baselines, the resolution, and the signal-to-noise ratio of these early near-IR interferometric studies were limited, and therefore failed to achieve well-constrained physical properties of the disc-emitting region.
The structure of the paper is as follows. The observations and data reduction are discussed in Sect. 2, and the results, modelling, discussion, and conclusions are presented in Sects. 3-6 respectively.

Observations and data reduction
51 Oph was observed using the VLTI with the K-band beam combiner GRAVITY (GRAVITY Collaboration 2017) on two dates, 29 May 2017 and 15 August 2017. The observations were performed using the 1.8 m auxiliary telescopes (ATs) with the configurations B2-K0-D0-J3 and A0-G1-J2-K0 in May and August, respectively. The resulting projected baselines and position angles (PAs) are shown in Table 1. The observations were performed in single-field high-combined mode of GRAVITY, in which the visibilities of the science target are recorded on both the fringe-tracker detector (at low spectral resolution, R = 30) and the science detector (at high spectral resolution, R = 4000). The detector integration time (DIT) on the fringe tracker and science detector was set to 0.85 and 30 ms, respectively, with a total exposure timne of 5 min per block on the science detector. This resulted in a total integration time on source of 2 and 1 h for the first and second night, respectively.
The data were reduced using the instrument pipeline (version 1.0.11). In order to calibrate the transfer function, the calibrator star HD 163955 was used. The wavelength calibration was refined using the atmospheric telluric absorption lines present in the spectra of both target and calibrator. A shift of 7 Å was found and taken into account. In addition, the spectrum of the calibrator was used to correct the target spectrum for the absorption telluric features and the instrumental response. Finally, the 2MASS K-band magnitude of K = 4.3 was adopted to flux calibrate the spectrum of 51 Oph.
For the continuum analysis, the fringe-tracker data were used, as the fringe tracker operates at higher frequencies, making it less sensitive to the atmospheric changes. There is no significant difference between the level of the calibrated visibilities measured by the fringe-tracker and the science detector.

Results
For each epoch, GRAVITY observations of 51 Oph provide us with a spectrum, six spectrally dispersed visibilities and differential phases, and four closure phases (Fig. 3). In both epochs, the K-band spectrum clearly shows the CO overtone emission and bright H I Brγ. Figure 1 shows the flux-calibrated spectrum of 51 Oph for the May dataset. The continuum K-band contribution is also overplotted (dashed red line) in the same figure. This continuum was used to derive the continuum-subtracted spectra in the text (for the May and August spectra) and was computed as the sum of a black body emission at the effective temperature of the star (T = 10 000 K, Gray & Corbally 1998), which contributes 80% to the total K-band flux, plus additional black body emission at a temperature of 1500 K, representing the K-band dust emission, and contributing 20% to the total K-band flux (Lazareff et al. 2017; GRAVITY Collaboration 2019, hereafter Paper I). No sign of variability is observed between epochs. However, the August dataset is about three times noisier than the May dataset (see Fig. 2). We focus this paper on the CO emission, and defer the analysis of Brγ to a following paper.

CO spectrum
The 51 Oph continuum-subtracted spectrum is shown in Fig. 2 for the May (black) and August (red) 2017 datasets. The A50, page 2 of 21 GRAVITY Collaboration: The GRAVITY young stellar object survey Fig. 1. K-band flux-calibrated spectrum of 51 Oph for the May 2017 dataset. The dashed red line shows the adopted continuum distribution, computed as the sum of a black body at T = 10 000 K, which contributes 80% to the total flux, and a 20% contribution from a black body at T = 1500 K, approximately the silicate sublimation temperature (Lazareff et al. 2017; GRAVITY Collaboration 2019). emission of the first four ro-vibrational transitions (υ = 2−0, 3−1, 4−2, and 5−3), and an indication of the fifth (υ = 6−4) is clearly observed.
The profile of the first three bandheads shows a sharp peak at similar intensity, while the last two have slightly lower peak values. The so-called blue shoulder is observed in the first CO bandhead even at the low spectral resolution of GRAVITY. The blue shoulder is very likely due to Keplerian broadening (e.g. Chandler et al. 1995), which causes a broadening towards shorter wavelengths of the otherwise intrinsic sharp flux increase of the bandhead. No individual or blended CO J components are observed in the band tails because the spectral resolution and signal-to-noise-ratio of our spectra are not high enough, and also because the many telluric atmospheric absorption lines in this region of the spectrum were imperfectly removed. As previously mentioned, the spectra at the two epochs are very similar, with marginal differences in the third and fourth bandheads that are well within the uncertainties.

Interferometric observables
Visibilities, differential and closure phases were measured in 51 Oph at 12 different baselines (6 for each epoch, see Fig. 3).
The visibilities in all the baselines show that the continuum is marginally resolved with values ranging from 0.7 to 0.9. In several baselines, the visibility at the peak of the bandheads is slightly higher than that of the continuum, suggesting a smaller emitting region.
A clear differential phase signal is detected in the first four CO bandheads for most baselines in both epochs, indicating a clear shift of the photo-centre of the CO emission with respect to that of the continuum.
No closure phase is detected within the calibration error (∼5 • ), indicating a symmetric circumstellar environment.

Circumstellar environment: the continuum emission
The observed continuum emission is the sum of the emission of the central star and that of surrounding circumstellar matter, which very likely originates from a disc. The relative flux and size of the two components are different, and the resulting observed visibility V c can be written as where F * and V * are the flux and visibility of the star, F disc and V disc are the flux and visibility of the disc, and F c = F * + F disc is A50, page 3 of 21 A&A 645, A50 (2021)  the total continuum flux. In 51 Oph, the dominant flux contribution in the K band is that of the star, which provides 80% of the total flux (see Sect. 3). The total continuum visibility (star plus disc) is fitted with a two-Gaussian model: a one-dimensional Gaussian for the star, with the stellar angular radius as the single free parameter; and a two-dimensional Gaussian for the disc, characterised by three free parameters: radius, inclination, and position angle of the major axis of the disc.
Our best-fitting model gives a stellar radius (half width at half maximum) of R * = 0.4 ± 0.1 mas and a disc radius of R disc = 4.0 ± 0.8 mas with an inclination (i) and major-axis PA of i = 52 • ± 1 • and PA = 116 • ± 1 • , respectively (see Table 2 and Fig. 4, solid black line). At the distance of 51 Oph, this translates into a stellar radius of R * = 0.05 ± 0.01 au, or R * = 10.6 ± 2.6 R , and a disc emission about ten times higher (i.e. R disc = 0.49 ± 0.10 au).
Our measurement of the stellar radius is in agreement with the measurements obtained at 6000 Å with the CHARA interferometer by Jamialahmadi et al. (2015), although the CHARA observations, thanks to the much longer baselines, reveal a highly flattened stellar photosphere of radius 0.42 ± 0.01 × 0.6 ± 0.05 mas. It should be mentioned that, as discussed in Jamialahmadi et al. (2015), the stellar luminosity computed from the interferometric size and the effective temperature derived from the spectral type of 51 Oph would be about 10 3 L , four times higher than the values derived from the observed SED (van den Ancker et al. 2001). This suggests that 51 Oph is a giant A50, page 4 of 21 GRAVITY Collaboration: The GRAVITY young stellar object survey star rotating close to break-up, either in the pre-main sequence and still in the contraction phase, or a more evolved star already in the post-main sequence evolutionary phase (see Jamialahmadi et al. 2015 for more details).
Our derived disc radius, inclination, and PA are also roughly in agreement with those found in Paper I from the same set of data using a variety of different geometrical models to account for the disc emission. The main difference between our values and those found in Paper I comes from the fact that in the latter, the star was assumed to be a point source (i.e. V * = 1).
As discussed in detail in Paper I, the R disc value is consistent within the uncertainties with the expected location of the silicate sublimation region. When a stellar luminosity of 250 L and a silicate sublimation temperature of 1500 K are assumed, the dust sublimation radius is located between 0.5 to 1.7 au, depending on whether small (cooling efficiency = 0.1) or large (cooling efficiency of ∼1) silicate grains are considered (see Paper I, Kama et al. 2009;Isella & Natta 2005;Dullemond et al. 2001,

Circumstellar environment: the CO emitting region
To compute the size of the CO line emitting region, we subtracted the contribution of the star and the disc and defined the continuum-subtracted visibility at any given wavelength within the CO bandheads V line as In this expression, F line is the line flux, F tot = F c + F line is the observed flux, V tot is the observed visibility, V c and F c are the continuum visibility and the continuum flux, and Φ is the observed differential phase signal (see Weigelt et al. 2011, for more details). The continuum visibility is measured in the nearest line-free wavelength interval, and the continuum flux is taken at each wavelength from the continuum fit shown in Fig. 1 (dashed red line). V line and F line depend on the observed continuum properties (visibility and line-to-continuum ratio), not on the continuum deconvolution in stellar + disc components. In our observations, the differential phase is very small (<5 • ). Thus Eq.
(2) can be approximated as The errors are calculated by propagating Eq.
(3). The continuum-subtracted visibilities at the peak of the first and second bandheads are computed by averaging over three consecutive frequencies at each bandhead. The results are shown in Fig. 4. The continuum-subtracted visibilities for the first and second bandheads are very similar, and in all cases, they are slightly higher than the continuum visibilities. A 2D elliptical Gaussian fitting to the continuum-subtracted CO visibilities (assuming the inclination and PA of the disc) gives a radius of the CO emitting region of R CO = 0.7 mas (i.e. 0.09 au at a distance of 123 pc). This value is much lower than the continuum disc radius. However, given that the CO is barely resolved, it should be considered as a rough estimate of the real size of the CO emitting region.
A much more accurate estimate of the geometrical properties of the CO emitting region can be derived from the differential phase data (Fig. 3, middle panel). As for the visibility, the continuum contribution can be removed from the observed values. The continuum-subtracted differential phase (∆φ) at any wavelength within the line is computed following the expression (Weigelt et al. 2007) with the same symbols as in Eq.
(2). The displacement of the photo-centre of the emission at any given wavelength (δ) is then A50, page 5 of 21 A&A 645, A50 (2021) computed as: where B is the length of the baseline.
The δ values computed for the first three bandheads as a function of the wavelength are shown in Figs. 5-7 for the 12 observed projected baselines (PBLs). Under the assumption that the continuum emission is centrally symmetric, the quantity δ provides a measurement of the photocentre displacement of the line-emitting region at the specific wavelength within the line with respect to the continuum. In the simplest case of an individual line, such as HI Brγ , emitted by a tilted rotating ring, δ will vary from a maximum value in the blue wing of the line to a minimum in the centre and again to a maximum of opposite sign in the red wing. This typical behaviour has been used to measure the size, inclination, and PA of the Brγ emitting region, as well as the rotation direction, in the Herbig Ae HD 163296 (e.g., Ellerbroek et al. 2015). The case of the CO bandheads is more complex, as at any given wavelength the emission is due to the overlap of various blue-and red-shifted individual J components. As a consequence, the simple association of wavelength and location of the emitting region is lost (see e.g. Koutoulaki et al. 2019, for more details) and a proper modelling of the bandhead emission is required to interpret the observed displacement of the CO emission with wavelength (see Sect. 4.2).

Modelling the CO spectrum
Following Koutoulaki et al. (2019), and Kraus et al. (2000), we computed synthetic spectra and intensity maps of the CO bandhead emission as a function of wavelength.
In brief, the CO emitting region was modelled as a narrow ring of gas with temperature (T), and column density (N CO ), rotating at a Keplerian velocity v rot , around the central star. In each vibrational band, we considered 100 rotational transitions; each J component has a Gaussian profile of width ∆v (the socalled intrinsic line width). We created a spectrum over the range covered by the observed bandheads, with very high spectral resolution (R = 100 000), so that each individual J transition profile is well sampled. We first computed the total optical depth at each wavelength, then the intensity given by I(ν, J) = B ν (T )(1 − e −τ ν,J ). The resulting spectrum was then convolved over the Keplerian velocity pattern, reduced to the instrumental spectral resolution of R = 4000, and normalised to the peak of the first bandhead. The bandhead profiles depend on four free parameters: the A50, page 6 of 21 GRAVITY Collaboration: The GRAVITY young stellar object survey temperature (T), the column density of CO (N CO ), the intrinsic line width ∆v, and the Keplerian velocity of the ring projected along the line of sight (v rot sin i).
Even at the GRAVITY spectral resolution of R = 4000, the shape of the bandheads provide some interesting constraints to the physical properties and location of the CO emission. To illustrate this, we computed spectra by varying the four free parameters over a wide range of values, and we estimated the quality of the fit by considering how well different specific spectral features are reproduced. A series of examples on how the CO bandhead spectrum changes as function of the different free parameters, and the corresponding comparison with the observed spectrum is shown and discussed in Appendix A. Our main findings are described below.
The first three bandheads observed in 51 Oph are due to the overlap of very optically thick individual J components (see Fig. A.1). As a consequence, the two parameters T and N CO have large uncertainties because an increase in T can be compensated for by a decrease in N CO , for instance (see Figs. A.5 and A.6). The fourth and fifth bandheads are optically thinner than the first three bandheads and therefore can give better constraints on T and N CO (see Fig. A.1). However, the line profiles of these two bandheads are more affected by the poor atmospheric transmission at the edge of the K band which makes their line profiles noisier and thus less reliable. We find that the observed spectrum is well reproduced by a range of models with pairs of temperature and column density values ranging from T = 1900-2800 K and N CO = 9 × 10 20 -9 × 10 21 (see Fig. A.4). The range of acceptable pairs of T and N CO was selected by simultaneously increasing (decreasing) the T and decreasing (increasing) the N CO values until the observed spectrum (peak and/or tail of the bandheads) differed by more than two σ from the synthetic spectrum. Outside the range of acceptable values (e.g. lower temperature values coupled with higher column densities or vice versa) no reasonable fit is found as the bandheads become less optically thick, and therefore they are more sensible to variations in T and N CO (see Figs. A.2, and A.3).
The observed spectrum is well reproduced by v rot sin i values in the range of 130 ± 30 km s −1 . Different values of v rot sin i change the width of the blue shoulder (especially for the first bandhead), as well as the location of the bandhead peaks (see Fig. A.9). The range of acceptable values of v rot sin i has been defined as that for which the synthetic and observed spectrum differ by less than two spectral channels.
The intrinsic line width ∆v of the single bandhead components is a very interesting parameter related to the turbulence level. However, this parameter is very difficult to measure even at high spectral resolution, as its value is typically much lower than the line Keplerian broadening (Lee et al. 2016). In the case of 51 Oph, we find acceptable values of ∆v in the range of 15 ± 5 km s −1 (see Fig. A.8). As previously discussed, the individual J components in the CO spectrum of 51 Oph are very optically thick in the bandhead tails as well. This means that increasing ∆v increases the emission in the line wings in such a way that, at the GRAVITY resolution, the bandhead tails look like a continuum of increasing value with respect to the peak. As done previously, the acceptable range of parameters was selected by modifying ∆v until the observed spectrum (peak and/or tail of the bandheads) differed by more than two σ from the synthetic spectrum. The best range of acceptable fitting parameters is shown in Fig. 8 and reported in Table 3.

Modelling the CO line displacement
To derive the synthetic values of the CO line displacement, δ, the azymuthal distribution of the CO intensity along the modelled ring was computed as a function of wavelength for the first three bandheads. Examples of intensity maps computed for an inclination of 70 • , a size of the CO emitting region of 0.1 au, and the physical properties of T = 2400 K, N CO = 4 × 10 21 cm −2 , ∆v = 15 km s −1 , and v rot sin i = 130 km s −1 are shown in Fig. 9 for the first three CO bandheads at a spectral resolution of R = 4000 and three velocities. As pointed out in Sect. 3.4, only the blue-shifted emission is confined to a specific region of the ring (as would be the case for an individual line), while the emission is much more homogeneous at all other wavelengths. It should be mentioned that as long as T and N CO vary in a correlated way, the intensity distribution does not vary significantly within the range of acceptable values reported in Table 3.
By inclining and rotating these intensity maps along the line of sight, synthetic δ values can be computed and used to further constrain the geometry and size of the CO emitting region. This was done by computing as a function of wavelength the photo-centre shifts with respect to the centre of the ring for the 12 observed PAs. In principle, the variation of δ with wavelength depends on the same physical parameters that reproduce the CO spectrum (T , N CO , and ∆v), and on the geometry (i.e. i and PA), kinematics (v rot ) and size of the emitting region. We then started by creating an intensity map that reproduces the CO spectrum (see Table 3), and exploring a range of sizes, inclinations, PA, and rotational velocities that produce a good match between our observations and the model.
The resulting synthetic δ values estimated from the intensity maps shown in Fig. 9 are plotted as solid black lines in Figs. 5-7 along with the δ values derived from the observations (see Sect. 3.4) for the first three bandheads. To interpret our results, we need to take into account that for an inclined ring in Keplerian rotation, the greatest shift in velocity is expected to occur along the major axis of the projected ring, while the smallest shift A50, page 8 of 21 GRAVITY Collaboration: The GRAVITY young stellar object survey 1900−2800 (0.9−9) × 10 21 15 ± 5 130 ± 30 (close to zero) is expected along the minor axis. This is exactly what is seen for the bluest observable wavelengths of each single bandhead as this wavelength range behaves similarly to an individual line (see Figs. 5-7). However, the mixture of blue-and red-shifted individual J components at the peak and the tails of the bandheads results in very small photo-centre asymmetries with respect to the continuum. The observed bluest δ values show maximum values along PA of 109 • and 135 • , indicating that the PA of the major axis of the CO emitting region must be oriented within this PA range. This is consistent with the PA = 116 • of the major axis of the dusty disk as derived from the continuum interferometric data. We therefore fixed the CO ring PA to that of the disk continuum emission, that is, we assumed that the CO and disk continuum emission have a similar PA within our uncertainties. Following the same principle, the size of the CO emitting region can be accurately derived from the maximum displacement measured along the major axis of the disk. δ always reaches its maximum value along the major axis. This value is independent of v rot and i variations (see Figs. B.1 and B.3). By varying the size of the CO ring and comparing the modelled and observed displacement, we derive a size of the CO emitting region of 0.10 ± 0.02 au (see Fig. B.7 and Table 4).
On the other hand, the displacement along orientations that are not aligned with the major axis PA can give important constraints on v rot and i. As shown in Fig. B.5, lower rotation velocities produce smaller displacements and change the shape of the displacement as a function of wavelength. Similarly, lower inclinations produce maximum displacements, which also modifies the shape of the displacement curve (Fig. B.6). Following Table 4. CO properties derived from the line displacement.
0.10 ± 0.02 70 ± 10 140 ± 30 this, we found that the observed displacements along orientations that are not aligned with the major axis PA are well reproduced assuming a v rot = 140 ± 30 km s −1 and i = 70 • ± 10 • (see Table 4).

Discussion
Our study of the CO overtone bandheads in the 2.3 µm region using GRAVITY shows that the CO properties are consistent with the emission of a ring of warm CO, located at a distance of ∼0.1 au from the star, well inside the dust sublimation radius. The emitting region is likely quite narrow. The nominal width, derived by comparing the observed flux to the model emissivity at the peak of the first bandhead is (∆R/R = 0.15). This may be a lower limit, and some contribution from a cooler, more extended region with slightly lower rotational velocity cannot be ruled out. However, this very simple model is able to reproduce the basic features of our interferometric observations within an acceptable range of physical parameters. A more complex disc model with a larger set of free parameters is therefore not required to reproduce our observations, especially at our moderate spectral resolution. 51 Oph is one of the few stars that exhibit bright rovibrational CO overtone emission, and it has been studied extensively in the past. Two spectroscopic studies using high spectral resolution data have been performed in the past (Thi et al. 2005, at R = 10 000, and Berthoud et al. 2007, at R = 25 000). The observed spectra show only small differences, mostly due to the different spectral resolution. Both teams modelled the CO emission as coming from a disc region with a radial A50, page 9 of 21 A&A 645, A50 (2021) Fig. 9. Model-predicted intensity maps of the blue-shifted (average over spectral channels at 2.2925, 2.3213, and 2.3516 µm), peak (average over spectral channels at 2.2945, 2.3237, and 2.3529 µm), and red-shifted (average over spectral channels at 2.2954, 2.3246, and 2.3543 µm) line components of the first (top panel), second (middle panel) and third (bottom) bandheads. Maps are normalised to the peak intensity and are displayed for an inclination of 70 • , and a radius of 0.1 au. The physical parameters of the model are the same as in Fig. 8. temperature and column density gradient, and fit only parts of the emission spectrum, focusing on the first bandhead. Both studies conclude that the CO is mostly emitted in a narrow, dust-free portion of the disc, very close to the central star. In particular, Thi et al. (2005) fit the spectrum with a CO disc extending between 0.15 and 0.35 au. The temperature and CO column densities are T = 2850 ± 500 K and N CO = (0.17-2.5) × 10 20 cm −2 at the inner radius, and T = 1500 K and N CO 2.3 times lower at the outer radius, respectively. The lower temperature and CO column density at the outer radius make a strong contribution of the outer region to the CO total emission very unlikely, even when the emitting area is increased. This means that the outer CO disc region might mostly affect the band tail emission. On the other hand, Berthoud et al. (2007) allowed for a much larger number of free parameters (10). Their results show an emitting region with a ratio of ∼4 between the outer and inner radii. The maximum temperature is very high (∼4000 K) and decreases steeply to ∼700 K at the outer disc radius. The column density at the inner radius is ∼3 × 10 20 cm −2 . In this model, different portions of the spectrum require rather different parameter values. However, about 80% of the flux comes from the inner half of the disc, with an average temperature of 2700 K and column density of 7.5 × 10 20 cm −2 . In both studies, the values reported at the inner CO disc radius, where most of the CO emission is most likely generated, are similar to those found in our ring modelling and confirm our main conclusions: most of the CO emission comes from a relatively warm and dense narrow region (see Table 3), very close to the star, well inside the dust sublimation radius.
Regarding the geometrical constraints, the disc inclination is poorly constrained by spectroscopic fittings alone, with a preference for very high values (∼88 • ). In a pioneering paper, Tatulli et al. (2008) used the visibilities from a single observation with the ESO VLTI instrument AMBER using the 8m UT telescopes to derive a size of the CO emission of 0.15 au and a disc inclination of 85 ± +5 −15 • . Our interferometric observations and the differential phase signal in particular give a more moderate value of i = 70 • ± 10 • , but this is still consistent within the uncertainties with all previous estimates.
The overtone CO emission provides one of the very few direct tests of the physical and chemical properties of the innermost gaseous discs. It is unfortunate that no detailed model exists for this region, and we can only be guided by considerations based on simple accretion disc models. Muzerolle et al. (2004) computed models for the inner disc around a star of T ∼ 9000 K and L = 30 L , following the procedure of solving the radiation transfer as outlined in Calvet et al. (1991). They used mean A50, page 10 of 21 GRAVITY Collaboration: The GRAVITY young stellar object survey opacities, computed assuming local thermal equilibrium, which are very likely not appropriate for the disc atmosphere, but provide a solution for the disc midplane temperature (Dullemond & Monnier 2010). At 0.10 au from the star, the midplane disc temperature for 51 Oph is around 2000 K assuming T eff = 10 000 K, somewhat lower than the CO temperature. The column density of gas at the distance R from the star can be written as a function of the midplane temperature, the mass accretion rate, and the viscosity parameter α visc (Lynden-Bell & Pringle 1974). The mass accretion rate of 51 Oph is not well known and is probably lower than <10 −7 M yr −1 (Mendigutía et al. 2011, and discussion therein). For a nominal value α visc = 0.01, the expected gas column density is therefore N H < 4.6 × 10 27 cm −2 . Assuming a CO/H of ∼10 −4 , the total column density of warm gas required by the observed CO emission is ∼4 × 10 25 cm −2 , and the CO emitting region is located above the midplane, at z ∼ 2.6H p , where H p is the pressure scale height. The density in this region is very high, about 3 × 10 15 cm −3 , sufficient to thermalize the CO vibrational levels (Thi et al. 2013b). For M acc = 10 −8 M yr −1 , the CO emitting region comes closer to the midplane, at z ∼ 1.7H p , with similar density. It is therefore likely that the CO is emitted in a column of gas slightly above and somewhat warmer than the midplane. Similar results were observed in the case of MYSO (see GRAVITY Collaboration 2020, for more details).
This column of warm CO in the intense radiation field that characterises the inner disc of 51 Oph depends on a number of chemical and physical processes, that need to be understood in detail. Recently, Bosman et al. (2019) suggested that CO is unlikely to exist in dust-depleted regions of Herbig Ae stars. However, the models do not extend to the innermost dust-free disc regions and have been computed for lower luminosity stars (∼30 L ; the luminosity of 51 Oph is 250 L ). The detections and characteristics of the vibrationally excited CO in 51 Oph (and of a few other HAeBe stars) indicate the need of computing appropriate models for the innermost dust-free regions and explore their robustness. Work in this direction is currently in progress (Gorti et al., in prep.).
On the other hand, as mentioned in the introduction, a detection of the overtone CO emission in Class II stars is rare. To the best of our knowledge, in addition to 51 Oph, it has been detected in 9 of the ∼91 objects surveyed by Ilee et al. (2014). Although the small number of detections prevents any statistical analysis, the lack of detections at luminosities below ∼20 L andṀ acc below ∼10 −8 M yr −1 is conspicuous. This is confirmed by the very low number of detections in T Tauri stars, where only three of the many objects searched so far show CO emission (including the whole Lupus sample studied by Alcalá et al. 2017; see Koutoulaki et al. 2019). Only one detection of the 15 objects has L > 16 000 L . Within this wide parameter range, the CO detections are distributed rather uniformly. These results indicate that a low stellar or accretion luminosity and/or a low mass-accretion rate prevent the formation of a sufficiently high column density of warm CO. This is confirmed by the detection of the overtone CO emission in some EXors during outbursts only (e.g. Aspin et al. 2010;Guo et al. 2020). High-mass young objects tend to have a larger fraction of CO detections than lower luminosity ones. Most recently, Pomohaci et al. (2017) detected the CO overtone emission in a subset of 7 out of 38 massive young stellar objects (MYSO). These detections seem to be confined to a narrow range ofṀ acc (10 −4.5 -10 −3.0 M yr −1 ), where they account for about 40% of the objects. As for lower luminosity objects, most upper limits to the CO bandhead luminosity are well below the detections, suggesting that the relatively small fraction of detections is not a sensitivity issue.
However, the detection of CO emission in only a few objects among many of similar properties remains puzzling. It is possible that overtone CO emission occurs only in a very limited range of physical conditions in the inner dust-free regions of discs because of a combination of factors, such as mass accretion rates, stellar and accretion luminosity, and X-ray emission from the central star that can only be understood through detailed nonlocal thermal equilibrium models. Even if rare, the detection and characterisation of this emission can give us the only constraints on the innermost gaseous regions of discs, in which accretion onto the star occurs and from which winds are launched.

Conclusion
We reported our results on the Herbig Ae/Be star, 51 Oph, using VLTI/GRAVITY interferometric observations in the K band at high spectral resolution. Our main conclusions are listed below.
-51 Oph shows prominent CO bandhead emission. The first four ro-vibrational transitions are clearly detected in the spectrum, as well as an indication of the fifth. Our GRAVITY observations show a clear increase in visibility signal within the first four bandheads, indicating that the CO emitting region is more compact than the total continuum emission. Clear differential phase signatures are observed at the positions of the first four bandheads. -By fitting the continuum visibilities, we derive a stellar radius of R * = 0.4 ± 0.1 mas (R * = 10.6 ± 2.6 R ) and a dusty disc with radius R disc = 4.0 ± 0.8 mas (R disc = 0.50 ± 0.01 au) at an inclination and PA of 63 • ± 1 • and 116 • ± 1 • . -By modelling the CO bandhead emission, we created intensity maps of the CO emission that were used to reproduce the observed spectrum and the differential phase signatures to broadly constrain the physical properties (T and N CO ), kinematics (v rot ) and geometry (size, i, PA) of the CO emitting region. -Our modelling shows that the bandhead CO emission in 51 Oph is due to an overlap of optically thick individual J components. This makes it difficult to accurately constrain the T and N CO , as a change in T can be compensated for by a change in N CO . Nevertheless, in order to reproduce the CO spectrum a warm (T ∼ 1900-2800 K) and dense (10 20 -10 21 cm −2 ) gas is needed in agreement with previous results (e.g. Thi et al. 2013a;Berthoud et al. 2007). -The synthetic intensity maps of the CO emission were used to reproduce the observed CO displacements with respect to the photocentre of the continuum. From the modelling and analysis of the observed displacement an estimate of the size, inclination, PA, and rotational velocity of the CO emitting region was derived. A CO ring of radius R CO ∼ 0.1 au reproduces the maximum displacement at the bluest wavelength of the three first bandheads well. In addition, we find that the CO PA and inclination is consistent with that of the dusty disk, indicating no major misalignment between the dusty disc and the CO emitting region.
In this section we illustrate with examples how the spectrum of the CO bandhead emission changes when different free parameters are varied. We start by showing that the CO bandhead emission in 51 Oph is due to an overlap of optically thick individual J components (Appendix A.1), followed by a description of the spectral variation when T and N CO are varied (Appendix A.2), and the intrinsic line width and v rot sin i values (Appendix A.3).

A.1. Optical depth
The CO bandhead emission in 51 Oph is due to an overlap of optically thick individual J components. The two last bandheads (υ = 5-3 and 6-4) are optically thinner than the first three bandheads (υ = 2-0, 3-1, and 4-2). This is illustrated in Fig. A.1, where the optical depth (τ) of the first five CO bandheads versus wavelength is plotted at a spectral resolution of R = 100 000. As a consequence, the temperature and CO column densities cannot be accurately constrained because an increase or decrease in T can be compensated for by an opposite change in column density. Figures A.2 and A.3 illustrate that an increase (decrease) in temperature can be compensated for by a decrease (increase) in CO column density so that most of the individual J components remain optically thick. When T is decreased and N CO is increased the optical depth increases from the most strongly red-shifted towards the most strongly blue-shifted wavelengths (Fig. A.2). Eventually, all the J-components reach the equivalent Planck function. In contrast, a gradual increase and decrease of T and N CO , respectively, decreases the optical depth from the highest to the lowest υ band components until they eventually become optically thin (Fig. A.3).
A.2. Temperature and column density  A.3. Intrinsic line width and v rot sin i In addition to the physical parameters, the CO bandhead profiles also depend on the assumed intrinsic line width of the individual J components. To better illustrate this dependence, we produced very high spectral resolution spectra (R = 100 000) of the individual J components contributing to the tail of the first bandhead in the wavelength range between 2.302 and 2.309 µm for two different ∆v values (see Fig. A.7). This figure shows that an increased ∆v value broadens the profile of the J components while the emitting area remains the same. As a consequence, an increase in ∆v produces a blend of J components. This effect translates into variations in line flux of the CO bandheads as shown in Fig. A.8. This figure shows the continuum-subtracted spectrum of 51 Oph (black dots) normalised to the peak of the first bandhead, along with three different synthetic spectra. The modelled spectra were computed for the same T , N CO , and v rot sin i, but the intrinsic line width was varied from 10 to 20 km s −1 . The increased blend of J components for higher values of ∆v results in a gradual increase in the bandhead flux.
On the other hand, changes in v rot sin i also produce changes in the CO banhdead profiles. At our moderate spectral resolution, changes in v rot sin i mostly affect the blue-shifted shoulder of the CO first bandhead and its peak position. Higher values of v rot sin i produce more pronounced blue-shifted shoulders, and shift the peak of the bandhead towards red-shifted wavelengths (see Fig. A.9).      .The data are smoothed at a resolution of R = 1000. When the data are smoothed the first bandhead becomes smaller because the line peak has only a few points, which causes the second and third bandhead to exceed than the unsmoothed spectra when we normalise.

Appendix B: CO line displacements
The CO line displacement with respect to the continuum photocentre depends on the physical properties of the CO (T, N CO , and ∆v), its geometry (i, PA), kinematics (v rot ), and on the size of the emitting region. For the particular case of 51 Oph and the moderate spectral resolution of our observations, the line displacement does not significantly change within the range of acceptable values of T , N CO , and ∆v, but it can help to better constrain the i, PA, and size of the CO emitting region. The main displacements are always observed at the most strongly blue-shifted wavelength, as the mixture of individual J-components at the peak and red-shifted tails of the bandheads results in small to zero photo-centre shifts. In an inclined ring in Keplerian rotation, the maximum and minimum displacements are achieved along the major and minor axis, respectively, reaching intermediate values. Along any fixed PA, the CO line displacement is independent of v rot , and it reaches its maximum value along the major axis (Figs. B.1, B.2, and B.5). By probing the line displacements along different PAs, an estimate of v rot can be found. On the other hand, for PAs other than that of the ring major axis, the line displacement depends on the ring inclination (see Figs. B.4 and B.6). This can be used to constrain the ring inclination by comparing the observed and synthetic line displacements along different PAs, assuming different values of i.
Similarly, the line displacement depends on the size of the emitting region. As before, the main displacements are observed for the most strongly blue-shifted wavelengths, especially along the major axis PA. For fixed values of i, and v rot , this can be used to derive the size of the CO emitting region (see Fig. B.7).