© The Authors 2023

1 Introduction

The impact massive stars (M_⋆ ≥ 8 M_⊙) have on their host galaxies is disproportionate to their numbers: they provide strong mechanical and radiative feedback, and they deposit new elements in their surroundings through powerful stellar winds and supernova ejecta which constitute the building blocks of planets and life. Understanding the way in which these stars form is therefore important, yet challenging. The formation events are rare, take place deep inside a dust obscured natal environment, and unfold on such short timescales that the stars arrive on the main sequence prior to the dispersal of their natal cloud. Notwithstanding these complications, both observational (e.g., Johnston et al. 2015; Ilee et al. 2016; Cesaroni et al. 2017; Beuther et al. 2018) and theoretical (e.g., Krumholz et al. 2009; Kuiper et al. 2010; Rosen et al. 2016; Meyer et al. 2018) evidence is mounting that disk-mediated accretion is crucial to the process.

Recent studies that reveal disks, large-scale disk-like structures, and/or outflows, probe the embedded stages of formation (e.g., Frost et al. 2019, 2021a; Maud et al. 2018; Pomohaci et al. 2017). However, much remains unclear about how massive stars arrive on the main sequence, what their stellar properties are when they do, and how and when their accretion is halted. One of the ways to study the latest stages of formation is to observe newly formed massive stars in their natal environment. By making a census of the young populations in massive star-forming regions and determining their stellar and circumstellar (e.g., multiplicity, disks) characteristics, one can constrain possible formation scenarios and provide robust initial conditions for population synthesis studies (e.g., Bik et al. 2006). Through one such study in the very young (≤1 Myr) Galactic cluster M17, Ramírez-Tannus et al. (2017, hereafter RT17) made the unique discovery of six pre-main-sequence (PMS) stars in the mass range ~6–12 M_⊙ that have observable photospheres and, at the same time, show spectroscopic and photometric features typical for more embedded massive young stellar objects (MYSOs).

RT17 obtained optical to near-infrared (NIR) spectra with VLT/X-shooter of 12 young OB stars with masses ranging between ~6 and 25 M_⊙. For most of these objects, the photo-spheric part of the spectrum was of sufficient quality to allow quantitative spectroscopic modeling to determine stellar properties. Combining this with available photometry, they confirmed the PMS nature of six objects (~6–12 M_⊙) in the sample, based on their position in the Hertzsprung-Russell diagram (HRD; see Sect. 5.2 and Fig. 6). Two more objects (B289 and B215) were found to be close to, but not on, the zero age main sequence, displaying only mid-infrared (MIR) excess and lacking the other disk signatures listed below. These objects were characterized by RT17 as young stars whose PMS nature could not be solidly confirmed. The objects that were confirmed to be PMS stars possess features that point to the presence of circumstellar disks: (1) double-peaked atomic emission lines, most notably HI, Ca II, and O I; (2) CO-bandhead emission; and (3) NIR to MIR excess. Of these, especially CO bandhead emission has been associated with high accretion rates and the earlier stages of formation in young stellar objects (YSOs; see Sect. 5.2 and Table 9 for references).

The ro-vibrational transitions of the CO molecule that result in CO bandhead emission have recurrently been used to study the innermost regions of circumstellar disks in YSOs of all masses. Especially for the higher-mass YSOs, where these regions can rarely be resolved, they constitute one of the few commonly used diagnostics. A few studies focus on a single object (e.g., GRAVITY Collaboration 2020; Fedriani et al. 2020), revealing the spatial distribution of the emission with interferometry or integral field units. Other studies have gathered samples of intermediate to high-mass YSOs to study the disk properties in a more statistical manner (Bik & Thi 2004; Wheelwright et al. 2010; Ilee et al. 2013, 2014). In all these cases, the emission has been modeled either as a flat disk or a ring of material in Keplerian rotation, with the general result that the emission likely originates from hot (2000–5000 K) and dense (N_CO ~ 10²⁰–10²² cm⁻²) disk regions inside of the dust sublimation radius.

This paper’s focus is to constrain the properties of the cir-cumstellar material for the five objects among the RT17 sample of PMS stars in M17 that show CO bandhead emission. As these objects also reveal NIR to MIR excess, we include this thermal emission in our analysis. Based on these features, we aim to determine properties of the inner gaseous disk (<1 AU) where the CO emission likely originates, and the somewhat further regions (≥3–5 AU) where the majority of the thermal NIR and MIR dust emission is produced. Our sample consists of relatively high-mass (6–12 M_⊙) PMS objects that additionally have well-constrained stellar properties. This allows us to investigate whether the disk properties can be linked to stellar properties, such as mass and evolutionary stage.

Though we took a similar modeling approach to Bik & Thi (2004), Wheelwright et al. (2010), Ilee et al. (2013, 2014), we refined it in several ways (see Sect. 3.4 for details). First, for the first time, we took the observed second overtone emission into account (Δυ = 3, from 1.54 µm), as opposed to only the first overtone (Δυ = 2, from 2.29 µm). All the included (first and second overtone) bandheads were fit together, rather than individually. Second, by studying the dust emission with the available photometry we obtained robust continuum estimates, which we used to normalize the modeled bandheads. Third, we took ¹³CO into account, which increased the quality of our fits and allowed us to probe the isotopologue ratio ¹³CO/¹²CO as a further test of the young nature of our sample.

The paper is organized as follows. In Sect. 2 we present the data and the previous analysis performed by RT17. The disk model and the treatment of gas and dust are explained in Sect. 3. In this section we also explore the sensitivity of the predicted line spectrum to model parameters. In Sect. 4 we describe our fitting approach and the resulting disk parameters. Section 5 places these results in the context of previous literature and of the evolutionary state of the M17 PMS stars. Section 6 contains a brief summary and our main findings.

2 Data and previous analysis

The data used in this work consist of near- to mid-infrared photometry (1–11 µm) and the NIR part (1–2.4 µm) of spectra obtained with VLT/X-shooter.

2.1 Spectra

The X-shooter spectra (Vernet et al. 2011) of the 5 young stars studied in this paper were first analyzed by RT17, who derived stellar parameters from quantitative spectroscopy and optical (λ ≲ 1 µm) spectral energy distribution (SED) fitting. For a detailed listing of the spectroscopic observations used in this study we refer to that paper. Only for B275 we used a different spectrum, observed on 2019-06-06 for Program 0103.D-0099 (P.I.: Ramírez-Tannus). In Table 1 we summarize the results from RT17 that we use directly in our analysis.

For the fits of the CO bandheads we only used the NIR part (1–2.4 µm, with spectral resolution R = 13 000) of the spectra, which contain both the first and the second overtone emission. We (re-)reduced these spectra using version 3.3.5 of the X-shooter pipeline (Modigliani et al. 2010), running under the ESO Reflex environment (Freudling et al. 2013) version 2.11.0. The standard settings of the pipeline resulted in flux- and wavelength-calibration issues at the location of the first overtone bandheads. At wavelengths beyond ~2.29 µm, the flux showed wave-like fluctuations for some objects, hindering normalization. Additionally, in this range, wavelengths were shifted by ~80 km s⁻¹. In our attempt to solve these problems, while at the same time obtaining an optimal telluric correction, we first reduced the objects without flux calibration and then used the MOLECFIT tool (Smette et al. 2015; Kausch et al. 2015), version 1.5.9, for a first telluric correction. Importantly, MOLECFIT also provided a correction of the faulty wavelength calibration. We then divided the object spectra by their respective telluric standard stars that were reduced and MOLECFIT corrected in the same way. This way both the instrument response function and the bulk of the residuals of the MOLECFIT telluric correction are divided out. The continuum of the resulting spectra is the ratio of the continua of the science object and the telluric standard star. From this, it is in principle possible to obtain a flux calibrated spectrum by multiplying by the SED of the telluric standard star and applying some correction for slit losses. However, due to uncertainties involved in that procedure, we deemed a more reliable and consistent estimate of the continuum could be obtained from fitting SEDs to the photometric data points (Sect. 2.2). Therefore, we chose to normalize the spectra and compare the observations to models that are normalized to continuum estimates obtained from fitting SEDs, as described in Sect. 3.2.

Finally, we correct for the hydrogen absorption features of the telluric standard stars in the wavelength range 1.55–1.8 µm where they interfere with the CO second overtone bandhead line fluxes. We do not correct for the stellar and disk (emission) features of the science spectra in the same wavelength range, but exclude the most contaminated regions from our fits¹.

2.2 Photometry

The photometric points for our objects were taken from online catalogs and literature. We used photometry from the DENIS (Fouqué et al. 2000), 2MASS (Skrutskie et al. 2006), Spitzer GLIMPSE (Reach et al. 2005), and WISE (Wright et al. 2010; Jarrett et al. 2011) point source catalogs. For B275 and B331 we could add the extra points around 10 µm available from Nielbock et al. (2001) and Kassis et al. (2002). Though a few points beyond these wavelengths were available from the same papers, we did not include them in the final fits because they likely contain envelope emission, which was not considered in our models. All adopted photometric data per object are listed in Appendix A.

3 Methods: Analytical disk model with dust and gas

3.1 Keplerian disk model

The disk model that was used to fit the CO bandheads and continuum constitutes a flat disk with a Keplerian rotation profile, containing dust and/or CO gas². The disk has a radial temperature and surface density profile for the dust and the gas, described by analytical power laws as follows:

$$T\left(r\right)=T_i\cdot {\left(\frac{r}{R_i}\right)}^p$$(1)

$$N_{\text{H}}\left(r\right)={\left(N_{\text{H}}\right)}_i\cdot {\left(\frac{r}{R_i}\right)}^q$$(2)

where T is the temperature in K, N_H the hydrogen (H₂) surface density in cm⁻², and T_i, (N_H)_i their initial values at the inner radius R_i. The hydrogen surface density is converted to a ¹²CO gas density using the canonical CO/H₂ abundance of 10⁻⁴ (Lacy et al. 1994), and to a dust particle density using a gas to dust mass ratio of 100 (Bohlin et al. 1978). For including ¹³CO in the models, several abundances were tested (see Sect. 3.4). For the final results a standard interstellar abundance ratio N(¹²CO)/N(¹³CO) of 89 was used, taken from the HITRAN database³, which reproduces the isotopologue lines well. To avoid confusion, all densities are reported as hydrogen column density N_H.

Independent treatment of gas and dust

In all models that were used for fitting, the gas and dust properties are treated independently of each other, that is each have their own values. The reason for this approach stems mainly from the fact that the CO emission originates from hot gas well within the dust sublimation radius, with temperatures well above the assumed dust sublimation temperature of 1500 K. We did test a fitting approach in which dust and gas were combined and coupled with a source function ${\text{S}}_{\lambda }^{\text{tot}}$ as follows:

$${\text{S}}_{\lambda }^{\text{tot}}=\frac{{\tau }_{\lambda }^{\text{CO}}*{\text{S}}_{\lambda }^{\text{CO}}+{\tau }_{\lambda }^{\text{dust}}*{\text{S}}_{\lambda }^{\text{dust}}}{{\tau }_{\lambda }^{\text{tot}}}$$(3)

with the total optical depth ${\tau }_{\lambda }^{\text{tot}}={\tau }_{\lambda }^{\text{CO}}+{\tau }_{\lambda }^{\text{dust}}$ the sum of the CO line and dust optical depths, the CO and dust source functions given by Planck curves according to the same temperature (so in this case ${\text{S}}_{\lambda }^{\text{CO}}={\text{S}}_{\lambda }^{\text{dust}}$) and where all quantities are a function of wavelength. This revealed that the CO emission and the dust emission responsible for the NIR-MIR excess are representative of different regions in the disk, with little or negligible overlap. Where the two do overlap the main effect of importance is the optical thickness of the dust “damping” the CO emission. However, since dust only exists at lower temperatures, where the CO emission is already significantly reduced, this effect is minimal.

Of course, even though the CO line and dust continuum emission likely do not originate from the same region, the dust continuum still has a significant influence on the strength of the observed bandheads: for the same absolute CO emission, a stronger continuum leads to a weaker bandhead signal relative to the continuum. Therefore, it remains important to determine and model the continuum. We opted to fit the SEDs and the CO bandheads separately with a dust-only and gas-dust disk respectively; where for the latter the dust parameters were set by the best fits obtained from the SEDs. For the fitting of the SEDs not only dust continuum, but also stellar continuum was taken into account. We now elaborate on the specifics of the respective modeling and fitting approaches for the continuum and CO bandhead emission.

3.2 Continuum treatment

The continuum flux has a stellar component and one that arises from thermal emission of dust near the star. The dust contribution is computed using the described analytic disk model including only dust. For the dust 0.1 µm astronomical silicate was used with opacities from Laor & Draine (1993). The stellar and extinction parameters as determined previously by RT17 (T_eff, log ɡ, and A_V; Table 1) are input in our calculations using the corresponding Kurucz models (Kurucz 1993; Castelli & Kurucz 2003) and reddening (Cardelli et al. 1989) the resulting SED after adding the disk. The near- and mid-infrared photometry points are then fit by means of a grid of disk models – stellar parameters remain fixed.

The outer radius of all the dust disk models was fixed to 500 AU, regardless of the temperature at this point. We do not expect to constrain an outer radius of the entire disk, because the photometric points available for our objects are only representative of the warm to hot dust (T ≳ 150 K) and emission from cooler dust does not change the models in the fitted wavelength ranges (1–12 µm).

For the inner radius R_i we used a radiative equilibrium approach (following Lamers & Cassinelli 1999) to calculate the radius at which a certain dust temperature is reached as a function of the previously determined T_eff and R_⋆ (Table 1). This R_i is therefore dependent on T_i for each object. This approach results in five free parameters: T_i·, p, (N_H)_i·, q from Eqs. (1), and the inclination i. The grid of parameter values used for fitting is given in Table 2. Only for B331 a good fit could not be obtained with R_i as function of T_i and stellar parameters. We therefore fitted the SED of this object using a separate grid in which also R_i is a free parameter and whose values are given in Table 3.

3.3 CO emission treatment

3.3.1 Combining dust and gas

In the models used to compute the normalized CO bandheads, the dust and CO gas are both included, but do not share the same temperatures and densities, as motivated above. Instead, the dust parameters were fixed by the parameters obtained from the SED fits, and the CO gas parameters were varied over in a grid of models. Only the inclination was shared between the dust and gas disk.

The outer radius of the gas disk model is taken to be the point where the temperature drops to 600 K. Below this temperature the gas does not substantially contribute to the bandhead emission anymore, because the excited vibrational states are not sufficiently populated. The outer radius for the dust disk is again 500 AU, as before (see Sect. 3.2).

Because the CO gas excitation temperatures are typically much higher than the dust temperatures and because the inner radii for the gas-disk are so close to the star (see Table 4) the gas and the dust do not necessarily overlap. In fact, they overlap only in the models where the CO gas temperature drops to 600 K beyond the dust sublimation radius (calculated as described in Sect. 3.2). In these cases the source function in the overlap region was calculated according to Eq. (3), where the CO and dust source functions are given by Planck curves according to their respective temperatures. In all regions and cases where there is no overlap, that is only gas and no dust or only dust and no gas, the gas and dust emission are calculated separately and added. For the CO modeling the main significance of the dust disk is its contribution to the continuum, to which the bandheads are normalized.

3.3.2 The gas disk and parameter grid

The vibrational-rotational level populations of CO are assumed to be in LTE, with Einstein A coefficients and level transitions taken from Li et al. (2015). The line profiles of the rotational transitions are assumed to be Gaussian with a width υ_G. This parameter accounts for intrinsic width of the lines as well as thermal and micro-turbulent velocities. The intensities are calculated per ring at radius r, with a dependence on the ring temperature and density. They are then convolved with a Keplerian velocity profile per ring and corrected for inclination, before integrating over the radial extent of the disk. Finally, the model spectra are convolved with the instrument resolution of R = 11 300 and normalized to the continuum fits presented in Sect. 4.1, adjusted only to match the inclination of the CO disk⁴.

As the masses of the stars are previously determined (Table 1; needed for Keplerian velocities), we have 7 free parameters for our model: T_i·, p, (N_H)_i·, q, R_i from Eqs. (1) and (2), the inclination i, and the Gaussian width υ_G. To fit the data we calculated a grid of models with parameter ranges given in Table 4. The parameter values for this final grid are not all regularly spaced and are based both on literature and on several other exploratory grids that were calculated to investigate the relevant parameter space for our objects. The inner gaseous disk radius R_i is given in units of stellar radii (Table 1, RT17). Both the initial hydrogen column density (N_H)_i and the R_i values are evenly spaced on a log scale. The positive value (+1) for the density exponent q was included, because of evidence for an inverse density profile of the innermost regions of the disk (Antonellini et al. 2020).

3.3.3 Correcting for pseudo-continuum

When calculating the final CO bandhead spectrum, each subsequent bandhead flux is added to the P-branch (ΔJ = −1) and higher R-branch (ΔJ = +1) rotational transitions of the previous one. This, when convolved with the instrument resolution, creates a kind of secondary continuum and influences the overall appearance of the bandhead spectrum. This is true for both overtones (see also Fig. 1). However, when normalizing the observed spectra in the second overtone wavelength regions it is impossible to recognize such effects, especially because the overlapping line wings of the hydrogen Brackett series also influences the continuum there. To account for this pseudo-continuum effect, we introduced a correction of the bandhead continuum of the modeled second overtone bandheads before fitting. We did this by fitting a 3rd degree spline through the minima just before the bandhead-onset wavelengths and dividing the modeled band-heads (already normalized to star and dust continuum) by this “bandhead continuum”. The normalization of the data in the first overtone region is much less problematic and the models there did not require similar corrections.

3.4 Model exploration

Several examples of CO bandhead fitting using analytical disk models are available in literature (e.g., Kraus et al. 2000; Bik & Thi 2004; Ilee et al. 2013). Similar to our approach, these works adopt an LTE treatment of the gas, Keplerian rotation and a 1D thin disk approximation. Here we list changes and improvements in our approach relative to these earlier efforts.

First, we introduce an internally consistent treatment of the continuum and normalization. In previous work, the bandheads are usually continuum subtracted and normalized to the first bandhead. Instead, we estimate the continuum from SED fitting and use this to normalize our models. With this approach the information contained in the strength of the bandheads is retained.

Second, we fit the first five bandheads of the first vibrational overtone (v = 2−0, 3−1, 4−2, 5−3, 6−4) at the same time. We also predict the sixth bandhead (υ = 7–5), but because it is at the very edge of the observed spectral range, the signal-to-noise ratio is insufficient to include it in the fits.

Third, for the first time, the second vibrational overtone bandheads (Δυ = 3, from 1.5 µm onward) are included.

Fourth, we include ¹³CO. These lines are weak and their diagnostic value in constraining disk properties is limited. However, including these lines allows for the isotopologue ratio¹³ CO/¹²CO to be constrained.

Finally, we probe higher temperatures than the commonly assumed dissociation temperature of 5000 K (Bosman et al. 2019); grid values run up to 8000 K. Such high temperatures should be taken with care as they may not be physical, since CO dissociation is not included in the models. Allowing for higher temperatures may mimic limitations in the LTE assumption, which can overestimate the temperature if noncollisional excitation mechanisms, such as UV pumping, contribute to setting the state of the gas. The maximum temperatures retrieved from our fits can then serve as an additional probe of the validity of the LTE assumption.

To facilitate a comparison with previous results of similar fitting efforts in literature we identify where and how these novel aspects in our approach might impact the parameter estimation. To do this we explore the effects of the free parameters on the appearance of the first and second overtone bandhead emission using Fig. 1. This figure was made by over-plotting several models, sequentially varying only one parameter and keeping all others fixed at a value fitting to the B275 bandheads. The fitted model is plotted in black in each panel to give an impression of where the data for this object are. In the very last panel we show the effects of including ¹³CO in varying abundances; we discuss the inclusion of ¹³CO at the very end of this section.

The amount of emission in both the first and second overtone bandheads is most sensitive to the temperature and the extent of the (hot) gaseous disk, as can been seen on the T_i, p, and R_i plots in Fig. 1. These parameters have very similar effects on the amount of flux of the bandheads. Since the temperature always drops going outward, a small initial radius will lead to a small area of the hottest gas, that is less flux. This can be compensated for by a shallower temperature drop (smaller |p|) resulting in a larger emitting surface. We note the strong emission in the T = 6000 K model, which is stronger than the commonly observed strength of these features. Again, CO dissociation is not included in our models.

Since the disk is assumed to be Keplerian, the velocities are set by the mass of the central star (kept constant for each object), the inner radius R_i and the inclination i⁵. These parameters, therefore, determine the onset wavelengths of the bandheads and influence their shape. High (more edge on) inclinations and smaller inner radii broaden the bandhead profile and shift the onset wavelength blue-ward. Both these parameters, however, also influence the total flux.

For the relative strength of the first and second overtone bandheads optical depth effects play an important role. The lines of the second overtone are intrinsically weaker and therefore remain optically thin for higher column densities. For this reason line broadening, that is higher υ_G, leads to higher flux only for the first overtone bandheads and not for those of the second overtone. The inclination has slightly less effect on the second overtone flux than on the first overtone flux: the radiating surface of an optically thick inclined disk is reduced, but the second overtone remains optically thin for longer sight-lines through the disk.

Finally, increasing the initial column density (N_H)_i, or decreasing |q|, beyond certain grid values has a larger impact on the second overtone flux, because these lines remain optically thin within the grid values. These changes to the second overtone flux are more detectable than the changes to the first overtone within crucial ranges of (N_H)_i. Therefore, including the second overtones in the fits especially helps toward constraining the column density. We now discuss the effects of each of the four modeling improvements.

Fig. 1 Overplotted models, showing for each parameter the effect of varying its value, while keeping all other parameters fixed at: Ti = 4500 K; p = −0.75; (NH)i = 8.3 × 1025 cm−2; q = −0.5; Ri = 1.6 R⋆; i = 30º; υG = 1 km s−1. The models are all without inclusion of13 CO, except in the bottom panel, where the effect of adding 13CO in varying abundances is shown. The continuum and stellar parameters used for these models are those of B275. The “fixed” model is colored black on each plot and is a fit to the B275 spectrum.

Fig. 2 Observed first overtone CO bandheads of B163 (black line) shown together with models varying in 12CO/13CO ratio – a decreasing ratio signifying more 13CO. The red line indicates the best fit model for this target, with no 13CO included, the green lines indicate models differing only in the amount of 13CO. Including 13CO with the ISM value (89) clearly improves the fit with respect to fitting 12CO only. Decreasing the ratio leads to progressively worse fits, indicating that 13CO is likely not enhanced in the circumstellar environment of this object.

3.4.1 Continuum normalization

It is important to keep in mind that the stellar parameters and thermal dust continuum that are kept fixed for each star in the fitting process, do strongly influence the appearance of the band-heads - and therefore the fit results. A thorough discussion of these effects is given in Appendix B. Here we mention the most important effect: a higher continuum originating from the star and dust disk will lead to weaker bandheads for the same gas disk model.

The most important constraint from including the bandhead strength relative to the continuum is on the ranges of values for p and υ_G. Shallow temperature gradients and high intrinsic line velocities lead to very high bandhead fluxes beyond observed ranges, that cannot be compensated for by other parameters. Thus, we limited the probed ranges in our grid accordingly (Table 4).

3.4.2 Fitting multiple first overtone bandheads

The effect described earlier (in Sect. 3.3) that the bandhead series, if strong enough, creates a secondary continuum, can be clearly seen on Fig. 1 for any parameter that influences the strength of the bandheads. The later bandheads (i.e., those with higher vibrational quantum numbers) are also more sensitive to temperature changes, likely due to excitation effects. These effects generally lead to better temperature constraints when fitting more bandheads together.

3.4.3 Including the second overtone bandheads

Following from our previous observations, the most important determinant of the detectability and (relative) strength of second overtone bandheads is the column density (profile) of the disk.

We tested this by comparing fit results with or without including second overtones. The exponent q is poorly constrained in general and even its selective effect on the second overtone did not change fit results significantly. However, including the second overtone did yield significant changes for the initial column density (N_H)_i.

3.4.4 Including ¹³CO

CO bandhead emission is a signature not unique to YSOs -it is evidence for the presence of a (hot and dense gaseous) disk, but as such it is also observed to originate in the (decre-tion) disks of evolved B[e] stars (Kraus et al. 2000). In these stars the stellar surface layers can be enriched in ¹³C, due to the evolution of the (massive) star and chemical mixing processes. Stellar winds then cause this anomalous isotopologue ratio to be carried into the circumstellar environment, where consequently the N(¹²CO)/N(¹³CO) abundance ratio drops and ¹³CO bandhead signals are enhanced. Kraus et al. (2020) show that therefore, including ¹³CO emission in the models, can provide a test to distinguish the evolved objects from the younger ones.

The objects in our sample are expected to be YSOs due to their young cluster environment. To assess the sensitivity of our results to the ¹³CO feature, we computed models for gradually increasing ¹³CO abundances. In the last panel of Fig. 1 we show the model with only ¹²CO and with ratios of ¹²CO/¹³CO down to 4 (the ~ 22-fold enhancement of the abundance found by Kraus et al. (2020) in a B[e] supergiant). Even for the standard (interstellar) ratio 89, ¹³CO leaves a detectable signal in the first overtone, which becomes very prominent for the highest abundances. In Fig. 2 we show the observed first overtone bandheads of B163, the object for which the ¹³CO feature is most pronounced, together with its best fit model and models with varying ¹³CO abundance. When included with the standard abundance, the ¹³CO feature clearly improves the fit; but increasing the abundance leads to progressively poorer results, limiting a possible enhancement of the abundance to ~30% in this case. We also find that including ¹³CO does not change the other disk parameters that are fitted. These findings are similar in the other objects; save for B331, where the hydrogen Pfund line series interferes with the ¹³CO bandhead signal and we cannot assess the impact of the isotopologue. We conclude that, although the inclusion of ¹³CO is of limited diagnostic value to constrain disk parameters, its detected features, consistent with an interstellar abundance, do provide a confirmation of the young nature of our objects. Therefore, when presenting and plotting our results we use the models including ¹³CO, except for B331.

Fig. 3 Best fit model SEDs and photometric data points. On each figure the data and results for B275 (in black) are plotted along with one other object (in blue) for comparison. Kurucz models representing the stellar spectrum are given in dashed lines; the combined stellar and disk model in solid lines. The scattered points mark the photometry. Also indicated are the wavelength ranges of the second and first overtone bandheads. Note how B331 is both more luminous and also lacks NIR excess, with the dust emission starting around 4–5 µm.

4 Results

4.1 Continuum fits

The aim of the SED fits is, first, to determine a reliable continuum for normalizing the modeled CO bandheads, and second, to obtain constraints on the inner parts of the dust disk.

Figure 3 visualizes the best fitting model SEDs for the five objects in our study, overplotted with the available photometric points. The case of B331 is different and we discuss it separately.

For the objects B163, B243, B268, and B275 the available photometry could be fitted well with optically thin dust emission. One of the consequences of this is that the inclinations remain unconstrained. We used a reduced χ² statistic, in addition to an inspection by eye, to determine the range over which a reasonable fit can be obtained for each parameter. These ranges, along with the best fits are reported in Table 5. The temperature T_i and column density (N_H)_i at the inner rim are reasonably well constrained, with T_i ≈ 1500 K and (N_H)_i ≈ 10²¹ cm⁻² for all four objects. The temperature and density exponents p and q are degenerate: a less extended emission region (higher |q|) can be compensated for with a slower temperature decline (lower |p|). The derived values of the exponents are different than for the gaseous disk (Sect. 4.2), with the values for p generally shallower and the values for q steeper than those obtained by fitting the CO bandheads. Apart from degeneracies, this may also reflect the difference between the disk regions probed, that is the gaseous part close to the star (see Sect. 4.2) versus the inner parts of the dust disk.

With respect to the best fit model, B163 has an excess in the J (~1.2 µm) and the H (~1.7 µm) bands. This could be because the disk is not actually optically thin (see also Sect. 5.4) or due to refractory dust grains with sublimation temperatures in excess of 1500 K. The normalization of the second overtone bandheads is only slightly affected; maximally by a factor of 1.3.

The SED of B331 could not be fit with the grid used for the other objects. It is the only object which does not have a NIR dust excess within the X-shooter spectral range (the excess begins at λ ~ 3.6 µm), signifying that the star is surrounded by cooler dust. Here too, the parameter that is best constrained is the dust temperature at the inner rim. However, when using the radiative equilibrium approach to calculate R_i the cool dust disk would start as far as ~140 AU, leading to poor fits. We opt to fit B331 with R_i as a free parameter and arrive at a dust disk with R_i ≈ 30 AU. For this object the modeled emission is optically thick, hence the best fit column density should be considered a lower limit. The inclination is again poorly constrained, this time because it is degenerate with the density decline |q|, that is the extent of the disk.

We use these best fits and their parameters when modeling the dust for the CO bandhead models and for normalizing those models, as described in Sect. 3.3. It is more challenging to use the results to truly constrain the parameters of the inner dust disk. Ultimately, because the data are simply too scarce to determine all the parameters, we find it more realistic to describe the results in terms of an amount (in the optically thin case) or surface (in the optically thick case) of dust in a ring around the star at a certain temperature. Because our photometry points only reach up to 12 µm, they are only representative of dust emission from a limited part of the dust disk, that is the hotter, inner parts. To understand and compare the characteristics of these parts we define the inner dust disk as the part of the disk where 99% of the dust emission between 1 and 12 µm comes from. Using the best fitting model we can then provide an inner and outer radius and a surface area for this emission region, and determine the temperatures and column densities at those radii. For B331 we calculate the inclined disk surface, because of the emission being optically thick. Integrating the column density between inner and outer radii provides a total mass for the emission region. The results are reported in Table 6.

The inner and outer rim column densities as well as the masses of the defined emission region are very similar for the first three objects, even though their outer temperatures (T_out) and their areas are rather different. In combination with identical inner rim temperatures (T_i), it follows that the best fitting models are not sensitive to the size of the emission region and its temperature structure, and that it is thus the hottest dust that dominates the emission. B275 clearly has more dust emission, which is reflected in a higher disk mass.

The dust emission from B331 is dominated by longer wavelength emission from cooler dust. The NIR emission seen in the other objects is absent (Fig. 3). The fact that the inner rim temperature is so low suggests a dust-free inner cavity in the disk. A higher density and a slower density decline lead to a much higher disk mass than for the other objects. The caveat here is, however, that it is possible that the three longest wavelength points are contaminated with envelope emission. The even longer 20 µm and 37 µm points from Lim et al. (2020) could not be fit with our disk model, most likely due to the presence of envelope emission which we did not include in our models.

4.2 CO bandhead fits

The best fits resulting from the grid of models described in Sect. 3.3, are plotted for each object in Fig. 4, together with the data and the fit regions. For B331 we fitted only two of the first overtone bandheads, due to strong emission in the hydrogen Pfund line series. All other spectra are fitted including five first overtone bandheads. The second overtone bandheads suffer from poor signal to noise due to a combination of hydrogen and other emission features and telluric contamination. We fitted only those bandheads that were more or less free from hydrogen emission. Even if they were not clearly detected, at least some second overtone bandhead regions were always included in each fit, as their absence also provides constraints.

Because of the degeneracy between different parameters and limits to the data quality and spectral resolution, it was nontrivial to understand and determine the accuracy of fits to the data. The most straightforward approach of quoting the minimal χ² best fit values was insufficient as “next-best” fits sometimes produced very different parameter values. In addition, especially the power law exponents p and q were poorly constrained within the probed ranges. To get better insight into the parameter space and to determine errors, we used marginalized likelihood distributions for each parameter, as described in detail in Appendix C. One caveat with this approach is that, in general, the errors and mean depend on the parameter space chosen. In this case, the parameter space is large and limited to physically reasonable values, so that we can assume that likelihoods for parameter values outside the grid are (close to) zero. To better understand the temperature, density and velocity information, we fixed the p, q, and υ_G parameters at their maximum likelihood values after fitting with the entire grid. We determined the best fits for the remaining parameters (T_i, (N_H)_i, R_i, i) from a “reduced” grid with the p, q, and υ_G parameters fixed.

The results are summarized in Table 7. For T_i, (N_H)_i, R_i, and i we determined both the most probable and the mean values, with their respective uncertainties, based on the marginalized probability distributions. As explained in Appendix C, the best fit, the most probable and the mean values can be rather different, but for almost all parameters the best fit values fall close to and within the error bars of either the mean or the most probable value or both - we list the one that is closest to the best fit value. In the table the most probable value is indicated as “max prob”, the mean value as “mean” and the best fit as “bf”. Only for the inner radius R_i of B243 does the best fit value in the reduced grid lie outside the error bars of the quoted mean. This is the object with the weakest bandheads and the poorest signal to noise.

In order to illustrate and compare these results in Table 7 a bit better we made a similar emission region analysis as for the inner dust disk (Table 6). In Table 8 we show the inner and outer radii of the part of the disk where, according to the best fitting model, the CO emission originates. By construction the outer radius R_out is given by the point at which the gas temperature drops to 600 K. We also show where the dust disk begins (R_dust) and provide the column densities at R_i and R_dust. Finally, we provide the total gas mass of the CO emitting region (since there is no dust there), based on the assumed CO abundance.

Fig. 4 Best fits (in red) overplotted with data (in black) for all objects. The data points included in the fits are marked in orange. Note the different y-axis scales for the first (right) and second (left) overtone. All fluxes are normalized. On the second overtone bandhead plots the hydrogen Brackett series (not marked) can also be seen, in absorption from the stellar photosphere and/or in emission from the disk.

4.3 Summary of the most important aspects of the results

Before we discuss our results we briefly bring together the different aspects that can help form a picture of the studied disks. Firstly, the SEDs are fit with a thin dusty disk, but the parameters remain poorly constrained apart from the temperature at the inner rim. Except for B331 all inner rim temperatures are at or very close to the assumed dust sublimation temperature, suggesting that these inner dust disks are undisrupted. B331 is an exception with cooler dust at larger radii, pointing to a perturbed inner dust disk (see Sect. 5.2 for further discussion). We do not further discuss the SED-derived densities for reasons explained in Sect. 5.4.

From Tables 7 and 8 we see that in all cases the CO band-heads originate in a relatively narrow ring close to the star. There does not appear to be a relation between R_i and T_i. Temperatures range between 2000 and 5000 K − the lower limit given by our grid, the upper likely being related to the temperatures at which CO abundances become negligible due to hot gas chemistry dissociating the CO molecules (Bosman et al. 2019), the local radiation field temperature being too low to cause photo-dissociation. Thus, despite probing higher temperatures, we retrieve the expected maxima in the best fit values, confirming that even if non-LTE mechanisms are at work, it is unlikely that they would drastically change our results. Taking into account the continuum and fitting multiple bandheads are important measures to ensure that the emission region and temperatures are reliably constrained.

The distribution of best-fit inclinations is consistent with random orientations. The Gaussian line widths are on the order of the thermal velocities at the CO temperatures, suggesting that turbulent or other stochastic motions are not significant in the emitting regions.

Though the column densities are similar for all objects, we see that they are slightly lower (by a factor of ~3−8) where the second overtone is not observed (see objects B243, B331 in Fig. 4). This is consistent with the second overtone being observed only when column densities are high enough. The fact that the emission originates from a narrow ring means that the density change over larger parts of the disk and hence q will be poorly constrained (also apparent on Fig. 1). Nonetheless, the most likely value for q was found to be identical for all objects. The last column in Table 8 shows that the total masses of the CO emitting region are surprisingly similar. This largely relates to, but is not fully explained by, the similarity in the derived densities.

In general we note that, given the accessible parameter space, the overall similarity in the emitting regions is striking. To illustrate this we show in Fig. 5 the range of temperatures and column densities that are in principle consistent with the normalized strength of the observed bandheads. Especially the range in column densities spans orders of magnitude, while the fit results are similar within one order of magnitude. We return to this observation in the discussion.

Fig. 5 Illustration of the temperature and column density ranges that, for fixed parameters Ri ≈ 0.25 AU, q = −1.5, i = 40°, and υG = 1 km s−1, lead to CO first overtone bandhead maximum flux values between 2 and 25% continuum, which represent the minimum detectable and the maximum observed normalized fluxes respectively. Left: the accessible parameter space for different values of the temperature exponent p (models are for B268). Right: the same, for the different stars as varying in mass and continuum emission (models are for p = −2).

5 Discussion

5.1 Comparison with previous CO bandhead modeling results

We first compare the inner gas disk parameters determined from the CO bandhead emission with those resulting from similar modeling efforts in literature - particularly those that pertain to YSOs in the higher mass ranges.

There are mainly two approaches when modeling CO bandhead emission toward YSOs. The first one assumes a thin disk model with free parameters similar to ours (e.g., Kraus et al. 2000; Wheelwright et al. 2010; Ilee et al. 2013, 2014). The second approach differs from the first mainly in that it does not fit a temperature or density structure, but models an isothermal ring at one density and v sin i, where the velocity can be converted to a distance assuming an inclination and Keplerian rotation (e.g., Bik & Thi 2004; Koutoulaki et al. 2019; Fedriani et al. 2020; GRAVITY Collaboration 2020). Both approaches have been used to successfully fit first overtone emission. We do note that all previous studies normalize the bandheads to the peak of the first bandhead (except Kraus et al. 2000; Wheelwright et al. 2010), thereby ignoring the strength of the emission. Furthermore, in our study, for the first time, the second overtone CO bandhead emission is reported and fitted, leading to better constraints on the column density. With these differences in approach in mind, we summarize the similarities and differences in the results.

Our derived CO column densities are all of order N_CO ∼ 10²¹ cm s⁻², which is in excellent agreement with literature values that are mostly within the range N_CO ~ 10²⁰−10²² cms⁻². There is, however, often a larger spread in column densities when a sample of objects is studied (e.g., Ilee et al. 2013, 2014). While the morphology and strength of our bandhead profiles is equally diverse as in these studies, the determined densities within our sample are remarkably similar.

The spread in temperature values is similar to those found in previous studies. All inner radii in our sample are well within 1 AU, which is on the lower end of values derived for this parameter by others. Ilee et al. (2013, 2014), for instance, find inner radii up to 6 AU, but mostly around 1−2 AU. They also generally find rather shallow temperature exponents (p 0.6 on average), leading to (very) large outer radii (determined by where the temperature drops below 1000 K) ranging from a few to hundreds of AU. In our study, the strengths of the bandheads are very sensitive to R_i and p (Fig. 1), leading to (mostly) steep temperature exponents that mimic the effect of modeling a ring of emission close to the star. This justifies the assumptions underlying ring models. In the only study to date where the CO bandhead emission of a MYSO⁶ is spatially resolved, the authors indeed find, by fitting the visibilities for each bandhead separately, that the CO emitting region must be relatively small ∆R/R ≤ 20%, where $R={0.58}_{-0.04}^{+0.04}$ AU is the ring radius (GRAVITY Collaboration 2020). This radial location is in good agreement with our findings and the ring thickness is comparable to and even slightly narrower than our “rings” of emission (Table 8).

Similarly to other studies confronting the disk model to a sample of systems (Ilee et al. 2013, 2014), the density exponent q remains poorly constrained. This is not surprising since the emission region is narrow. This parameter is therefore also of little consequence to other results or statements about the disk (e.g., the total mass).

For the inclinations we find a spread in values, consistent with random orientations. Ilee et al. (2014) find a preference for higher (closer to edge-on) inclinations and suggest that this may be due to a geometric selection effect. This is interesting in the context of detection rates of CO bandhead emission (see Sect. 5.2). Within our small sample, we do not find evidence for such a selection effect, in agreement with the larger sample of Ilee et al. (2013). In fact, we predict (see Fig. 1) that higher inclinations make the (optically thick) CO emission rather less detectable, due to a decrease in effective emitting surface.

Caratti o Garatti et al. (2017) and Fedriani et al. (2020) report CO emission toward two deeply embedded MYSOs to be the reflection of light from the inner gaseous disk onto the perpendicular outflow cavity wall. In these cases the disk is almost fully edge on and the central source with its inner gas disk is too enshrouded for the CO emission to be detected directly. However, because the outflow cavity “sees” the inner disk face-on, the observed CO emission is rather fitted with low velocities (or a face-on orientation). It is highly unlikely that the CO emission toward our objects is indirect considering that all of them have detected photospheres. Thus, we do not expect that our inclinations are “artificially” low. But these observations can potentially help explain why detection rates of CO bandhead emission remain low overall: when looking at (embedded) MYSOs at high(er) inclinations one might not observe the innermost disk regions.

Finally, the intrinsic line width v_G (also denoted as ∆v in literature) is often taken to be a measure of turbulent motions of the gas. In our case, the main constraint on this parameter comes from line fluxes of the optically thick first overtone, that become very high for higher widths (>5 km s⁻¹ ) that exceed thermal velocities of the gas at the given temperatures (Sect. 3.4 and Fig. 1). Studies that fit the bandheads normalized to the first bandhead peak lack the sensitivity to this effect of intrinsic line widths. Such studies commonly yield higher values up to 10–20 km s⁻¹.

In summary, our results point to very similar conditions for the origin of the CO emission as previous studies and, if anything, place even stronger constraints. Despite the fact that our objects represent different masses and (likely) different evolutionary stages (see Sect. 5.3) than (M)YSOs investigated in the studies mentioned so far, it seems that the emission originates in disks that at least locally look very much alike. Even toward low-mass YSOs very similar conditions are derived, where the CO bandhead emission is known to signal a temperature inversion, that is a heated disk atmosphere by stellar irradiation, winds or magnetospheric accretion heating (e.g., Najita et al. 2007, see also Sect. 5.2). Thus, it is plausible to assume that in many if not all YSOs the emission originates in disk regions with similar conditions, if not (locally) similar disks.

5.2 CO bandhead detection rates and the link with accretion and age

CO first overtone bandhead emission (in this section referred to as CO emission) has been observed in YSOs over a large mass range, as well as in Be and B[e] stars, which can be YSOs, but also objects of a more evolved nature such as B[e] supergiants (Lamers et al. 1998; Kraus et al. 2020) or classical Be stars (Cochetti et al. 2021). Here we focus on findings for YSOs.

We summarize the detection rates among YSOs of different masses and formation stages in Table 9. Detection rates, especially among low-mass YSOs, appear to be higher for earlier formation stages. In general, CO emission toward lower luminosity YSOs is associated with properties that point to the earliest stages of formation and/or high accretion rates, for instance deep embedding, energetic outflows, veiling and/or measured high accretion rates (Najita et al. 2007). In the intermediate to highmass YSOs detection rates are consistent with this trend. The lower rates among Herbig Ae/Be stars can perhaps be understood in terms of these objects representing a later stage in formation than most MYSOs, yet earlier than T Tauri’s (see also Sect. 5.3).

Apart from high occurrence rates in objects undergoing an accretion burst (e.g., Contreras Peña et al. 2017), the link with accretion is supported by correlations of CO emission (strength) with other features that signal accretion, such as high veiling (Connelley & Greene 2010) and Br-γ luminosity (Ilee et al. 2014; Pomohaci et al. 2017). Caratti o Garatti et al. (2017) report on an accretion burst in a MYSO of ~20 M_⊙ where the CO emission is detected only during the burst. Despite all this evidence for a link with accretion, actual reported accretion rates span a large range of values: from ~10⁻⁸ M_⊙ yr⁻¹ in the lowest mass objects (Koutoulaki et al. 2019) to ~5 × 10⁻³ M_⊙ yr⁻¹ in the highest mass objects (Caratti o Garatti et al. 2017), with Herbig Ae/Be stars in between with ∼10⁻⁷ −10⁻⁶ M_⊙ yr⁻¹ (Ilee et al. 2014; Contreras Peña et al. 2017). On the other hand, a theoretical study by Ilee et al. (2018b) predicts CO emission to be most prominently observable with moderately high accretion rates of ∼10⁻⁴−10⁻⁵ M_⊙ yr⁻¹.

Finally, it should be noted that CO emission detection rates remain relatively low over all mass ranges, also in accreting objects. In MYSOs it is definitely less common than other emission features, such as Br-γ emission at 2.16 µm, also in the K-band, and which is detected toward 70–90% of objects in the quoted studies. Moreover, CO emission is not thought to be a direct probe of accretion, that is it does not originate in the accretion flow itself, where it is likely destroyed by UV light from the surface. As such it seems reasonable to follow Koutoulaki et al. (2019) in their conclusion that accretion seems to be a necessary but not sufficient condition for CO emission to be observed. The low detection rates could be due to geometric effects (as discussed in Sect. 5.1), veiling due to high continuum emission or simply varying conditions in disks around accreting objects. It remains remarkable, however, that such vastly different accretion rates around objects with a wide range in mass could lead to emission that consistently points to very similar disk conditions (see Sect. 5.1).

Our own sample was taken from a total of six YSOs, or eight, if we include the two young stars with only MIR excess from RT17 (B289 and B215 on Fig. 6; see also Sect. 1). It is surprising to find that five (63–83%) of these objects show CO bandhead emission. This raises the question as to the evolutionary state of these objects, which we address in the following section.

5.3 Evolutionary stage(s) of the intermediate to high-mass YSOs in M17

The mass range of the YSOs in our sample (6− 12 M_⊙) spans the upper mass ranges of Herbig Be stars and the lower-mass ranges of MYSOs. To better understand the formation stage of the studied sample, we briefly discuss the most important characteristics of the two categories and make comparisons to objects in each category.

Massive YSOs are minimally understood to be objects that are in the process of forming a star (or more in a multiple system) that is at least 8 Mo when reaching the zero age main sequence (ZAMS). The widely used MYSO catalog Red MSX survey (RMS; Lumsden et al. 2013) selects objects based on their MIR brightness and high luminosities (≳10⁴ L_⊙), which are hence defining properties for studies in the near- and mid-infrared (e.g., Ilee et al. 2013; Frost et al. 2019). Some of these studies explicitly exclude from their definition objects that have started to ionize their surroundings to form an HII region (e.g., Wheelwright et al. 2010; Pomohaci et al. 2017). However, this specification is not included in studies in (sub)millimeter to radio wavelengths (e.g., Beltrán & de Wit 2016; Maud et al. 2018; Johnston et al. 2015; Ilee et al. 2018a), that reveal large scale (∼10²−10⁴ AU) molecular outflows and/or disk-like structures. Takami et al. (2012) find that these mm-bright objects can be embedded at MIR wavelengths. Studies that explicitly include both wavelength regimes are rare (Hsieh et al. 2021), which makes it hard to compare objects in either categories. It appears that the term MYSO is used for objects with a range of masses and evolutionary stages, from the near molecular core phase up to objects for which (some) photospheric features are observed, covering thereby a wealth of objects that differ in their characteristics. However, we note that the overwhelming majority of studies into MYSOs pertains to objects that are highly embedded in gas-dust envelopes (A_V ≳ 20–30) and have (estimated) masses of at least 10–15 M_≳ with no detectable photospheric features.

Of the M17 PMS stars, B331 with 12 M_≳ best fits the category MYSO, while also having a well detected photosphere. In a detailed study of a MYSO of ∼25 M_⊙, Frost et al. (2019) find evidence for a 60 AU inner clearing of the dust disk, while also reporting on the detection of CO bandhead emission. This combination of features resembles our findings for B331. Notably the authors refer to transition disks in low-mass young stars which show dust clearing while still having a small gaseous disk where accretion could be ongoing (Wyatt et al. 2015). They tentatively suggest that the object of their study could represent such a stage in a MYSO and propose photo-evaporation or a companion as the most likely dust clearing mechanisms. These conclusions are corroborated in their follow-up study of a sample of eight MYSOs, from which they find evidence for an evolutionary sequence, the later stages of which are indeed characterized by inner dust clearing by photo-evaporation (Frost et al. 2021a,b). According to this interpretation the disk of B331 is in a state of transition. For our other objects this suggests an earlier stage of formation where they have not yet started to clear their inner dust disk. This scenario is consistent with the masses and luminosities of the respective objects. The two young stars from RT17 close to the ZAMS, with masses of 20 M_⊙ and 10 M_⊙ (B289 and B215 on Fig. 6), also support this picture as they lack all inner (gaseous or dust) disk signatures, but do show MIR excess similar to B331.

Recent studies that identify and characterize large samples of Herbig Ae/Be stars describe these objects as optically revealed PMS stars in the mass ranges 2−10 M_⊙, characterized by emission lines and NIR excess attributed to disks (Vioque et al. 2018, 2022; Guzmán-Díaz et al. 2021). Because the M17 PMS stars fit this description well, we highlight the recent study of Vioque et al. (2022; hereafter V22) for comparison. V22 analyze a sample of 128 Herbig Ae/Be stars with stellar masses up to ~20 M_⊙. Based on optical spectra, Gaia parallaxes and photometry they derive stellar parameters and accretion rates.

The V22 Herbig stars with masses >4 M_⊙, most of which are spectral type B, are shown on the HRD in Fig. 6 together with the M17 PMS stars. In the figure we include B337, the one PMS star studied by RT17 without CO bandhead emission, but with otherwise very similar characteristics. We also include the two objects from that study (B289 and B215; in brown), that show no signatures of inner (gaseous or dust) disks in their X-shooter spectra. The more extincted M17 sources (A_V ~ 13) are shown in purple. From the HRD plot, it appears that the M17 PMS objects are all younger than the V22 sources. This is consistent with the derived extinctions for the respective samples; with A_V = 2.0−6.7, and mean $\overline{A_{\text{V}}}$ for the V22 sub-sample (M_*, > 4 M_⊙), and A_V ~ 7−13 for the M17 sources (RT17 and Table 1). Furthermore, V22 report photo-evaporation of inner dust disks for sources above 7 M_⊙ based on a decline in NIR/MIR excess, similar to B289 and B215. In this context, the youth of the M17 PMS sources is again confirmed by the presence of inner dust and gaseous disks, with the most luminous and highest mass object (B331) showing signs of inner dust disk clearing (in line with the comparison to the MYSO from Frost et al. 2019).

V22 do not mention CO bandhead emission, as the necessary wavelength ranges were not included in their study. The main reference for the occurrence rate of CO bandhead emission amongst Herbig Ae/Be stars therefore remains Ilee et al. (2014). Of the 23 objects in their sample with masses ≳4 M_⊙, 13−17% exhibit CO bandhead emission (masses from Fairlamb et al. 2015). This is on the low side compared to detections toward MYSOs, which, as proposed earlier, could be because these objects are typically in a later stage of formation than most MYSOs. Though we did not perform a full census of PMS stars in M17, the presence of CO bandhead emission in five out of six them is yet another pointer toward a young age.

Fig. 6 Hertzsprung-Russell diagram displaying the MIST pre-main-sequence tracks (Dotter 2016) for stellar masses between 6 and 20 M⊙, with the solid black line on the left denoting the ZAMS and the gray dashed line on the right indicating the birthline. The objects in this study, along with the one M17 PMS star from RT17 without CO band-head emission (B337), are shown together with the >4 M⊙ Herbig stars from the Vioque et al. (2022) sample. Also included are the two RT17 objects (B289 and B215) close to the ZAMS without inner dust or gaseous disk signatures, but with MIR excess. The M17 PMS stars appear to represent a younger stage of formation, consistent with the high detection rate of CO bandhead emission in their spectra and their high extinction. For comparison, the estimated CO detection rate representative for >4 M⊙ Herbig stars is 13–17% (see text).

5.4 Limitations and open questions

In this final section we discuss the limitations of this work and make suggestions for further research.

The column densities derived from the SED fitting (Table 5) are orders of magnitude lower than those derived from the CO bandhead fitting (Table 7). The explanation for this lies in the fact that the disk model used lacks vertical disk structure. This means that the relatively high dust temperatures associated with the NIR emission are likely valid only on the disk surface, which can hide a cold and dense mid-plane underneath. Thus the derived column densities from the dust are representative only of the top disk layers and it is highly unlikely that the thermal emission is truly optically thin throughout the vertical disk extent as the best fitting models suggest.

Similar statements can be made for the CO bandhead modeling: we have derived the properties of the region in the disk where the emission originates, which, as we have concluded, is relatively small. Therefore, it is difficult to make statements about the disk as a whole. Longer wavelength data (MIR and (sub)mm) are needed to constrain the masses and sizes of these disks. This, in turn could give further information on the disk evolution around these PMS stars - for instance whether the disks are (also) eroded outside-inward.

The question whether and how much these stars are accreting is also an open one. V22 estimate accretion rates of ~10⁻⁴ −10⁻⁶ M_⊙ yr⁻¹ for objects in the mass range 6−12 M_⊙, based on line luminosities and extrapolations of magnetospheric accretion models. Since the M17 PMS stars appear younger than the V22 sample and show CO bandhead emission, they could well be still accreting with rates on the higher end. However, since the similarity in disk properties from CO emission seems hard to reconcile with the disparity in reported accretion rates, the link with accretion remains ambiguous.

With regard to the possible accretion mechanisms, that is boundary layer (BL) vs. magnetospheric accretion, we note that although CO emission probes regions close to the star, for most of our objects it still originates from outside the coro-tation radius. With v sin i values taken from RT17 and using Eq. (2) from Ilee et al. (2014) we find corotation radii of ≤0.081,0.19 and 0.14 AU for B243, B268, and B275 respectively. Though we have no constraint on v sin i values for B163 and B331, the CO emission in those objects also (mostly) originates outside these radii. Thus, only for B268 (and perhaps B163) the CO emitting region covers the corotation radius. This could be an indication that the disk extends close enough to the stellar surface for BL accretion to take place, but for most of our objects the CO emission cannot probe this, in agreement with Ilee et al. (2014).

6 Summary and conclusions

In this work we have studied the CO bandhead and thermal infrared emission from the disks around five intermediate to massive PMS stars in M17, with the aim to constrain their inner dust and gaseous disk properties and further our understanding of the end stages of massive star formation. To this end, we developed an LTE disk model that accounts for dust and CO overtone emission. For the first time, we fit the second overtone. We fit normalized bandheads, using models for stellar and dust continuum. We compared and discussed the derived disk properties in the context of the previously derived PMS evolutionary states of the central stars (RT17) and CO bandhead emission detection (rates) among YSOs in the literature. We arrived at the following conclusions:

1. Taking into account the breadth of the parameter space and the diversity in line morphology, we have arrived at surprisingly similar inner gaseous disk characteristics for the CO emitting region, which is consistent with results in the literature for YSOs of different masses. Therefore, it seems that this emission is typical of certain (more or less) specific physical conditions, that is to say high densities (N_CO ~ 10²¹ cm⁻²) and temperatures (2000–5000 K).

2. The (overall low) detection rates reported in the literature suggest that CO overtone bandhead emission is preferentially observed toward objects with high(er) accretion rates and in earlier stages of their formation, across all mass ranges. The high detection rate among the M17 PMS stars, in combination with their position on the HR diagram, suggests that these objects are exceptionally young for their mass range.

3. The previous conclusion is supported by the SED analysis, which points to undisrupted inner dust disks, except in the case of the most luminous object B331, for which the disk is likely in transition.

4. Though it is likely that these objects are accreting, the CO emission is not suited to probe the rate and mechanism of accretion. Results in this work point to CO bandhead emission regions close to the star, but outside the corotation radius (save for one object) within which boundary-layer versus magnetospheric accretion can be probed (in agreement with Ilee et al. 2014). Analysis of hydrogen lines for these objects is more suited to probe accretion regions, as these lines originate even closer to the stellar surface (Backs et al. 2023).

5. Taking into account the strength of the bandheads as well as the inclusion of the second overtone emission provides more stringent constraints on the emission region and column densities than previous studies lacking these diagnostics.

6. We clearly detect the ¹³CO feature, but find no evidence for its enhancement in the spectra, consistent with the PMS nature of the M17 objects. Though of low diagnostic value to constrain disk parameters, we corroborate that ¹³CO abundances can be useful in determining the evolutionary state of objects with CO overtone emission, as suggested by Kraus et al. (2020).

7. The M17 intermediate to high-mass YSOs appear unique in their combination of having both observable photospheres while also spectral and SED properties typical of younger, more enshrouded objects. Multiwavelength analysis applied to the same objects is crucial for a better understanding of these objects in particular and the different stages of massive star formation in general.

Our analysis has demonstrated that the studied pre-main-sequence stars in M17 are relatively young, and that they are excellent candidates for studies at longer wavelengths, for example with ALMA, to further constrain the final stages of the formation process of intermediate to high-mass stars.

Acknowledgements

We express our gratitude to the anonymous referee for their helpful comments and insights toward improving this manuscript. J. P. acknowledges support from the Dutch Research Council (NWO)-FAPESP grant for Advanced Instrumentation (P.I.: L. Kaper). This work is based on observations collected at the European Organization for Astronomical Research in the Southern Hemisphere under ESO programs 089.C-0874(A) and 103.D-0099. We thank SURF (www.surf.nl) for support in using the Lisa Compute Cluster. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. This work depended for a major part on the use of the Python programming language, in particular the packages NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), and Matplotlib (Hunter 2007).

Appendix A Photometry tables

Tables A.1 to A.3 list the photometry points used for the fits.

Appendix B Continuum normalization and impact of stellar parameters

The effects stellar and continuum characteristics have on the appearance of the CO bandhead emission are visualized in Figure B.1. Here we overplot the bandheads resulting from the exact same disk model for the 5 different objects in our sample. There are three main effects on the model fluxes.

Firstly, the continuum to which the bandheads are normalized affects their appearence. This can been seen in Figure B.1 when comparing the first overtones of B243, B163, B268, and B275. The first three objects have similar SED fluxes in the relevant wavelength range (see Figure 3), while B275 has a much stronger continuum and therefore the weakest bandheads for the same model. In the same way, the relative strength of the continuum in the first and second overtone wavelength regions also influences the relative strength of the respective set of band-heads. For B331 the slope of the SED is much steeper than for the other objects, causing the second overtone for this object to be relatively weak.

Secondly, the stellar radius is the measuring unit for R_i, which is the main reason that the first overtone of B331 is so much stronger: its stellar radius is simply much larger than that of the other objects (Table 1), leading to a larger emitting surface for the same parameters. It also influences the Keplerian velocities which decrease as the square root of the stellar radius.

Finally, the Keplerian velocities are proportional to the square root of the stellar mass. This effect, however, is largely compensated for by the previous point, as stars with a higher mass generally also have larger radii.

Fig. B.1: Comparison of a model with the same disk parameters, changing only the object, that is the continuum, stellar mass, and stellar radius. The model in black is a good fit to the B268 spectrum. Note the difference in scale for the first and second overtone. The used disk parameters are: Ti = 5000 K; p = −0.75; (NH)i = 8.3 × 1025 cm−2; q = −1.5; Ri = 1.1 R*, ; i = 50°; vG = 2 km s−1 ; and 13CO is not included.

Appendix C Fits and error calculations

In order to get insight into the parameter space and to determine errors on the fit results, we made marginalized likelihood distributions for each parameter as follows. We converted all reduced χ² values to a likelihood $L=e^{-{\chi }_{red}^2/2}$ and established the likelihood of a value a for a parameter p_i by summing over the likelihood of all other parameter value combinations:

$$L\left(p_i=a\right)={\displaystyle \sum _{B_i}L\left(a;B_i\right)}$$(C.1)

where B_i is the set of all possible outcomes for all parameters except p_i . This likelihood was then normalized by the area under the curve (see Figure C.1 for visualization):

$$P_{p_i}\left(a\right)=\frac{1}{{\displaystyle {\int }_{-\infty }^{\infty }{L}_{p_i}}\left({p^{\prime }}_i\right)d{p^{\prime }}_i}{L}_{p_i}\left(a\right)$$(C.2)

where the integration runs over the grid value range. The obtained 1D probability distributions P_pi have a mean:

$$\text{μ=}{\displaystyle {\int }_{-\infty }^{\infty }{p^{\prime }}_i{P}_{p_i}}\left({p^{\prime }}_i\right)d{p^{\prime }}_i$$(C.3)

and standard deviation:

$${\sigma }^2={\displaystyle {\int }_{-\infty }^{\infty }{\left({p^{\prime }}_i-\text{μ}\right)}^2{P}_{p_i}\left({p^{\prime }}_i\right)}d{p^{\prime }}_i$$(C.4)

where the integrations run over grid value range. Only if the resulting distribution is Gaussian, the best-fitting value and the maximum likelihood (or most probable) value are equal to the mean, which for most of our fits and parameters is not the case. An example of the resulting distributions is shown on the diagonal of Figure C.1. The standard deviations on the most probable values (indicated as “max”) are calculated as in Equation (C.4), replacing µ with the maximum likelihood value. The confidence intervals indicated as “conf” are the 1-σ confidence intervals obtained by integrating between µ − σ and µ + σ. All errors were calculated such that no off-grid values fall within the 1-σ confidence interval. Because of this and because the resulting distributions are not Gaussian, the 1-σ confidence intervals does not generally result in (the Gaussian) 0.68. Most intervals have values between 0.55 and 0.75. On the off-diagonal panels in the figure the 2D equivalents of the described probability distributions are shown for each pair of parameters. The probabilities are color-coded with lighter colors indicating higher probability for a pair of values. These plots give insight into degeneracies in the parameter space and show how well a parameter is constrained.

The value that is quoted as the final result in Table 7 in this work, is the most probable or the mean value, depending which was closest to the best fit.

Fig. C.1: Off the diagonal: 2D marginalized likelihood distributions (see text), for the model grid fitting results on fitting five first and six second overtone bandheads for B275. On the diagonal: the 1D probability distribution functions (PDF) for each parameter resulting from normalizing the marginalized likelihood distributions to the area under the plotted line. The meaning of the mean, max, and conf values is explained in the text. The patterns in the 2D distributions give insight into the degeneracies between parameters.

References

All Tables

All Figures

Fig. 1 Overplotted models, showing for each parameter the effect of varying its value, while keeping all other parameters fixed at: Ti = 4500 K; p = −0.75; (NH)i = 8.3 × 1025 cm−2; q = −0.5; Ri = 1.6 R⋆; i = 30º; υG = 1 km s−1. The models are all without inclusion of13 CO, except in the bottom panel, where the effect of adding 13CO in varying abundances is shown. The continuum and stellar parameters used for these models are those of B275. The “fixed” model is colored black on each plot and is a fit to the B275 spectrum.

In the text

Fig. 2 Observed first overtone CO bandheads of B163 (black line) shown together with models varying in 12CO/13CO ratio – a decreasing ratio signifying more 13CO. The red line indicates the best fit model for this target, with no 13CO included, the green lines indicate models differing only in the amount of 13CO. Including 13CO with the ISM value (89) clearly improves the fit with respect to fitting 12CO only. Decreasing the ratio leads to progressively worse fits, indicating that 13CO is likely not enhanced in the circumstellar environment of this object.

In the text

Fig. 3 Best fit model SEDs and photometric data points. On each figure the data and results for B275 (in black) are plotted along with one other object (in blue) for comparison. Kurucz models representing the stellar spectrum are given in dashed lines; the combined stellar and disk model in solid lines. The scattered points mark the photometry. Also indicated are the wavelength ranges of the second and first overtone bandheads. Note how B331 is both more luminous and also lacks NIR excess, with the dust emission starting around 4–5 µm.

In the text

Fig. 4 Best fits (in red) overplotted with data (in black) for all objects. The data points included in the fits are marked in orange. Note the different y-axis scales for the first (right) and second (left) overtone. All fluxes are normalized. On the second overtone bandhead plots the hydrogen Brackett series (not marked) can also be seen, in absorption from the stellar photosphere and/or in emission from the disk.

In the text

Fig. 5 Illustration of the temperature and column density ranges that, for fixed parameters Ri ≈ 0.25 AU, q = −1.5, i = 40°, and υG = 1 km s−1, lead to CO first overtone bandhead maximum flux values between 2 and 25% continuum, which represent the minimum detectable and the maximum observed normalized fluxes respectively. Left: the accessible parameter space for different values of the temperature exponent p (models are for B268). Right: the same, for the different stars as varying in mass and continuum emission (models are for p = −2).

In the text

Fig. 6 Hertzsprung-Russell diagram displaying the MIST pre-main-sequence tracks (Dotter 2016) for stellar masses between 6 and 20 M⊙, with the solid black line on the left denoting the ZAMS and the gray dashed line on the right indicating the birthline. The objects in this study, along with the one M17 PMS star from RT17 without CO band-head emission (B337), are shown together with the >4 M⊙ Herbig stars from the Vioque et al. (2022) sample. Also included are the two RT17 objects (B289 and B215) close to the ZAMS without inner dust or gaseous disk signatures, but with MIR excess. The M17 PMS stars appear to represent a younger stage of formation, consistent with the high detection rate of CO bandhead emission in their spectra and their high extinction. For comparison, the estimated CO detection rate representative for >4 M⊙ Herbig stars is 13–17% (see text).

In the text

Fig. B.1: Comparison of a model with the same disk parameters, changing only the object, that is the continuum, stellar mass, and stellar radius. The model in black is a good fit to the B268 spectrum. Note the difference in scale for the first and second overtone. The used disk parameters are: Ti = 5000 K; p = −0.75; (NH)i = 8.3 × 1025 cm−2; q = −1.5; Ri = 1.1 R*, ; i = 50°; vG = 2 km s−1 ; and 13CO is not included.

In the text

Fig. C.1: Off the diagonal: 2D marginalized likelihood distributions (see text), for the model grid fitting results on fitting five first and six second overtone bandheads for B275. On the diagonal: the 1D probability distribution functions (PDF) for each parameter resulting from normalizing the marginalized likelihood distributions to the area under the plotted line. The meaning of the mean, max, and conf values is explained in the text. The patterns in the 2D distributions give insight into the degeneracies between parameters.

In the text