K-band High-resolution Spectroscopy of Embedded High-mass Protostars

Figure 2 shows the Br γ line profiles for 11 sources. In Table 2 we list the equivalent widths of the emission/absorption lines, and by Gaussian fitting, the measured FWHMs. Line emission is observed toward eight sources over all the categories of H ii region associations, including "No H ii region" (IRAS 18151-1108 and S255 IR SMA2), HC H ii region (IRAS 04579+4703 and IRAS 23151+5912), and UC H ii region (AFGL 490, W33A-MM1 Main, G29.96-0.02, and NGC 7538 IRS1). All but G29.96-0.02 show a triangle-like profile with a full width half maximum (FWHM) of ∼70–400 kms⁻¹, while G29.96-0.02 shows a significantly narrower linewidth (FWHM ∼ 28 kms⁻¹). Out of the remaining three sources, W3 IRS 5 shows an absorption line with a FWHM of ∼37 kms⁻¹.

Zoom In Zoom Out Reset image size

Table 2. Properties of Brγ Emission/Absorption

Name	EW	ΔV	AV (JH)	LBrγ	Lacc	${\dot{M}}_{\mathrm{acc},\mathrm{model}}$	${\dot{M}}_{\mathrm{acc},\mathrm{cluster}}$
	(Å)	(kms−1)		(L⊙)	(L⊙)	(M⊙ yr−1)	(M⊙ yr−1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
W3 IRS 5	0.66 ± 0.03	36.6 ± 2.1	19.4 ± 7.2	⋯	⋯	⋯	⋯
AFGL 490	−3.39 ± 0.02	160.8 ± 1.6	27.9 ± 0.5	0.21 ± 0.01	847 ± 28	$\left(1.26\pm 0.04\right)\times {10}^{-5}$	$\left(1.33\pm 0.04\right)\times {10}^{-5}$
IRAS 04579+4703	−5.37 ± 0.05	67.9 ± 0.9	21.0 ± 1.7	0.011 ± 0.001	60 ± 7	$\left(8.2\pm 1.0\right)\times {10}^{-7}$	$\left(8.9\pm 1.1\right)\times {10}^{-7}$
S255IR SMA1	<0.06	⋯	24.3 ± 2.8	<0.0001	<0.8	< 8.7 × 10−9	< 9.7 × 10−9
S255IR SMA2	−1.119 ± 0.004	200.3 ± 7.2	24.3 ± 2.8	0.0019 ± 0.0004	11.7 ± 2.3	$\left(1.25\pm 0.25\right)\times {10}^{-7}$	$\left(1.39\pm 0.27\right)\times {10}^{-7}$
W33A-MM1 Main	−3.17 ± 0.05	180.2 ± 2.9	21.0 ± 0.9	0.066 ± 0.005	300 ± 20	$\left(2.36\pm 0.16\right)\times {10}^{-6}$	$\left(2.77\pm 0.18\right)\times {10}^{-6}$
IRAS 18151-1208	−3.43 ± 0.04	148.9 ± 2.0	23.3 ± 5.2	0.051 ± 0.020	$\left(2.3\pm 0.9\right)\times {10}^{2}$	$\left(2.5\pm 0.9\right)\times {10}^{-6}$	$\left(2.8\pm 1.0\right)\times {10}^{-6}$
G29.96-0.02	−89.76 ± 0.05	28.1 ± 0.1	26.8 ± 1.2	9.6 ± 0.9	− a	− a	− a
AFGL 2591 VLA3	<0.03	⋯	34.9 ± 5.1	<0.001	<9.1	< 9.7 × 10−8	< 1.1 × 10−7
NGC 7538 IRS1	−2.83 ± 0.05	80.9 ± 1.7	28.8 ± 7.2	0.1 ± 0.1	$\left(4.5\pm 2.3\right)\times {10}^{2}$	$\left(3.8\pm 1.9\right)\times {10}^{-6}$	$\left(4.3\pm 2.2\right)\times {10}^{-6}$
IRAS 23151+5912	−3.52 ± 0.11	404.6 ± 11.1	24.0 ± 7.2	0.3 ± 0.2	$\left(1.3\pm 0.6\right)\times {10}^{3}$	$\left(1.0\pm 0.5\right)\times {10}^{-5}$	$\left(1.2\pm 0.6\right)\times {10}^{-5}$

Note. Col. (2) Equivalent widths. Col. (3) Full width half-maximum velocity obtained from the Gaussian fitting. Col. (7) Disk accretion rate estimated using ${\dot{M}}_{\mathrm{acc}}=\tfrac{{L}_{\mathrm{acc}}{R}_{\star }}{{{GM}}_{\star }}$ with M_⋆ as a function of L_bol from the model in Davies et al. (2011). Col. (8) Disk accretion rate estimated using ${\dot{M}}_{\mathrm{acc}}=\tfrac{{L}_{\mathrm{acc}}{R}_{\star }}{{{GM}}_{\star }}$ with M_star = M_⋆,cluster

^a The emission is probably associated with a UC H ii region, which is not related to the parameters for mass accretion. See text for details.

Download table as:  ASCIITypeset image

The measured equivalent widths of the emission lines range from ∼−1.1 to −5.4 Å except for G29.96-0.02, which has an equivalent width larger by a factor of 20–80 (∼90 Å). For non-detected sources we estimated an upper limit of the equivalent width with a S/N < 3 assuming a line width of 50 km s⁻¹. The equivalent widths and the upper limits reported here are broadly consistent with those measured by previous observations at low resolutions (Appendix B).

Figure 3 shows the spectra at 2.288–2.334 μm, which covers the CO v = 2 − 0 and 3 − 1 bandhead emission. Both bandheads are seen in emission in W33A MM1 Main, and the v = 2 − 0 bandhead is marginally detected in IRAS 18151-1208.

Zoom In Zoom Out Reset image size

To measure the equivalent width of the v = 2 − 0 bandhead emission, we calculated the continuum baseline by employing a 3rd-order polynomial fit in the line-free regions of 2.190–2.291μm and 2.300–2.320μm and subtracted it. We then integrated the spectrum over 2.293–2.298μm. The equivalent widths measured for W33A MM1 Main and IRAS 18151-1208 are −2.33 ± 0.69 and −0.67 ± 0.55 Å, respectively. For the remaining sources, we derived a 3σ upper limit for the CO v = 2 − 0 emission of −1.0 to −4.4 Å. Near-IR spectroscopy at low-to-medium resolutions has already been done for most of the sources, and the equivalent widths and the upper limits we measured are broadly consistent with the previous observations (Appendix C).

We do not derive the equivalent width of the v = 3 − 1 bandhead emission observed toward W33A MM1 Main for the following reasons: (1) narrow v = 2 − 0 absorption lines cover the same wavelength range (see Figure 3 and Section 4.3); and (2) our spectra do not cover the continuum at the longer wavelengths, making the continuum fitting inaccurate.

In Figure 3 all sources except G29.96-0.02 show narrow CO v = 2 − 0 absorption lines at 2.32–2.33 μm. In Figure 4 we show the R(7)–R(12) profiles in detail. For these profiles the velocities are shown with respect to the source velocity measured using the millimeter observations. Table 3 shows their equivalent widths.

Zoom In Zoom Out Reset image size

To increase signal-to-noise ratios, we averaged the R(7)-R(10) profiles (i.e., those with the largest absorption) for all the sources but AFGL 490 and IRAS 04579+4703. These line profiles are shown at the bottom of each panel. For the latter two sources we use R(7) and R(8) only for the averaged spectra due to relatively low signal-to-noise of the upper transitions.

The absorption line profiles show a rich variety between the sources. Given the systemic velocity from sub/millimiter observations (Table 1), we found that, for S255IR SMA1 the absorption is deepest at ΔV ∼ 15 km s⁻¹, with blueshifted wing absorption toward ΔV ∼ − 15 km s⁻¹. S255IR SMA2 shows a similar line profile but with an absorption maximum at ΔV ∼ 20 km s⁻¹ and a blueshifted wing down to ΔV ∼ 0 km s⁻¹.

AFGL 2591 VLA3 shows two separated absorption components, with absorption maxima at ΔV ∼ − 30 and −5 km s⁻¹, respectively. We measure FWHM widths of 18 and 10 km s⁻¹, respectively, using Gaussian fitting. The line profiles observed toward W3 IRS 5 can be also attributed to two components with Gaussian profiles, with peak velocities of −13 and −3 km s⁻¹, respectively. The Gaussian fitting indicates that the individual absorption components are not clearly resolved at the given velocity resolution of 8.6 km s⁻¹.

The remaining sources show a single absorption component line with a Gaussian profile. The absorption maximum is blueshifted by 3–8 km s⁻¹ for AFGL 490, W33-MM1 Main, IRAS 18151-1208 and NGC 7538 IRS1; and redshifted by ∼2 km s⁻¹ for IRAS 23151+5912. The absorption component for IRAS 04579+4703 does not show any clear evidence for a blueshift or redshift. Those for AFGL 490, IRAS 18151-1208 and IRAS 23151+5912 have a Gaussian-fitted FWHM width of 8–12 km⁻¹. The remaining profiles are not resolved at the given velocity resolutions of 4.3 and 8.6 km s⁻¹.

We estimate extinction toward each source using the method adopted by Cooper et al. (2013) as described below. We calculate two extinctions; A_V (JH) and A_V (HK), using the JHK magnitude with the following equation:

$\begin{eqnarray}&&{A}_{V}=\displaystyle \frac{{m}_{1}-{m}_{2}+{c}_{\mathrm{int}}}{{0.55}^{1.75}({\lambda }_{1}^{-1.75}-{\lambda }_{2}^{-1.75})}\end{eqnarray} \tag{ 1 }$

where m₁ and m₂ are the magnitudes, λ₁ and λ₂ are the corresponding wavelengths, c_int is the intrinsic color (J − H = − 0.12 or H − K = − 0.05) of a B0 star from Koornneef (1983).

In Figure 5, we plot A_V (JH) versus A_V (HK) of our sources together with those of the Cooper et al. (2013) sample. The sources observed by Cooper et al. (2013) show approximately linear correlation between A_V (JH) and A_V (HK), but the latter is systematically higher by a factor of typically ∼1.5. In contrast, our sources have A_V (HK) significantly larger than A_V (JH), up to a factor of ∼4, due to the K-band excess shown in Section 2. Therefore, we will adopt A_V (JH) to correct extinction in Section 5.2. These values are tabulated in Table 2. We provide more complete information, including the 2MASS JHK magnitudes and A_V (HK), in Appendix A.

Zoom In Zoom Out Reset image size

One might be concerned that we cannot apply the above conventional extinction law, in particular if we observe the above fluxes via scattering on the disk or in an outflow cavity. Actually, for some HPMOs well-studied using millimeter interferometry, we do not seem to be able to directly observe light from the protostars in the infrared (e.g., Takami et al. 2012). In this case, the target emission must be extended. However, only a few of our targets were marginally resolved during our observations: the remaining targets are not resolved at the given seeing (06–08, corresponding to 700–3700 au). Observations at a higher angular resolutions (e.g., by using James Web Space Telescope) would allow us to better investigate this issue.

As shown in Section 4.1, the observed Br γ emission shows a triangular profile with a FWHM of 70–400 km s⁻¹, with most cases ranging from ∼100–200 kms⁻¹. These are similar to the Br γ profiles observed toward low-mass PMS stars (Folha & Emerson 2001). For low-mass PMS stars, these line luminosities are well correlated with the disk accretion rates (Muzerolle et al. 1998), and the emission is highly likely to be associated with magnetospheric accretion columns (Calvet et al. 2000).

One might think that the Br γ emission for intermediate protostars to HMPOs could have a different origin. Tambovtseva et al. (2016) presented modeled profiles of Br γ emission for Herbig AeBe stars associated with the following: (1) a magneto-centrifugal disk wind, (2) an accretion disk heated by the illumination of the central source (hereafter "a passive disk"), (3) a disk heated by accretion viscosity (hereafter "an active disk"), and (4) magnetospheric accretion. However, these line profiles cannot explain the observed profiles adequately, as described below. For three of the above four cases, a symmetric line profile with a single peak would require a very specific viewing angle: i ≪ 20° (nearly pole-on) for the passive disks and the disk winds; and ∼60° for magnetospheric accretion. These trends do not match the the following observations: among the sources with Br-γ emission, disk inclination angles of 30°–60° were measured for three of them (AFGL 490, W33A-MM1 Main and IRAS 18151-1208) using millimeter interferometry (Table 1).

The active disk models would, on the other hand, yield a symmetric line profile with a single peak over a large range of viewing angles (i < 40°). However, their FWHMs are up to ∼50 km s⁻¹, significantly smaller than the observations. Furthermore, these models would yield a triangular profile over a very small range of inclination angles (i ≪ 20°): the remaining line profiles resemble a Gaussian profile with blueshifted and redshifted wing components. This, again, is not consistent with the observations of the disk inclination angles described above.

Based on the line profiles similar to those of low-mass PMS stars, we assume that the Br γ emission is associated with magnetospheric accretion. Then, we estimate the accretion luminosities and disk accretion rates as follows. We first derive the Brγ luminosity using the observed equivalent width and the 2MASS K magnitude, correcting with extinction derived using the JH magnitudes (i.e., A_v (JH); see Section 5.1). The accretion luminosity is then estimated using the same method as Pomohaci et al. (2017). We first use the empirical relationship between L_Brγ and the accretion luminosity L_acc for Herbig Ae and Be stars (Mendigutía et al. 2011):

$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{L}_{\mathrm{acc}}}{{L}_{\odot }}\right)=(3.55\pm 0.80)+(0.91\pm 0.27)\mathrm{log}\left(\displaystyle \frac{{L}_{\mathrm{Br}\gamma }}{{L}_{\odot }}\right)\end{eqnarray} \tag{ 2 }$

Then the disk accretion rate ${\dot{M}}_{\mathrm{acc}}$ is obtained using the following equation:

$\begin{eqnarray}&&{\dot{M}}_{\mathrm{acc}}=\displaystyle \frac{{L}_{\mathrm{acc}}{R}_{\star }}{{{GM}}_{\star }}\end{eqnarray} \tag{ 3 }$

where R_⋆ is the stellar radius, M_⋆ is the stellar mass, and G is the gravitational constant. Here we adopt the stellar mass from the models of the zero-age-main-sequence-star in Davies et al. (2011) given the L_bol or observational estimates (Table 1). The stellar radius R_⋆ is then obtained from the same models. It is noteworthy that the stellar radius might increase by a factor of a few to ∼100 for a massive protostar at an accretion phase (Hosokawa et al. 2010). If this is the case, the mass accretion rate is underestimated, which might explain the mass accretion rate problem (see below).

Figure 6 shows that the disk accretion rate in comparison with a possible evolutionary scheme: i.e., No H ii region → HC H ii → UC H ii (Beuther & Shepherd 2005). No clear correlation is found in between the mass accretion rate and evolution. The derived disk accretion rates ranges between ≲ 10⁻⁸ and ∼10⁻⁴ M_⊙ yr⁻¹. Most of these accretion rates are not sufficient for forming massive O stars (M_* ≳ 10M_⊙) in their formation timescale of ∼10⁵ yr (McKee & Tan 2003). A possible explanation is that these stars have already accreted most of their mass. If this is the case, our selected sources, as bright in infrared, have completed their main-accretion phase. This suggests that the aforementioned evolutionary scheme is not plausible for all HMPOs. Alternatively, HM stars are associated with accretion bursts, which may be responsible for low-mass protostars accreting a significant fraction of their final mass (Audard et al. 2014; Hsieh et al. 2018, 2019). In fact, a variation of mass accretion implied by a luminosity outburst was observed toward the HM protostar S255IR-SMA1 (Caratti o Garatti et al. 2016; Liu et al. 2018).

Zoom In Zoom Out Reset image size

Liu et al. (2020) estimated an envelope mass infalling rate of ∼10⁻⁴ M_⊙yr⁻¹ at r = 200–100 au for S255IR SMA1. This is significantly larger than the upper limit of the disk accretion rate of ∼10⁻⁸ M_⊙yr⁻¹ we measured using the Br-γ emission. Such discrepancies between the envelope and disk accretion rates have been recognized for low-mass protostellar evolution, and support an evolutionary scenario with accretion outbursts triggered by an enhanced disk mass (Calvet et al. 2000).

As shown in Section 4.1, the Br γ emission toward G29.96-0.02, the most massive source in our target list, exhibits a large equivalent width (∼−89 Å) and a narrow linewidth (28kms⁻¹). Martín-Hernández et al. (2003) found extended Br γ emission and suggested that such emission is caused by an ionizing star with a spectral type of O5−O6 V. Of our sources, G29.96-0.02 is the only source with a detection in the Hydrogen Pfund series (Figure 3) and is the most massive among sources associated with an UC H ii region (Table 1). These peculiarities in the observed spectra suggest that this is the unique source for which the Br γ emission is associated with an extended H ii region rather than a compact ionizing region in accretion flows.

As shown in Figure 2, Brγ absorption is detected toward only W3 IRS 5 among our sources (Figure 2), which was not identified by the previous median-resolution observation (R = 4000) in Bik et al. (2012). No clear emission feature was found by either Bik et al. (2012) or this work, implying that the accretion has stopped (or temporarily halted). A Br γ absorption with a narrow linewidth were found to originate from the photospheric absorption in main-sequence stars (Hanson et al. 1996; Wallace & Hinkle 1997). However, it is found that the linewidth of the absorption usually exceeds 100 km s⁻¹ toward massive (O-type) main sequence stars (Hanson et al. 2005), much larger than that in W3 IRS5 (ΔV ∼ 37kms⁻¹). Therefore, the origin of the absorption in W3 IRS5 is still unclear.

There is growing consensus that, for low-mass YSOs, the CO bandhead emission at 2.3–2.4 μm is associated with the the inner disk region where dust grains are sublimated (Dullemond & Monnier 2010). High-resolution spectroscopy of the band emission profiles for some objects show double-peak line profiles and/or a blue-shoulder that are well explained by the disk models (see Najita et al. 2007, for review). Such studies have also been extended toward young stars with intermediate-to-high masses (Ilee et al. 2013, 2014; Murakawa et al. 2013). Medium-to-high resolution spectra of some HMPOs also show band emission with a distinct profile indicating disk rotation (Ilee et al. 2013; Murakawa et al. 2013).

In both low-mass and high-mass cases, the observed CO fluxes or luminosities appear to be correlated with those of Br γ when both emission features are detected (Connelley & Greene 2010; Ilee et al. 2014; Pomohaci et al. 2017). However, the CO bandhead emission usually has a detection rate significantly lower than the Br γ emission: ∼20% (CO) and ∼80% (Br γ) for observations of low-mass Class I protostars by Connelley & Greene (2010), 7% and 70% for observations of Herbig AeBe stars by Ilee et al. (2014), and 17%–34% and 74%–97% for observations of massive protostars by Cooper et al. (2013), and Pomohaci et al. (2017).

In our observations of 11 massive YSOs, the CO bandhead emission v = 2 − 0 was detected only in W33A-MM1 Main and was marginally seen in IRAS 18151-1208, yielding a detection rate of only 9%–18%. A detection rate lower than Cooper et al. (2013), Pomohaci et al. (2017) is probably due to the lower signal-to-noise of our observations, which are partially or fully attributed to the following: (1) a spectral resolution significantly higher than previous observations (R = 35,000–70,000, 500 and 7000 for this work, Cooper et al. 2013; Pomohaci et al. 2017, respectively); and (2) a relatively large continuum excess emission (Section 5.1). For the former, the emission lines are significantly blended at the bandhead, making the feature broad (>0.002 μm; e.g., Ilee et al. 2013; see also the spectrum of W33A-MM1 Main in Figure 3), and a very high angular resolution splits this emission into many spectral pixels, yielding a low signal-to-noise. We smoothed the spectra to increase signal-to-noise, but did not find clear detection toward more sources.

We measure an extinction-corrected CO luminosity of W33A-MM1 Main of 0.04 ± 0.01 L_⊙. The CO/Br-γ luminosity ratio of ∼0.6 roughly agrees with the empirical correlation for Herbig AeBe stars and massive protostars shown in Pomohaci et al. (2017).

Ilee et al. (2018) used LTE analytic disk models to explain the low detection rates of the CO overtone emission. These authors demonstrated that the CO overtone emission would be observed only over a specific range of disk accretion rates. The emission is considered to be associated with the gaseous inner disk where dust is sublimated. A low disk accretion rate ($\dot{M}\lesssim {10}^{-6}\,{M}_{\odot }$ yr⁻¹) would result in a small dust-sublimation region in the inner disk, a correspondingly smaller emitting area for CO, and therefore fainter CO emission with respect to the continuum. On the other hand, a high disk accretion rate ($\dot{M}\gtrsim {10}^{-4}\,{M}_{\odot }$ yr⁻¹) would yield not only a large dust-sublimation region in the inner disk (and therefore the CO emitting region), but also intense thermal dust continuum emission, making the CO emission lower with respect to the continuum.

As derived in Section 5.2, the disk accretion rates inferred from the Br-γ emission range from <10⁻⁸ to ∼10⁻⁴ M_⊙ yr⁻¹ (Table 2, Figure 6). Among these objects, intermediate disk accretion rates (2-3×10⁻⁶ M_⊙ yr⁻¹) were measured for W33A-MM1 Main and are marginally seen in IRAS 18151-1208, i.e., those with detection of the CO bandhead emission. In this context, our detection would be broadly consistent with the model predictions by Ilee et al. (2018); the CO bandhead EWs are at maximam values for ${\dot{M}}_{\mathrm{acc}}\sim {10}^{-5}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. However, a few objects have similar disk accretion rates (e.g., 2–6 ×10⁻⁶ M_⊙ yr⁻¹ for NGC 7538 IRS1, Table 2) but without CO detection. We speculate that their negative detection can be attributed to a disk inclination angle close to edge-on, making the CO emission region undetectable. Even so, we detected CO emission in only two of seven objects with disk accretion rates between ∼10⁻⁶ and ∼10⁻⁴ M_⊙ yr⁻¹, for which Ilee et al. (2018) predicted relatively large band-to-continuum ratios. A large sample with better signal-to-noise ratio is required to further address this issue.

As shown in Section 4.3, the CO absorption observed toward many of our sources is either blueshifted or redshifted, suggesting that the absorption is associated with a relatively cool outflow or inflow. This is in contrast to the discussion in previous studies of a few HM stars, for which the authors attributed the absorption to a circumatellar disk. Davies et al. (2010); Murakawa et al. (2013) observed CO absorption toward a few protostars using integral field spectroscopy with a modest spectral resolution (R ∼ 5500, yielding a velocity resolution of ∼55 km s⁻¹). These authors analyzed the spatial distribution of the velocity centroids of the absorption lines, and revealed kinematics similar to disk rotation. Assuming that the observed kinematics of Δv = 5–10 km s⁻¹ is due to Keplerian rotation, they derived central masses of these objects of 15–30 M_⊙.

In Section 4.3 we measured the absorption equivalent widths of the individual transitions associated with most of our target sources. We apply a simple model to derive the temperature and column density, and discuss the origin of this cool gas component. For a given transition, the absorption line profile I_ν is expressed as

$\begin{eqnarray}&&{I}_{\nu }={I}_{{\rm{c}}}\exp (-s{\phi }_{\nu }{N}_{v,J}),\end{eqnarray} \tag{ 4 }$

where I_c is the continuum flux, s is the integrated cross section, ϕ_ν is the normalized Gaussian function and N_v,J is the column density in the lower energy state (Smith et al. 2009). Under optically thin conditions, the equation can be re-written as

$\begin{eqnarray}&&\displaystyle \frac{{EW}}{\lambda }={{sN}}_{v,J},\end{eqnarray} \tag{ 5 }$

where EW and λ are the equivalent width and wavelength, respectively (Goto et al. 2003).

We derived N_v,J using Equation (5) for individual sources and transitions, and plot them in Figure 7. Using scipy.optmize.curve_fit, we fit these level populations assuming that the CO absorption is associated with an optically-thin LTE gas layer with a single temperature (T). These indicate a temperature for the CO absorption region of 100–200 K for all the sources but W3 IRS 5 and S255IR SMA2, for which the derived temperatures have large uncertainties; these large uncertainties come from the relatively high noise level of the measured EWs (Table 3).

Zoom In Zoom Out Reset image size

Our successful fitting with an LTE model indicates that the gas components observed at these transitions are well thermalized with a density higher than the critical density ${n}_{{{\rm{H}}}_{2}}\sim 5\times {10}^{5}\,{\mathrm{cm}}^{-3}$ of CO at v = 0, J = 8–15 (Yang et al. 2010). Given CO column densities of ∼(0.5–2.9) × 10¹⁸cm⁻², the thicknesses of the warm gas components are less than 700–3900 au taking the canonical CO/H₂ abundance of ∼10⁻⁴. If the warm gas component has a higher density, e.g., ∼10⁷cm⁻³, the absorption would trace an inner region within 30–200 au.

Given the temperature of these gas components, we estimate their distance to the central heating source. Given the source luminosity (Table 1), the dust temperature can be expressed as a function of radius if we assume that the power radiated is equal to that absorbed in a spherical symmetric cloud (i.e., ${T}_{{\rm{d}}}\propto \tfrac{{L}^{1/5}}{{R}^{2/5}}$, see Equation (31) in Goldreich & Kwan (1974)). As a result, we estimate that the absorption occurred inside a region with a typical radius of a few hundred aus ( ≲ 200–600 au. depending on source), except for W3 IRS 5; with the large uncertainty in temperature, the radius of the CO component is not constrained well at ∼10–3000 au.

In practice, absorption components at different velocities would be associated with gas at different radii. Figure 4 shows that S255 SMA1 and SMA2 are associated with redshifted CO absorption up to ∼30 km s⁻¹. One would estimate that the absorption component at the highest velocity would be associated with r ∼ 6 au, assuming freefall motion toward the HPMO with a stellar mass of 20M_⊙. Observations at a higher signal-to-noise would allow us to better investigate temperature distributions of CO gas at different radii .

The equivalent widths of Br γ line emission have been measured by the surveys of Cooper et al. (2013) and/or Pomohaci et al. (2017). The measured equivalent widths are −4.8, −1.4, −0.3, −5.1 and −4.0 Å for AFGL 490, IRAS 18151-1208, AFGL 2591 VLA3, NGC 7538 IRS1, and IRAS 23151+5912, respectively. These measurements are consistent with ours within a factor of ∼2.5 except for AFGL 2591 VLA3. While Pomohaci et al. (2017) found an equivalent width of −0.3 Å toward AFGL 2591 VLA3, we give an upper limit of 0.03 Å (Figure 2).

Bik et al. (2012) conducted a K_s -Band spectroscopic survey of young stars and YSOs in the W3-Main star-forming regions in 2009. These authors did not detect Br-γ emission associated with W3 IRS 5. In our high-resolution spectrum, we found a Br γ absorption feature at the systemic velocity with an equivalent width of 0.58 ± 0.04 Å and a narrow linewidth of 36.6kms⁻¹.

Ishii et al. (2001) measured the equivalent width of Brγ to be −12.6 Å which is larger than our derived value by a factor of ∼2.3 (Table 2).

As mentioned in C.4 S255IR SMA1 experienced a recent accretion burst between 2007 and 2016 (Caratti o Garatti et al. 2016; Liu et al. 2018). Due to high extinction and/or strong veiling, no clear Br γ emission was detected toward the source center (S255IR IRS3 in Caratti o Garatti et al. 2016), consistent with our result in S255IR SMA1.

The equivalent width of Brγ is measured to be −5 Å and −3.17 Å in Davies et al. (2010) and this work, respectively. This small difference might come from the different spectral resolutions or time variability.

G29.96-0.02 shows distinct features from our survey sample. It has an unusual strong Br γ line emission with an equilvalent width of ∼−90 Å (Table 2). Such strong Br γ emission was also observed in Martín-Hernández et al. (2003) with the median-resolution observation using VLT. In addition, several helium recombination lines were detected in the literature. Pfund series emission is also detected in G29.96-0.02 (Figure 3). These results suggest that G29.96-0.02 contains a large amount of highly excited material. The peculiarity of this source is probably related to the fact that it is the brightest (L = 8 × 10⁵ L_⊙) and the most massive (M_star ∼ 33 − 56M_⊙) of our sources (Table 1).

Ks-band spectroscopic observations covering CO v = 0–2 were taken by Cooper et al. (2013) and/or Pomohaci et al. (2017) in these sources, and thus we discuss them concurrently. CO bandhead emission is not detected in these sources in the literature. Cooper et al. (2013) and Pomohaci et al. (2017) executed spectroscopy for these sources, but they did not detect the CO bandhead emission. These non-detections are mostly consistent with our work (except for IRAS 18151-1208, see C.6) as our observations provide a higher spectral resolution but a lower sensitivity.

Bik et al. (2012) conducted Ks-Band spectroscopic observations toward the sources in the W3-Main star-forming region. Among their targets, 15 candidates were identified as OB stars through spectral classification. The photospheres of OB stars are detected from the more evolved diffuse H ii regions, indicating that most of these OB stars have finished their star formation. Among their targets, W3 IRS 5 is the only source that is found to be still accreting. CO bandhead emission was not detected in Bik et al. (2012) or in our work (with an rms noise level of 3 Å integrating over 2.289–2.313 μm).

Ishii et al. (2001) did not detect CO bandhead emission. While Ishii et al. (2001) measured the upper limit of the CO v = 2 − 0 bandhead emission to be > −6.2 Å, our data gives a constraint of > −6.4 Å given the integrated range of 2.289–2.313 μm.

S255IR SMA1 has been found to have experienced an accretion burst on the basis of near-infrared spectroscopic observations in 2007 and 2016 (Caratti o Garatti et al. 2016). Soon after that, Liu et al. (2018) reported that the submillimeter continuum emission dimmed from 2016 to 2017. Caratti o Garatti et al. (2016) show the spectra of the source and the outflow cavity in the pre-burst phase (2007) and the outburst phase (2016). As a result, the CO bandhead emission is only seen toward the outflow cavity in the pre-burst phase. This is consistent with our result of non-detection toward the source as also reported in Cooper et al. (2013) from observation in 2005.

CO bandhead emission is detected in Davies et al. (2010) and Ilee et al. (2013) as well as in our observations. With high-resolution observation (R ∼ 30, 000), Ilee et al. (2013) find that CO bandhead emission is well fitted by a model of a Keplerian rotation disk; the CO emission traces the disk from the inner region from a few au to a few tens of au. In our spectra we measured an equivalent width of the CO v = 2–0 bandhead emission of −2.33 ± 0.69 Å.

The v = 2 − 0 and 3 − 1 CO bandhead emission was observed by Murakawa et al. (2013). At a lower spectral resolution (R∼5500), Davis et al. (2004) detected strong CO bandhead emission between v = 2 − 0 and v = 7−5. This observations were taken with a much lower resolution and the emission is interpreted as inner disk surface or outer accreting funnel-flow. In Cooper et al. (2013), the CO emission is not detected, likely due to the sensitivity.

As shown in Section 4.2, we have marginally detected the v = 2 − 0 bandhead emission. It is difficult to discuss detection of the CO v = 3 − 1 bandhead emission since the continuum is not well defined (see Table 3).