Open Access
Issue
A&A
Volume 690, October 2024
Article Number A209
Number of page(s) 9
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202449641
Published online 08 October 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Massive stars (M > 8 M) are among the most influential objects of a galaxy. They exert significant influence over the motion of surrounding material and shape a galaxy’s chemical composition by generating heavier elements. Despite their importance, their formation and evolution are still not well understood, particularly in terms of how these massive objects sustain mass accretion in the face of substantial radiation pressure.

Massive stars are believed to form following the accretionejection paradigm (e.g., Rosen et al. 2019; Klassen et al. 2016; Hosokawa et al. 2010). The short lifetimes of discs surrounding massive young stars (<105 yr; e.g., Kuiper & Hosokawa 2018) seem to be coupled with the strong photoevaporation wind driving most of the disc dispersal (e.g., Owen et al. 2011; Ercolano & Rosotti 2015).

Ilee et al. (2013) performed kinematic modelling of the CO bandhead emission towards 20 massive young stellar objects (MYSOs) and found evidence of gaseous discs that ‘live’ at distances as close as a few (sub-)au from the central protostars (see also, Bik & Thi 2004; Chandler et al. 1993, 1995). Direct spatial information of such small-scale hot discs has been limited to a couple of studies of individual objects via hot dust emission or the CO bandhead (Kraus et al. 2010; GRAVITY Collaboration 2020a). Koumpia et al. (2021) added a sample of six MYSOs to this dedicated effort. More recently, the Na I doublet and Helium emission have been found to also originate from such small scales in the interface between the dusty disc and protostellar objects (Koumpia et al. 2023).

The CO first overtone (or ‘bandhead’) emission is not commonly observed in the spectra of MYSOs at low to medium resolutions. Instead, it is detected in only about 25% of their spectra (Ilee et al. 2018). Until now, there has not been a study combining kinematic modelling of the molecular gaseous component (CO) using high spectral resolution observations (R∼30 000) with direct spatial information. The only case that used a similar approach was limited to a spectral resolution of R ∼4000 (GRAVITY Collaboration 2020a). Yet, it is crucial to provide the ultimate proof that not only do small-scale gaseous discs surrounding MYSOs exist, but they are the pathways to active accretion.

In addition, observationally, all MYSOs appear to show Brγ emission in their spectrum (Bunn et al. 1995; Frost et al. 2021). The origin of the Brγ emission has been previously attributed to ‘disc-wind’ (Drew et al. 1998) or jet (Caratti o Garatti et al. 2016). Still, distinguishing between the two is not always straightforward, especially when the spectral resolution is not adequate to do so (Koumpia et al. 2021). Given the influence of both mechanisms on the final mass of the central object, the combination of high spectral and spatial resolution is crucial. Interferometry and the unique information provided by the differential phase of GRAVITY in the HR mode are powerful tools for telling those processes apart, as they enable astrometric solutions at a high angular resolution.

In this paper, we focus on the MYSO G033.3891 (IRAS 18490+0026), which with a lower-mass limit of M > 12 M and luminosity of L~1.3×104 L is a good example of this class of objects. Moreover, G033.3891 is one of the very few MYSOs known to host an accreting gaseous disc, as traced via the CO bandhead emission (Wheelwright et al. 2010), and it is bright enough for high angular resolution observations in the K-band (K=7.m2$K = 7\mathop .\limits^m 2$) with GRAVITY on the Very Large Telescope Interferometer (VLTI). Ilee et al. (2013) modelled the observed Cryogenic high-resolution Infrared Echelle spectrograph (CRIRES) spectrum of the source around the bandheads (R~30 000) and found that the best solution was achieved for a disc with an inner radius of only ~2 au. The inclination of the disc was determined to be approximately ~40º.

This work presents the first spatially resolved observations of the hot dust as traced with a 2.2 µm continuum and the ionised (Brγ) and molecular (CO) gas content on the source using near- IR interferometry. Our observations allowed us to apply detailed geometrical models and to constrain the origin of the Brγ emission spatially using astrometry. For the first time, we constrain the origin of the CO bandhead emission towards G033.3891 both spatially and kinematically.

Table 1

Coordinates, bolometric luminosity, K-band magnitude, and distance of G033.3891 as taken from the RMS database. The mass and effective temperature of the source are also listed.

2 Methods and observations

2.1 GRAVITY observations and data reduction

The MYSO MSX6C G033.3891+00.1989 (IRAS 18490+0026; RA = 18h51m33.8s, Dec = +00º29′51.″1 [J2000]; Table 1) was observed on the 30 May 2021 with the GRAVITY interferometric instrument (GRAVITY Collaboration 2017; Eisenhauer et al. 2011) on the VLTI using the four 8.2-m Unit Telescopes (UTs). GRAVITY covers the K-band (1.99–2.45 µm). The deeply embedded nature of the source and the absence of suitable guiding stars on-axis required the use of the Coudé CIAO system for natural star guiding off-axis. The observations were taken using the high spectral resolution in combined polarisation mode (HR; R~4000, Δv~75 km s−1). The observing run1 was executed under good atmospheric conditions, with a seeing varying between 0.6″and 0.8″ and τcor ~ 5 ms. The obtained uv-coverage is presented in Figure 1. The fringes did not show consistent stable behaviour throughout the entire run, yet a high fringe tracking ratio of 81–97% could still be obtained.

HD 174606 (RA = 18h51m 16.67s, Dec = +00º49′59.″4 [J2000]) was observed as a standard calibrator under good conditions similar to the science object and showed stable fringes with a fringe tracking ratio of 100%. HD 174606 is characterised by a K-band magnitude of 7m, a spectral type of F3III, and a uniform disc size in the K-band of 0.157 mas (i.e. spatially unresolved, JMMC SearchCal; Bonneau et al. 2011).

The observations were reduced and calibrated using the standard GRAVITY pipeline recipes as provided by ESO (version 1.5.4). The instrumental transmission in the spectrum was corrected via the pipeline. The spectra were also corrected for tellurics using the HITRAN models (Rothman et al. 2009) via PMOIRED2 (similar to Molecfit, Smette et al. 2015; Mérand 2022). To address continuum irregularities and normalise the spectrum, we applied up to a third-order polynomial fitting.

thumbnail Fig. 1

uv-coverage of G033.3891 obtained using GRAVITY/VLTI on UTs. The position angle (in degrees) and projected baseline length (in metres) of each baseline pair are indicated in the legend.

thumbnail Fig. 2

Normalised spectrum for G033.3891 as observed with GRAVITY. The telluric corrected spectrum shows the Brγ and the CO bandheads in emission.

2.2 Interferometric observables

The interferometric observables towards G033.3891 delivered by GRAVITY consist of spectral information, calibrated visibilities, and differential and closure phases in the K-band (2.0–2.4 µm). The K-band spectrum (Figure 2) shows a prominent Brγ emission line at 2.167 µm and the CO bandhead emission around 2.3–2.4 µm. Based on the presence of Na I in the low-resolution spectra of the source (Pomohaci et al. 2017), the GRAVITY spectra were also inspected for the presence of the Na I doublet emission but without a clear detection. The root-mean-square of the spectrum across the full spectral window is 0.034. An inspection of the calibrated visibilities and phases around the spectral lines of interest allowed us to extract some qualitative information regarding the relative size and symmetric or asymmetric nature of the line emitting regions (Brγ, CO) compared to the 2.2 µm continuum emission.

The absolute visibility – also referred to as visibility amplitude V(f,λ) – at the peak intensity of the Brγ line shows an increase at all six observed baselines compared to the visibilities of the surrounding continuum (Figure 3). At the longer three baselines (~90–130 m), this increase is greater than 5% and reaches up to ~13% at the longest one (U4U1; 127 m). This is higher than the typical uncertainties introduced by calibrations and the transfer function. The observed increase for the other three baselines is comparable to typical uncertainties (2–5%; excluding residual calibration effects). The peak of the line emission coincides with the velocity location corresponding to the peak of the visibility, suggesting that the observed variations at all baselines are due to a real geometrical effect. Larger visibilities along emission lines are suggestive of a smaller emitting area compared to that of the dust continuum. Therefore, the Brγ emitting region is more compact compared to the 2.2 µm continuum.

The differential phases as obtained with GRAVITY show variations around the Brγ emission up to 2.8°. The observed changes suggest an offset between the photocentre of the 2.2 µm continuum emission and that of the Brγ emitting region. On the other hand, no variations were observed in the observed closure phases around the line emission, suggesting an axisymmetric nature of the ionised gas emitting region.

Examining the calibrated visibilities around the CO bandhead emission, there is an increase compared to the surrounding continuum at two baselines (U4U1, U4U2), which is approximately 5–7% at the longest baseline. The consistency in the location of the peak visibility with that of the emitting components of the CO bandheads demonstrates that the observed variation is real. The absolute visibilities indicate that the CO bandhead arises from a more compact area compared to the dust continuum emission. In Sect. 3.1, we employ detailed geometric modelling to quantitatively constrain the sizes of the continuum and the Brγ and CO bandhead emission.

There are no significant variations in the phases (closure or differential) around the CO bandhead emission compared to the continuum, suggesting that the molecular gas component shares the same photocentre as that of the continuum and that the emitting region is symmetric. This is only the third case of an MYSO known in the literature with spatially resolved observations of the CO bandhead emission (NGC 2024 IRS 2, IRAS 13481-6124; GRAVITY Collaboration 2020a; Koumpia et al. 2023).

3 Results

3.1 Geometric modelling

We employed geometric modelling of the GRAVITY observations of G033.3891 to characterise the size and morphology of the hot dust emission, as traced by the 2.2 µm continuum, the Brγ, and the CO bandhead emission, all captured within the K- band spectrum. The present analysis utilised the specialised tool PMOIRED (see also Sect. 2.1) designed for parametric polychromatic modelling of spectro-interferometric data. The code enables the simulation of interferometric observables for diverse models of standard brightness distribution blocks (e.g. Gaussian, disc, ring) and any combination thereof. The PMOIRED model composition assumes that the three-dimensional flux can be expressed as F(x,y,λ)=I(x,y)×S(λ),$F(x,y,\lambda ) = I(x,y) \times S(\lambda ),$(1)

where x and y represent the angular coordinates on the sky and λ is the wavelength. Here, I represents the achromatic (modeldependent; wavelength-independent) spatial distribution of the component, while S represents its spectrum (flat or smooth function for the continuum emission; wavelength and line width - dependent for the line emission). This formulation allows for the separate consideration of spatial structure (I) and spectral characteristics (S) when modelling astrophysical objects (see also, Monnier 2003).

Subsequently, the user can determine the model that best aligns with the observations by evaluating the optimal fit through a least-square approximation and a bootstrapping estimation of the uncertainties, as retrieved from the covariance matrix. The correlation matrix is crucial for ensuring that the parameters are independent of each other. A high correlation (e.g. greater than 90%) suggests that the parameters are degenerate, either due to an inadequate model parametrisation or because the data cannot distinguish between the two parameters. Our best-fit results, as presented below, do not suffer from degenerate solutions.

Our visibility models follow the standard approach for modelling YSOs, including a point source representing the unresolved star in the 2.2 µm emission together with an extended emission component (e.g., Kraus et al. 2008; Tatulli et al. 2008). This approach ensures a comprehensive representation of the circumstellar continuum emission around an MYSO.

3.1.1 K-band continuum

We first constrained the size and structure of the K-band emission. Our fitting process encompasses various combinations of brightness distributions, including Gaussian, disc, and ring profiles (see Table 2). To limit the degree of complexity during geometric modelling, in the case of a ring or a disc, a face-on orientation was adopted.

The optimal fit to the 2.2 µm continuum emission is characterised by a reduced χ2 value of 1.2, and it was achieved by adopting a ring of an inner diameter of 1.95±0.08 mas and an outer diameter of 2.8±0.06 mas. The reported size represents the best-fit result obtained from fitting the entire K-band continuum, assuming that its size and geometry exhibit minimal variation across the observed spectral range. At a source distance of 5 kpc (determined kinematically in RMS survey; Lumsden et al. 2013), the inner and outer radius of the ring were determined to be 5 au and 7 au respectively. The inner ring radius is in perfect alignment with the predicted dust sublimation radius (Rdust) of the source, which was computed to be 5 au (or else 1 mas at 5 kpc) for a source ~1.3×104 L, and after adopting a dust sublimation temperature of 1300 K (e.g. olivines, silicates; see Isella & Natta 2005; Kama et al. 2009; McClure et al. 2013). The dust sublimation radius (Rdust) was estimated from the dust sublimation temperature, Tsub, and the stellar luminosity, L*, using Rdust =L*16πσSBTsub4,${R_{{\rm{dust }}}} = \sqrt {{{L*} \over {16\pi {\sigma _{SB}}T_{sub}^4}}} ,$(2)

where σSB is the Stefan-Boltzmann constant.

The simultaneous fit of the visibilities and the closure phases (−2.5° < CP < 3°; typical uncertainties <1.5°) of the continuum emission required the introduction of an asymmetric emission. In PMOIRED, it is possible to simulate deformed or asymmetric components analytically by the introduction of a ‘slant’ in a specific direction. The geometric model that best fits the interferometric observables of the continuum required the inclusion of a slant of 1.2 at a position angle (PA) of ~70±6° (Fig. 3).

We conclude that the 2.2 µm emission traces the annulus corresponding to the inner part of the hot dusty disc surrounding this MYSO. This conclusion is based on the temperatures of the material that the 2.2 µm emission is sensitive to (~1300 K), the ring-shaped morphology of its emitting region derived from our geometrical models, and the observation that the determined inner radius is consistent with the predicted dust sublimation radius of this system.

thumbnail Fig. 3

Model image of G033.3891 of the 2.2 µm continuum emission (thick ring) combined with the Brγ (compact Gaussian) and CO bandhead molecular emission (thin bright ring). The modelled image of the CO emission is the average of the individual images of the CO bandheads extracted at their central wavelengths (2.29, 2.32, 2.35, and 2.38 µm). The total geometric model (red line) fits well all interferometric observables (closure phases, T3PHI; differential phases, DPHI; visibilities, V; normalised flux, NFLUX) of the continuum, Brγ, and CO bandheads. PMOIRED normalises the flux from both the data and the model using the same algorithm. The grey shades represent the systematic errors of the interferometric observables, taking into account the calibration and transfer function uncertainties. The colour scheme of the different baselines follows Figure 1.

Table 2

Parameters of the best fit to the interferometric observables of the 2.2µm, Brγ, and CO emission towards G033.3891.

3.1.2 Brγ emission

We employed a similar approach for constraining the brightness distribution of the Brγ emission as for the continuum emission. Specifically, we explored diverse combinations of brightness distributions (Gaussian, disc, ring). To limit the number of free parameters during the fitting process of both the line and continuum emissions, we constrained the size and geometry of the continuum emission to be identical to the best-fit result for 2.2 µm described in Sect. 3.1.1.

To derive the size of the line emitting region, the actual visibilities of the line needed to be extracted by subtracting the continuum contributions. This was done by using Eq. (3), which describes the total visibility (Vline+cont) in terms of continuum and line visibilities (V) and fluxes (F) for multi-component sources, solving for (Vline) as follows: Vline + cont =Vcont ×Fcont +Vline ×Fline Fcont +Fline .${V_{{\rm{line }} + {\rm{ cont }}}} = {{{V_{{\rm{cont }}}} \times {F_{{\rm{cont }}}} + {V_{{\rm{line }}}} \times {F_{{\rm{line }}}}} \over {{F_{{\rm{cont }}}} + {F_{{\rm{line }}}}}}.$(3)

The best-fit for the Brγ emission corresponds to a Gaussian brightness distribution (with a reduced χ2 of approximately 3.2), resulting in a full width half maximum (FWHM) of 0.66±0.06 mas (equivalent to 3.3 au). To accurately reproduce the observed changes in differential phases along the Brγ emission, we allowed the central position of the Gaussian distribution to be a free parameter, rather than fixing it to the same location as the photocentre of the continuum. Introducing an offset between the photocentre of the Brγ emission and that of the continuum is necessary to accurately reproduce the observed variations in differential phases, as shown in the top panels of Figure 3.

The best fit results in an offset between the photocentre of the Brγ emission and that of the compact continuum of about 0.4±0.1 mas (or 2 au) towards the north-east at PA~41.4±3.3°. The lack of closure-phase variations along the Brγ emission suggests that there is no need to introduce an asymmetry in its brightness distribution.

The geometrical models demonstrate that the Brγ emission comes from an area that is significantly more compact compared to the continuum and well within the inner radius of the hot dusty disc. The observed offset between the Brγ photocentre and that of the continuum is found in the same direction as the brighter part of the disc and is indicative of some degree of diversion from a face-on geometry.

3.1.3 Molecular gas emission

To model the emission of the CO bandheads (2.29, 2.32, 2.35, and 2.38 µm), we followed a similar approach to that used for Brγ and the continuum. To minimise the number of free parameters during the fitting, we adopted the size and geometry of the continuum emission to match the best-fit result obtained for the 2.2 µm emission described in Sect. 3.1.1. We then tested the three standard brightness distribution models (Gaussian, uniform disc, ring) to fit the spectro-interferometric observables around the CO emission lines.

The best-fit model (reduced χ2 ~ 3.0) was achieved for a ring brightness distribution with an inner diameter of 1.10±0.12 mas (5.5±0.6 au) or else an inner radius of 0.55±0.06 mas (~2.8±0.3 au). The preference for the ring model is also justified based on physical expectations. The CO bandheads are generally understood to originate from a circumstellar disc with an inner cavity (e.g., Carr et al. 1993; Tatulli et al. 2008; Ilee et al. 2013), naturally forming a ring-like structure.

Constraining the outer radius of the ring appeared to be the most difficult parameter to implement in the fitting process. When considering the ‘thickness’ of the ring (in PMOIRED, this reflects the radial extension of the ring) as a free parameter, convergence of the code led to degenerate solutions, oscillating between an infinitely thin (model parameter ‘thick’ < 0) or infinitely thick (‘thick’ > 1) ring. To address this issue and proceed with the fitting process, we ran multiple fits by constraining the thickness of the ring to multiples of 0.1 (up to one) at a time. Upon fixing the thickness at 0.2 (corresponding to an outer radius of 3.4 au), the best fit was achieved. By finally fixing the thickness of the ring at this value, we could further facilitate a robust fit to all observables, leaving the rest of the parameters free. This approach resulted in a final geometry that does not violate the mathematical solution of a physical system.

We note that the goodness of the fit is not as strong for three baselines (U3U2, U4U3 and U2U1) compared to the others. These baselines correspond to shorter baseline lengths (46– 55 m), which trace the larger spatial scales, and it likely reflects the difficulty encountered during the fitting process in accurately modelling the outer radius of the ring. No significant changes in the differential or closure phases along the CO bandheads were observed, implying that the CO emission shares the same photocentre as the continuum and that there is no need to implement an asymmetry on the flux distribution of the emitting area.

The final geometrical model is shown in Figure 3, and it combines the continuum and the Brγ and the CO bandhead emission (see also Table 2). We conclude that the CO emission is similar in size to the Brγ emitting region, and both are significantly more compact compared to the dust sublimation radius, therefore tracing the star-disc interface. While our models have successfully constrained the inner radius of the CO ring structure, the robust determination of the molecular content’s outer extent remains challenging due to instrumental limitations, such as the spatial resolution, the limited interferometric field of view, and the wavelength coverage of GRAVITY.

3.2 Kinematic modelling of the CO bandhead

In addition to the direct spatial information we retrieved from the GRAVITY dataset in Sect. 3.1.3 regarding the CO size and distribution, we obtained an independent measurement via kinematic modelling of high spectral resolution spectra. The CO first overtone emission has been observed towards G033.3891 with VLT/CRIRES (Ilee et al. 2013), providing a relatively high spectral resolution view of the bandhead features (R ~ 300 000). In that study, the v = 2–0 bandhead was well reproduced by a model of a thin Keplerian disc at an inclination of 40° with an inner radius of 2.1 au3 (Figure 4).

We note that the stellar parameters describing G033.3891 used in the previous work to fit the CRIRES spectra were taken from a preliminary collation of information in the RMS survey (Lumsden et al. 2013). Since that work, the improved estimate of the bolometric luminosity of the source has increased from 1.0 × 104 L to 1.3 × 104 L. The expected change in stellar mass resulting from this difference is of the order 0.5 M (note the estimated mass of 12.3 M; Ilee et al. 2013), which would correspond to changes in orbital velocity of less than 1.5 km s−1 beyond 2 au. Since this is far less than the spectral resolution of the observations (Δv ~ 10 km s−1), this slight change does not affect the final results of the inner radius of the gaseous disc.

The inner size of the gaseous disc traced with CO bandhead emission when derived independently using spatial interferometric observations (2.8 au) is in alignment with the size derived when using kinematic models applied to high spectral resolution spectra (2.1 au). The inner radius derived from both methods is independent of the inclination angle, as it represents an absolute distance from the central star rather than a projected value. While different orientations can influence detailed interpretations, the primary focus of this study is the consistency between the inner radius of the CO obtained from kinematic modelling and that derived from geometric modelling. This key result is robust and independent of the assumed inclination.

The outer radius of the CO gaseous region could not be accurately constrained by either kinematic or geometric modelling. The kinematic modelling set the outer radius as far out as 3200 au. This value is a result of extrapolation using the temperature exponent, which in turn is poorly defined when fitting only the 2–1 bandhead. Combined with the inability of the models to set an uncertainty, we conclude that there is no real evidence that the CO gaseous ring extends that far out. The geometric models applied to GRAVITY observations suffered from a similar uncertainty regarding the fit of the outer radius.

Fitting the rest of the bandheads is crucial for constraining the outer radius of the disc and confirming its true extent, especially since the first overtone emission is known to be optically thick (e.g., GRAVITY Collaboration 2021). The GRAVITY spectrum could aid in this task, but the combination of it being heavily affected by the atmosphere and being of a lower resolution compared to the CRIRES spectrum makes it very challenging. Probing the fundamental CO emission in the mid-infrared wavelengths (e.g. with MATISSE) in the future is therefore important to constraining the outer radius of the disc.

Our study highlights the powerful impact of combining spatial and kinematic data. This approach helped us pinpoint with great accuracy where the CO bandhead emission originates from, revealing the closest region of gas accretion around this MYSO, well within the sublimation radius. Notably, our analysis of G033.3891 marks the second instance of such an investigation into MYSOs, contributing significantly to the understanding of how these massive young stars gather mass.

thumbnail Fig. 4

VLT/CRIRES observation of the CO bandhead in black (R ~ 30 000) overlaid with the best fitting disc model in red (adapted from Ilee et al. 2013).

thumbnail Fig. 5

Brγ photocentre shift displacements for different velocity channels calculated from the differential phases after correcting for the continuum contribution. The black contour corresponds to the inner ring size of the continuum emission.

3.3 Spectro-astrometry

For the Brγ emission, for which clear variations along the differential phases were observed, it is possible to obtain an astrometric solution. An astrometric solution allows for the precise measurement of both the position and the motion of the region emitting the Brγ line, informing us about its spatial distribution and dynamics. To constrain the displacement of the photocentre for each velocity channel of the Brγ emission, we used the calibrated differential phases, following the equations described in Caratti o Garatti et al. (2016). To ensure that the photocentre offsets correspond to the line-emitting region, the continuum contributions were subtracted from both the visibilities and the differential phases. The displacement of the photocentre of the emission at any given wavelength was then computed as δ=ΔΦλ2πB,$\delta = - \Delta \Phi {\lambda \over {2\pi B}},$(4)

where ∆Φ is the continuum subtracted differential phase and B is the length of the baseline (Lachaume 2003; Le Bouquin et al. 2009). The displacements were calculated for all available baselines and for each spectral channel over the line profile where line-to-continuum ratio is greater than 10%. The velocities were calculated with respect to the local standard of rest, which is 26.85 km/s.

Figure 5 shows that the displacement of the photocentre of the blueshift and redshift velocities generally follows a distribution where the blueshifted components are mostly located towards the north-east and the redshifted components are mostly located towards the south-west. The exception to this trend is the three redshifted data points towards the east and one blueshifted data point towards the west. We note that the observed blend of red- and blueshifted velocity components could be due to the spectral resolution. The reported velocities range between −75 and 125 km/s. If the observed velocity discrepancies are real, they add an extra complexity to interpreting the physical origin of the emission, which altogether seems to follow the asymmetric direction observed in the continuum emission. We found the PAs of the continuum asymmetry and the offset from the photocentre of the Brγ emission to be 70° and 41°respectively.

The global astrometric solution follows the geometry we obtained in Sect. 3.1.2, where the brightness distribution of the continuum and the displacement of the Brγ are both stronger towards the north-east. This solution for the Brγ emission suggests that if it is attributed to a jet, the jet deviates from a perfect edge-on geometry. An inclined disc could explain the observed asymmetry of the assumed jet, as well as the decrease in the brightness of the dusty disc in the south-west at PA of 250° (i.e. an obscured inclined disc).

In conclusion, our findings suggest that there is a mix of low and high velocity moving ionised gas towards the east-north, west-south orientation, following an orientation similar to that of the dusty disc. The lack of a systematic trend in the velocity field indicates that the Brγ emission is probably a result of both a disc-wind and a jet mechanism within the inner 0.5 mas, while the emission becomes more collimated outwards (i.e. outflow or jet origin).

4 Discussion

4.1 Brγ emission

The origin of the Brγ line in (M)YSOs is a matter of long-term debate, with various proposed physical mechanisms that could give rise to the emission. These mechanisms include the accretion of matter onto the star, outflowing material originating from a disc wind, an extended wind or outflow, and a jet (Eisner 2007; Garcia Lopez et al. 2015; Stecklum et al. 2012; Koumpia et al. 2021).

Our geometric modelling combined with the astrometry obtained from the HR GRAVITY/VLTI observations of G033.3891 suggest that the main component of the Brγ emission has a bipolar origin. The mix of low and high velocity components and blue- and redshifted emission at similar distances from the central star indicate that the mechanism responsible for the Brγ emission is rather a combination of disc-wind and a jet. We report velocities that exceed ~100 km/s while noting that jets are expected to reach typical velocities as high as several hundred kilometres per second (e.g., Torrelles et al. 2011).

The general observation of increasing velocities with increasing distance from the central source can exclude a pure accreting origin of the emission (e.g. magnetospheric accretion). The estimated size of the Brγ emission is also larger by a factor between 4 and 40 compared to the Alfvén radius, assuming an approximation of the magnetospheric radius to be in the range of 2–20 R* (see, Hartmann et al. 2016; GRAVITY Collaboration 2020b; Zhu et al. 2024). The size, the velocity information, and the geometry of the Brγ emission does not support magnetospheric accretion being its sole responsible underlying mechanism.

The visibility and differential phases observed at high velocities suggest that the outflowing material spans up to ~3 au. Therefore, the emission is quite compact. We note that for a jet, one would expect material that is more extended, reaching a few tens of au (Caratti o Garatti et al. 2016), which is within the range of scales that GRAVITY can trace. Initial jet opening angles of ~30 degrees are seen in collimated material around low-mass YSOs. The astrometry of Brγ around G033.3891 indicates a blend of both collimated and more dispersed emissions, further supporting a joint origin of jet and a disc-wind as the underlying cause of the Brγ emission.

4.2 Molecular gas emission

The work in this paper reports the third case of an MYSO showing spatially resolved CO bandhead emission GRAVITY Collaboration (2020a); Koumpia et al. (2023). It is also the second study ever to combine geometric and kinematic modelling to trace the emission’s exact origin.

Geometric models have revealed that the CO emission may stem from a ring tighter to the central protostar compared to the dust. The CO is somehow protected from undergoing photo-dissociation, a process that would typically occur due to exposure to intense ultraviolet (UV) photons by the central stellar object at such a proximity. This is probably due to the CO gas emitting disc’s self-shielding.

The CO ring that we found has an inner radius of around 2.8 au, which within uncertainties matches our results from the kinematic modelling of the CRIRES spectra (R~300 000) towards the source (2.1 au, see also, Ilee et al. 2013). The kinematic fit assumed a circumstellar disc in Keplerian rotation and found an inclination of 40° for G033.3891. This could cause the shadowing we observed in the dusty disc via the interferometric fitting. To achieve more robust conclusions, it would be necessary to include the inclination and position angle of the disc in the fitting process. This could also refine the size and error estimates, but a fuller uv-plane is necessary to aid these efforts. Constraining the outer radius of the gaseous disc via the CO bandhead emission requires the observation of optically thin emission, which is better traced with higher-order CO bandheads (e.g. third-order traced by CRIRES+).

We note that the most important finding of the current study concerns the constraint of the inner radius of the gaseous CO disc, as it proves active accretion via Keplerian motion, at a location only a few aus away from the central MYSO.

5 Conclusions

In this paper, we have presented the geometric modelling of the first K -band interferometric observations of the MYSO G033.2891 (GRAVITY/VLTI) combined with kinematic modelling of its former high spectral resolution CRIRES observations of the CO bandhead emission. The MYSO G033.2891 is the third known to have spatially resolved CO (together with IRAS 13481-6124 and NGC 2024 IRS 2), adding an object to the sparse sample of such observations. The present CO analysis marks only the second investigation of its kind within MYSOs, reinforcing the disc-like nature of accretion at MYSOs at the smallest scales. In addition, the differential phases around the Brγ emission allowed us to perform spectro-astrometric analysis and to further investigate its physical origin.

The main conclusions concern the sizes of the 2.2 µm, the CO, and the Brγ emissions and are summarised as follows:

  1. The 2.2 µm continuum emission arises from a ring with an inner radius consistent with the computed dust sublimation radius (~2 au).

  2. The inner radius of the CO emission as determined spatially arises from an area significantly smaller (~2.8 au) than the dusty disc (~5 au) and, within the errors, matches the inner radius determined independently from the kinematic modelling assuming Keplerian motions (~2.1 au).

  3. The astrometric solution of the Brγ emission suggests a mixture of collimated and scattered emissions without apparent systematics in the velocity field. This finding suggests that both a jet and a disc-wind are jointly responsible for generating the observed Brγ emission.

The importance of CO bandhead emission in understanding accretion and ejection in MYSOs is increasingly being acknowledged. To fully grasp its extent, combining spatial data with high-resolution spectra covering the optically thin bandheads is crucial. Including species such as Na I, CO, and Brγ can refine kinematic information by constraining the velocity gradient of the neutral, ionised, and molecular content of the star-disc interface. As individual MYSO studies multiply, emphasising larger, statistically significant sample sizes will become apparent (e.g. CRIRES+, GRAVITY+). Spectro-astrometry is a powerful tool for constraining the physical origin of emission lines, making it essential for future studies to leverage this technique for a deeper understanding of the accretion and ejection processes around MYSOs.

Acknowledgements

We thank the anonymous referee for their detailed feedback, which greatly improved the quality of this manuscript. We also acknowledge Antoine Mérand for his assistance with PMOIRED. J.D.I. acknowledges support from an STFC Ernest Rutherford Fellowship (ST/W004119/1) and a University Academic Fellowship from the University of Leeds. Based on observations collected at the European Southern Observatory under ESO programme 0105.C-0686(D) (GRAVITY).

References

  1. Bik, A., & Thi, W. F. 2004, A&A, 427, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Bonneau, D., Delfosse, X., Mourard, D., et al. 2011, A&A, 535, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Bunn, J. C., Hoare, M. G., & Drew, J. E. 1995, MNRAS, 272, 346 [NASA ADS] [Google Scholar]
  4. Caratti o Garatti, A., Stecklum, B., Weigelt, G., et al. 2016, A&A, 589, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Carr, J. S., Tokunaga, A. T., Najita, J., Shu, F. H., & Glassgold, A. E. 1993, ApJ, 411, L37 [NASA ADS] [CrossRef] [Google Scholar]
  6. Chandler, C. J., Carlstrom, J. E., Scoville, N. Z., Dent, W. R. F., & Geballe, T. R. 1993, ApJ, 412, L71 [NASA ADS] [CrossRef] [Google Scholar]
  7. Chandler, C. J., Carlstrom, J. E., & Scoville, N. Z. 1995, ApJ, 446, 793 [NASA ADS] [CrossRef] [Google Scholar]
  8. Drew, J. E., Proga, D., & Stone, J. M. 1998, MNRAS, 296, L6 [NASA ADS] [CrossRef] [Google Scholar]
  9. Eisenhauer, F., Perrin, G., Brandner, W., et al. 2011, The Messenger, 143, 16 [NASA ADS] [Google Scholar]
  10. Eisner, J. A. 2007, Nature, 447, 562 [NASA ADS] [CrossRef] [Google Scholar]
  11. Ercolano, B., & Rosotti, G. 2015, MNRAS, 450, 3008 [Google Scholar]
  12. Fairlamb, J. R., Oudmaijer, R. D., Mendigutía, I., Ilee, J. D., & van den Ancker, M. E. 2015, MNRAS, 453, 976 [Google Scholar]
  13. Frost, A. J., Oudmaijer, R. D., Lumsden, S. L., & de Wit, W. J. 2021, ApJ, 920, 48 [CrossRef] [Google Scholar]
  14. Garcia Lopez, R., Tambovtseva, L. V., Schertl, D., et al. 2015, A&A, 576, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. GRAVITY Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. GRAVITY Collaboration (Koutoulaki, M., et al.) 2021, A&A, 645, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. GRAVITY Collaboration (Caratti o Garatti, A., et al.) 2020a, A&A, 635, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. GRAVITY Collaboration (Garcia Lopez, R., et al.) 2020b, Nature, 584, 547 [Google Scholar]
  19. Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135 [Google Scholar]
  20. Hosokawa, T., Yorke, H. W., & Omukai, K. 2010, ApJ, 721, 478 [Google Scholar]
  21. Ilee, J. D., Wheelwright, H. E., Oudmaijer, R. D., et al. 2013, MNRAS, 429, 2960 [Google Scholar]
  22. Ilee, J. D., Oudmaijer, R. D., Wheelwright, H. E., & Pomohaci, R. 2018, MNRAS, 477, 3360 [CrossRef] [Google Scholar]
  23. Isella, A., & Natta, A. 2005, A&A, 438, 899 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Kama, M., Min, M., & Dominik, C. 2009, A&A, 506, 1199 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Klassen, M., Pudritz, R. E., Kuiper, R., Peters, T., & Banerjee, R. 2016, ApJ, 823, 28 [NASA ADS] [CrossRef] [Google Scholar]
  26. Koumpia, E., de Wit, W. J., Oudmaijer, R. D., et al. 2021, A&A, 654, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Koumpia, E., Koutoulaki, M., de Wit, W. J., et al. 2023, MNRAS, 519, L51 [Google Scholar]
  28. Kraus, S., Preibisch, T., & Ohnaka, K. 2008, ApJ, 676, 490 [NASA ADS] [CrossRef] [Google Scholar]
  29. Kraus, S., Hofmann, K.-H., Menten, K. M., et al. 2010, Nature, 466, 339 [Google Scholar]
  30. Kuiper, R., & Hosokawa, T. 2018, A&A, 616, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Lachaume, R. 2003, A&A, 400, 795 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Le Bouquin, J. B., Absil, O., Benisty, M., et al. 2009, A&A, 498, L41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Lumsden, S. L., Hoare, M. G., Urquhart, J. S., et al. 2013, ApJS, 208, 11 [Google Scholar]
  34. Martins, F., Schaerer, D., & Hillier, D. J. 2005, A&A, 436, 1049 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. McClure, M. K., D’Alessio, P., Calvet, N., et al. 2013, ApJ, 775, 114 [NASA ADS] [CrossRef] [Google Scholar]
  36. Mérand, A. 2022, SPIE Conf. Ser., 12183, 121831N [Google Scholar]
  37. Monnier, J. D. 2003, Rep. Prog. Phys., 66, 789 [NASA ADS] [CrossRef] [Google Scholar]
  38. Mottram, J. C., Hoare, M. G., Urquhart, J. S., et al. 2011, A&A, 525, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 412, 13 [Google Scholar]
  40. Pomohaci, R., Oudmaijer, R. D., Lumsden, S. L., Hoare, M. G., & Mendigutía, I. 2017, MNRAS, 472, 3624 [Google Scholar]
  41. Rosen, A. L., Li, P. S., Zhang, Q., & Burkhart, B. 2019, ApJ, 887, 108 [NASA ADS] [CrossRef] [Google Scholar]
  42. Rothman, L. S., Gordon, I. E., Barbe, A., et al. 2009, J. Quant. Spec. Radiat. Transf., 110, 533 [NASA ADS] [CrossRef] [Google Scholar]
  43. Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Stecklum, B., Caratti o Garatti, A., & Linz, H. 2012, ASP Conf. Ser., 464, 369 [NASA ADS] [Google Scholar]
  45. Tatulli, E., Malbet, F., Ménard, F., et al. 2008, A&A, 489, 1151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Torrelles, J. M., Patel, N. A., Curiel, S., et al. 2011, MNRAS, 410, 627 [NASA ADS] [CrossRef] [Google Scholar]
  47. Wheelwright, H. E., Oudmaijer, R. D., de Wit, W. J., et al. 2010, MNRAS, 408, 1840 [NASA ADS] [CrossRef] [Google Scholar]
  48. Zhu, Z., Stone, J. M., & Calvet, N. 2024, MNRAS, 528, 2883 [NASA ADS] [CrossRef] [Google Scholar]

1

DIT: 30 s; NDIT (on science): 12; Total Exposure Time of the SCI- CAL sequence: 36 minutes.

2

PMOIRED is a Python3 module and can be downloaded at https://github.com/amerand/PMOIRED

3

Due to the limited wavelength range able to be fit with the model, the uncertainties on these parameters were in some cases undefined. See Ilee et al. (2013), their Section 3.3 for further details.

All Tables

Table 1

Coordinates, bolometric luminosity, K-band magnitude, and distance of G033.3891 as taken from the RMS database. The mass and effective temperature of the source are also listed.

Table 2

Parameters of the best fit to the interferometric observables of the 2.2µm, Brγ, and CO emission towards G033.3891.

All Figures

thumbnail Fig. 1

uv-coverage of G033.3891 obtained using GRAVITY/VLTI on UTs. The position angle (in degrees) and projected baseline length (in metres) of each baseline pair are indicated in the legend.

In the text
thumbnail Fig. 2

Normalised spectrum for G033.3891 as observed with GRAVITY. The telluric corrected spectrum shows the Brγ and the CO bandheads in emission.

In the text
thumbnail Fig. 3

Model image of G033.3891 of the 2.2 µm continuum emission (thick ring) combined with the Brγ (compact Gaussian) and CO bandhead molecular emission (thin bright ring). The modelled image of the CO emission is the average of the individual images of the CO bandheads extracted at their central wavelengths (2.29, 2.32, 2.35, and 2.38 µm). The total geometric model (red line) fits well all interferometric observables (closure phases, T3PHI; differential phases, DPHI; visibilities, V; normalised flux, NFLUX) of the continuum, Brγ, and CO bandheads. PMOIRED normalises the flux from both the data and the model using the same algorithm. The grey shades represent the systematic errors of the interferometric observables, taking into account the calibration and transfer function uncertainties. The colour scheme of the different baselines follows Figure 1.

In the text
thumbnail Fig. 4

VLT/CRIRES observation of the CO bandhead in black (R ~ 30 000) overlaid with the best fitting disc model in red (adapted from Ilee et al. 2013).

In the text
thumbnail Fig. 5

Brγ photocentre shift displacements for different velocity channels calculated from the differential phases after correcting for the continuum contribution. The black contour corresponds to the inner ring size of the continuum emission.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.