Massive pre-main sequence stars in M17

TL;DR

This study uses optical to near-infrared spectroscopy to analyze massive young stellar objects in M17, confirming their pre-main sequence status, characterizing their disks, and linking their properties to evolutionary models of massive star formation.

Contribution

First spectroscopic confirmation and detailed characterization of pre-main sequence massive stars in M17, linking observations to evolutionary models and disk structures.

Findings

Six PMS stars confirmed with contracting radii

Disks show structured velocity profiles from emission lines

Massive stars are on track to reach the main sequence with high accretion rates

Abstract

The formation process of massive stars is still poorly understood. Massive young stellar objects (mYSOs) are deeply embedded in their parental clouds, they are rare and thus typically distant, and their reddened spectra usually preclude the determination of their photospheric parameters. M17 is one of the best studied HII regions in the sky, is relatively nearby, and hosts a young stellar population. With X-shooter on the ESO Very Large Telescope we have obtained optical to near-infrared spectra of candidate mYSOs, identified by Hanson et al. (1997), and a few OB stars in this region. The large wavelength coverage enables a detailed spectroscopic analysis of their photospheres and circumstellar disks. We confirm the pre-main sequence (PMS) nature of six of the stars and characterise the O stars. The PMS stars have radii consistent with being contracting towards the main sequence and are surrounded by a remnant accretion disk. The observed infrared excess and the (double-peaked) emission lines provide the opportunity to measure structured velocity profiles in the disks. We compare the observed properties of this unique sample of young massive stars with evolutionary tracks of massive protostars by Hosokawa & Omukai (2009), and propose that these mYSOs near the western edge of the HII region are on their way to become main-sequence stars ($\sim 6 - 20$ $M_{\odot}$) after having undergone high mass-accretion rates (${\dot{M}_{\rm acc}} \sim 10^{-4} - 10^{-3}$ $M_{\odot}$ $\rm yr^{-1}$). Their spin distribution upon arrival at the zero age main sequence (ZAMS) is consistent with that observed for young B stars, assuming conservation of angular momentum and homologous contraction.

Figures

Tables

Table 1. Table 2: Summary of the properties of the sources observed with X-shooter.

Source	Emission	CO	NIR	Classification
	lines	bandhead	excess
B111	-	-	-	O star
B163	-	+	+	mYSO
B164	-	-	-	O star
B215	-	-	*555NIR excess only longwards of 3 $\mu$m, not in the X-shooter wavelength range.	B star?
B243	+	+	+	mYSO
B253	-	-	-	B star
B268	+	+	+	mYSO
B275	+	+	+	mYSO
B289	-	-	*5	O star?
B311	-	-	-	O star
B331	+	+	*5	mYSO
B337	+	-	+	mYSO

Table 2. Table 3: Stellar properties derived from best-fit FASTWIND parameters and the SED fitting. The values for the stellar radius ( R ⋆ subscript R ⋆ \rm R_{\star} ) cited here are the ones obtained via SED fitting. The sources classified as PMS stars in this work are in bold face, and the two sources for which we cannot confirm the PMS nature are underlined.

Name	Sp. Type	Sp. Type	$T_{\mathrm{eff}}$	$\log g$	$v\sin i$	$R_V$	$A_V$	$\log L/L_{\odot }$	$R_{\star }$	$M_{\mathrm{spec}}$	RV
	Hanson	This work	K	cm $s^{-2}$	km $s^{-1}$	(SED)	mag		$R_{\odot }$	$M_{\odot }$	km $s^{-1}$
B111	O5 V	O4.5 V	${42850}_{-1600}^{+3050}$	${3.82}_{-0.27}^{+0.13}$	${170}_{-28}^{+42}$	${3.7}_{-0.3}^{+0.4}$	${5.4}_{-0.3}^{+0.4}$	${5.71}_{-0.04}^{+0.08}$	${13.0}_{-0.6}^{+0.4}$	${42}_{-13}^{+12}$	3.4$\pm$0.8
B163	?	kA5	$8200$	$-$	$-$	${4.0}_{\downarrow }^{\uparrow }$	${13.2}_{\downarrow }^{\uparrow }$	${2.95}_{\downarrow }^{\uparrow }$888Values from SED fitting	${10.1}_{\downarrow }^{\uparrow }$	$-$	$-$
B164	O7-O8	O6 Vz	${39100}_{-2350}^{+1150}$	${4.01}_{-0.27}^{+0.16}$	${108}_{-22}^{+30}$	${3.8}_{-0.3}^{+0.3}$	${8.3}_{-0.5}^{+0.4}$	${5.09}_{-0.07}^{+0.03}$	${9.3}_{-0.2}^{+0.4}$	${22}_{-8}^{+10}$	-1.8$\pm$2.3
B215	B0.5 V	B0-B1 V	${23500}_{-3450}^{+6100}$	${3.67}_{\downarrow }^{+0.5}$	${208}_{-38}^{+40}$	${4.1}_{-0.1}^{+0.2}$	${7.6}_{-0.1}^{+0.3}$	${3.98}_{-0.15}^{+0.23}$	${6.7}_{-1.04}^{+0.94}$	$6_{-2}^{+9}$	12.3$\pm$1.0
B243	B3 V	B8 V	${13500}_{-1250}^{+1350}$	${4.34}_{-0.3}^{\uparrow }$	${110}_{\downarrow }^{+106}$	${4.7}_{-0.8}^{\uparrow }$	${8.5}_{-1.0}^{\uparrow }$	${3.21}_{-0.06}^{+0.07}$	${7.5}_{-0.8}^{+1.0}$	${44}_{-22}^{+21}$	20.2$\pm$1.5
B253	B3V	B3-B5 III	${13000}_{-1200}^{+2100}$	${3.09}_{\downarrow }^{+0.52}$	${318}_{\downarrow }^{\uparrow }$	${4.1}_{-0.2}^{+0.1}$	${6.5}_{-0.2}^{+0.1}$	${3.21}_{-0.06}^{+0.1}$	${7.2}_{-1.3}^{+1.1}$	$3_{-1}^{+5}$	9.5$\pm$0.8
B268	B2 V	B9-A0	${12250}_{-1000}^{+850}$	${3.99}_{-0.38}^{\uparrow }$	${36}_{\downarrow }^{+126}$	${4.6}_{-0.8}^{\uparrow }$	${8.1}_{-1.0}^{\uparrow }$	${3.24}_{-0.05}^{+0.04}$	${8.8}_{-0.8}^{+1.2}$	${31}_{-15}^{+1}$	4.3$\pm$1.3
B275	late-O999From SED in H97	B7 III	${12950}_{-650}^{+550}$	${3.39}_{-0.11}^{+0.06}$	${88}_{\downarrow }^{+68}$	${3.8}_{-0.8}^{+0.7}$	${6.7}_{-1.0}^{+0.8}$	${3.37}_{-0.03}^{+0.02}$	${11.7}_{-0.5}^{+0.67}$	$8_{-2}^{+1}$	-11.2$\pm$1.4
B289	O9.5 V	O9.7V	${33800}_{-2150}^{+2500}$	${3.93}_{\downarrow }^{+0.35}$	${154}_{-34}^{+38}$	${3.7}_{-0.5}^{+0.4}$	${8.3}_{-0.8}^{+0.6}$	${4.84}_{-0.07}^{+0.08}$	${8.3}_{-0.4}^{+0.4}$	${19}_{-5}^{+18}$	-2.3$\pm$2.0
B311	O9-B2 V	O8.5 Vz	${35950}_{-1750}^{+2400}$	${3.95}_{\downarrow }^{+0.36}$	${42}_{-26}^{+12}$	${3.3}_{-0.3}^{+0.3}$	${6.1}_{-0.4}^{+0.4}$	${4.79}_{-0.06}^{+0.07}$	${7.6}_{-0.3}^{+0.2}$	${14}_{-4}^{+16}$	4.2$\pm$0.4
B331	B2 V101010Hoffmeister et al. (2008)	late-B	$13000$	$-$	$-$	${4.6}_{-0.5}^{+0.5}$	${13.3}_{-0.9}^{+0.9}$	${4.10}_{\downarrow }^{+0.37}$8	${21.8}_{-7.2}^{+9.6}$	$-$	13.9$\pm$1.5
B337	mid-B9	late-B	$13000$	$-$	$-$	${3.7}_{-0.7}^{+0.9}$	${13.6}_{-0.9}^{+0.9}$	${3.38}_{\downarrow }^{+0.02}$8	${8.7}_{-3.3}^{+5.0}$	$-$	7.5$\pm$0.6

Table 3. Table 4: Projected disk rotational velocity measured from half the peak to peak separation of the lines originating in the circumstellar disk.

	$v_{\mathrm{disk}}$ [km/s]
Line	B243	B268	B275	B331	B337
OI-844.6	76	-	38	-	-
CaII-849.8	-	112	59	-	48
Pa-16	86	-	-	83	-
CaII-854.2	-	119	57	-	41
Pa-15	91	-	-	85	-
Pa-14	86	-	40	81	-
CaII-866.2	-	122	54	-	35
Pa-13	83	-	-	80	-
Pa-12	84	-	43	-	-
Pa-11	79	-	43	64	-
Pa-9	78	-	45	63	-
Pa$\gamma$	63	63	29	-	-
Pa$\beta$	62	58	28	-	-
Br-16	90	-	-	74	-
Br-14	92	58	-	70	-
Br-12	73	-	44	64	-
Br-11	-	-	40	53	-
Br-10	70	72	-	48	-
Br$\gamma$	77	62	-	-	-

Table 4. Table 5: Age of our targets estimated by comparing their position in the HRD with PMS tracks from Hosokawa & Omukai ( 2009 ) (third column) and main sequence tracks from Brott et al. ( 2011 ) (fourth column).

star	PMS	PMS	MS
	track	lifetime	lifetime
	[$M_{\odot }$]	[Myr]	[Myr]
B111	25	$-$	$0-2$
B163	6	0.14	$-$
B164	25	0.01	$0-2$
B215	10	0.03	$-$
B243	6	0.20	$-$
B253	6	0.20	$-$
B268	6	0.20	$-$
B275	8	0.04	$-$
B289	20	0.02	$-$
B311	20	0.02	$0-2$
B331	12	0.02	$-$
B337	8	0.04	$-$

Full text

1 Introduction

In the past decades significant progress has been made regarding the understanding of star formation. However, the formation of massive stars – the only mode of star formation observable in external galaxies – remains a key unsolved problem (Zinnecker & Yorke, 2007; Beuther et al., 2007; Tan et al., 2014). One of the main reasons is that observations of the earliest phases in the life of a massive star are scarce. This is explained by the expected short formation timescale ($\sim 10^{4}-10^{5}$ yr) of massive stars, and the severe extinction ($A_{V}\sim 10-100$ mag) by which the surrounding gas and dust obscure their birth places from view. Additionally, massive stars are rare and, as a consequence, located at relatively large distances.

Observational and theoretical evidence is accumulating that the formation process of massive stars is through disk accretion, similar to low-mass stars. This despite the physical processes acting near hot, massive stars (e.g., radiation pressure, ionisation) that counteract the accretion flow onto the forming star (e.g., Wolfire & Cassinelli, 1987; Keto & Wood, 2006; Krumholz et al., 2009; Klassen et al., 2016). In the accretion scenario, one expects that soon after the gravitational collapse of the parental cloud the pre-stellar core is surrounded by an extended accretion disk (size $\sim 1000$ AU; e.g., Beltrán et al., 2004; Cesaroni et al., 2005; Sandell & Wright, 2010; Beltrán & de Wit, 2016). Apparently this is not a very efficient process, as most of these high-mass protostars seem to drive outflows (e.g., Zhang et al., 2001). Already before arriving on the main sequence, the young massive star is expected to produce a strong ultraviolet (UV) radiation field ionising its surroundings (Bernasconi & Maeder, 1996). At this stage the object will become detectable at radio and infrared wavelengths through the heated dust and recombination of hydrogen in an expanding hyper- or ultra-compact H ii (UCHii) region (e.g., Churchwell, 2002). Unfortunately, only very little is known about the physical properties of the central (massive) stars at this stage of formation (e.g. Martín-Hernández et al., 2003).

Given the short formation time of massive stars, the accretion rate must be high (up to $10^{-3}$ $M_{\odot}$ yr*-1*, Hosokawa et al., 2010). At these accretion rates, a massive young stellar object (mYSO) is expected to bloat up to about 100 $R_{\odot}$, resulting in a relatively low effective temperature and modest UV luminosity. The accretion process is unlikely to be constant with time; the blob structure observed in Herbig Haro outflows indicates that strong variations occur in the mass in- and outflow rate of young intermediate-mass stars (e.g. Reipurth & Bally, 2001; Ellerbroek et al., 2013b). Simulations of different accretion models applied to mYSOs show that the accretion rate increases as the mYSO grows in mass (Davies et al., 2011). Once the accretion rate diminishes, the “bloated” pre-main-sequence (PMS) star contracts to the main sequence on the Kelvin-Helmholtz timescale. Recently, candidates for such bloated, massive PMS stars have been spectroscopically confirmed (Ochsendorf et al., 2011; Bik et al., 2012; Ellerbroek et al., 2013a; Chen et al., 2015).

The environment (multiplicity, clustering) and the large distances to massive (proto)stars make their circumstellar disks very difficult to resolve (Beltrán & de Wit, 2016). In contrast, disks around intermediate-mass stars (2 ¡ $M_{\star}$ ¡ 8 $M_{\odot}$) have been characterised using (sub)millimeter as well as optical and near-infrared spectroscopy (see Bik & Thi, 2004; Wheelwright et al., 2010; Ellerbroek et al., 2011; Alcalá et al., 2014). The rotational structure of the disk provides information about the role of magnetic fields that can slow down the rotation rate below pure Keplerian. Super-Keplerian rotation of the disk could indicate that the inner disk contributes significantly to the gravitational potential of the system (Beltrán & de Wit, 2016). In the mass range from 20–30 $M_{\odot}$ circumstellar disks have been detected, but there is no consensus about their general properties.

M17, located in the Carina-Sagittarius spiral arm of the Galaxy is one of the best studied giant H ii regions ($L=3.6\times 10^{6}L_{\odot}$, Povich et al., 2007). Its distance has been accurately determined by measuring the parallax of the $\rm CH_{3}OH$ maser source G15.03-0.68: $d=1.98^{+0.14}_{-0.12}$ kpc (Xu et al., 2011). The bright blister H ii region is embedded in a giant molecular cloud complex, and divides the molecular cloud into two components: M17 South and M17 North, containing a total gas mass (molecular and atomic) of about $6\times 10^{4}$ $M_{\odot}$ (Povich et al., 2009). The centre of the H ii region hosts the cluster NGC 6618 including 16 O stars and over 100 B stars (Chini et al., 1980; Hoffmeister et al., 2008), providing the ionising power for the H ii region. Many of the OB stars are suspected to be binaries (Hoffmeister et al., 2008; Povich et al., 2009), which would explain why they are more luminous than expected from their spectral type. NGC 6618 is a young cluster, its age is estimated at about 1 Myr (Hanson et al., 1997; Broos et al., 2007; Povich et al., 2009), while the surrounding molecular cloud hosts pre-main-sequence stars. The photo-dissociation region to the southwest of NGC 6618 includes several candidate mYSOs: the hyper-compact H ii region M17-UC1 surrounded by a circumstellar disk (Nielbock et al., 2007), and IRS5, a young possibly quadruple system of which the primary star, IRS5A (B3–B7 V/III), is a high-mass protostellar object ($\sim$9 $M_{\odot}$, $\sim 10^{5}$ yr, Chen et al., 2015). Hanson et al. (1997) (hereafter H97, ) discovered a sample of high-mass (5–20 $M_{\odot}$) young stellar objects in the same area. Their SEDs display a near-IR excess, some show CO bandhead and (double-peaked) Pa $\delta$ emission, likely due to the presence of a circumstellar disk. These objects are the subject of study in this paper.

We are now on the verge of being able to quantitatively test the predictions of massive star formation models. Several samples of mYSOs have been presented (e.g., Bik et al., 2006; Urquhart et al., 2011; Ochsendorf et al., 2016). However, so far only for a handful of mYSOs the photospheric spectrum has been detected. We have obtained high-quality spectra with X-shooter on the ESO Very Large Telescope (VLT) of candidate mYSOs, first identified by H97, and a few other O and B stars in M17. The large wavelength coverage of the X-shooter spectra (and additional infrared photometry) allows for a detailed analysis of the spectral energy distribution, including an infrared excess produced by the dust component in the accretion disk. The observed (double-peaked) emission lines provide the opportunity to study the dynamical structure of the gaseous component of the disk.

This paper is organised as follows. In Sect. 2 we present our target sample, observations and data reduction procedure. In Sect. 3 we provide the spectral classification of the (pre-) main-sequence OB stars in our sample. After that (Sect. 4) we present their spectral energy distributions. To further constrain the properties of the mYSOs we modeled the optical spectra using FASTWIND. These results are shown in Sect. 5, where we also place the PMS objects in the Hertzsprung-Russell diagram and compare their position with PMS tracks and isochrones. Additional evidence for the presence of circumstellar disks is obtained from near-infrared spectroscopy (Sect. 6). We discuss the age distribution of the massive PMS population in M17, the extinction towards this region, the presence of circumstellar disks, and the spin properties of the young massive stars in Sect. 7. We summarise our conclusions in Sect. 8.

2 VLT/X-shooter spectroscopy

We obtained optical to near-infrared (300–2500 nm) spectra of candidate mYSOs, and some OB stars in M17. We used the X-shooter spectrograph mounted on UT2 of the VLT (Vernet et al., 2011). The spectra of B275 and B358 were obtained as part of the X-shooter science verification program; B275 has been subject of study in Ochsendorf et al. (2011). For some objects multiple spectra have been taken. A log of the observations is listed in Table 1.

2.1 Target sample

We selected our targets based on the study of H97 who obtained $K$-band spectra of the young massive star population of M17, complemented by optical and 1 $\mu$m spectroscopy for some of the moderately reddened objects. Several of the targets have also been studied by Hoffmeister et al. (2008) and some are included in the studies of Nielbock et al. (2001), Broos et al. (2007), and Povich et al. (2009). H97 reported that B111, B164, B253, and B311 do not present any NIR excess nor emission lines in their spectra and they are likely OB-star members of the central cluster NGC 6618. It is likely that all our targets belong to NGC 6681, at least the (massive) main sequence stars.

Based on the NIR excess derived from their position in the NIR colour-colour diagram (CCD, see also Sect. 2.2) and displayed by the spectral energy distribution (SED, see also Sect. 4), and the apparent presence of CO bandhead and/or double-peaked emission lines, the candidate mYSOs according to H97 are B163, B243, B268, B275, B289, B331, and B337. We also observed B215 (IRS15): this mYSO candidate has been studied by Nielbock et al. (2001). The spectral types from H97 and Hoffmeister et al. (2008) are listed in Table 3. We obtained an X-shooter spectrum of B358, but according to Chen et al. (2013) and H97 the spectrum is that of an early-/mid-G supergiant. It has not been detected in X-rays (Broos et al., 2007) and it is likely a background post-AGB star. From examining the X-shooter spectrum we support this conclusion.

Figure 1 (left panel) shows a colour composite image of M17 based on 2MASS $J$ (blue), $H$ (green) and $K$ (red) observations (Skrutskie et al., 2006). The stars for which we collected X-shooter spectra are labeled. The area includes the bright nebular emission produced by the H ii region and the central cluster NGC 6618. The mYSO candidates are predominantly located near the excited rim of the surrounding molecular cloud. This observation suggests that those objects could be the product of a later phase in the star formation process. In the right panel a colour composite is shown based on images taken with Spitzer (GLIMPSE survey; IRAC1 - 3.6 $\mu$m, green: IRAC4 - 8.0 $\mu$m, Benjamin et al., 2003) and Herschel (red: PACS - 70 $\mu$m, Pilbratt et al., 2010; Poglitsch et al., 2010). This image provides insight into the (dusty) surroundings of the H ii region. The 3.6 $\mu$m image traces the illuminated edges of the cloud irradiated by the massive stars, whereas the 8 $\mu$m emission is usually associated with emission features produced by polycyclic aromatic hydrocarbons (PAHs) and/or warm (100-200 K) dust that emits at these wavelengths when heated by ultraviolet radiation (Povich et al., 2009).

Throughout this paper we have color-coded the labels to present our results: the objects whose PMS nature is confirmed are plotted in blue, the OB stars in grey (black in figure 1), the objects for which we cannot confirm the PMS nature in green and the post-AGB star in red.

2.2 Colour-colour and colour-magnitude diagram of X-shooter targets in M17

The ($K$, $J-K$) colour-magnitude diagram is shown in Fig 2; the $J$, $H$, and $K$-band magnitudes are from the 2MASS database (Cutri et al., 2003). Similar diagrams have been presented by H97 and Hoffmeister et al. (2008), but not using the recent distance estimate of 1.98 kpc (Xu et al., 2011). The zero-age main sequence (ZAMS) is marked with a drawn line; the dashed line represents the reddening vector. The ZAMS properties for spectral types later than B9 were taken from Cox (2000) and for the earlier type stars from H97. We over-plot the properties for the O dwarf stars from Martins & Plez (2006) for reference. We used the Cardelli et al. (1989) extinction law in order to plot the reddening lines. The candidate mYSOs have a large $(J-K)$ colour and a bright $K$-band magnitude.

Figure 3 presents the near-infrared colour-colour diagram including the OB-type stars for which an X-shooter spectrum has been obtained. In the colour-colour diagram a near-infrared excess becomes very apparent, and thus the possible presence of circumstellar material. The solid line represents the unreddened main sequence. The dashed lines indicate the reddening lines for an O3 dwarf and an M0 giant. Most of the mYSO candidates are located to the right of the reddening line of an O3 V star, consistent with the expectation that these objects host a circumstellar disk. B163, B289, and B331 are candidate mYSOs according to H97; however, they do not exhibit a near-infrared excess in the colour-colour diagram. Nevertheless, B331 and B163 are mYSO candidates according to our criteria: they present double-peaked emission lines, CO bandhead emission and/or infrared excess beyond the $K$-band. We will discuss the spectral energy distribution and the possible presence of an infrared excess in Sect. 4, also at mid-infrared wavelengths.

2.3 Reduction VLT/X-shooter spectra

The X-shooter spectra were obtained under good weather conditions, with seeing ranging from 0$\aas@@fstack{\prime\prime}$5 and 1″ and clear sky. With the exception of the 2012 B289 spectrum and the 2009 B275 science verification spectrum, the spectrograph slit widths used were 1″ (UVB, 300 – 590 nm), 0$\aas@@fstack{\prime\prime}$9 (VIS, 550 – 1020 nm), and 0$\aas@@fstack{\prime\prime}$4 (NIR, 1000 – 2480 nm), resulting in a spectral resolving power of 5100, 8800, and 11300, respectively. The slit widths for the 2010 B275 observations were 1$\aas@@fstack{\prime\prime}$6, 0$\aas@@fstack{\prime\prime}$9, and 0$\aas@@fstack{\prime\prime}$9 resulting in a resolving power of 3300, 8800, and 5600, respectively. For the 2012 B289 observations we used the 0$\aas@@fstack{\prime\prime}$8, 0$\aas@@fstack{\prime\prime}$7, and 0$\aas@@fstack{\prime\prime}$4 slits corresponding to a resolving power of 6200, 11000, and 11300 for the UVB, VIS and NIR arms, respectively. The spectra were taken in nodding mode and reduced using the X-shooter pipeline (Modigliani et al., 2010) version 2.7.1 running under the ESO Reflex environment (Freudling et al., 2013) version 2.8.4.

The flux calibration was obtained using spectrophotometric standards from the ESO database. We then scaled the NIR flux to match the absolutely calibrated VIS spectrum. The telluric correction was performed using the software tool molecfit v1.2.0444http://www.eso.org/sci/software/pipelines/skytools/molecfit (Smette et al., 2015; Kausch et al., 2015). Parts of the spectra are shown in Figs. 4 to 6 and the full spectrum is available in the on-line material of this paper. We used the xsh_scired_slit_nod recipe to reduce the data, meaning that the sky subtraction is performed by subtracting the two different nodding positions. We note that due to the spatial variation of the nebular lines along the slit, residuals from the sky subtraction are present in some of the the reduced spectra (e.g., [N ii] 6548,6584 Å, [S ii] 6716,6731 Å and [S iii] 9069,9532 Å, and the H lines).

3 Spectra and spectral classification

Traditionally, the blue spectral region (400 – 500 nm) is used to perform the spectral classification of OB-type stars. In the case of the M17 sources studied in this paper the blue region is severely affected by extinction, but is for most of our objects detected in the X-shooter spectrum (Fig. 4). Many targets show the full hydrogen Balmer, Paschen and Brackett (up to Br$\gamma$) series, He i and He ii lines (the latter in the O-type stars), the Ca triplet lines, and diffuse interstellar bands. Figs. 5 and 6 provide an overview of the Ca triplet region (820–880 nm) and the $K$-band (2100–2400 nm, including Br$\gamma$ and the CO bandheads) for all targets. With the exception of B163, B331 and B337 which are strongly affected by interstellar absorption in the blue ($A_{V}>10$), we classified the stars based on the strength of the H, He and metal lines in the UVB spectrum. For the O stars we visually compared our spectra with the previously catalogued and published spectra by the Galactic O Star Spectroscopic Survey (GOSSS; Sota et al., 2011, 2014; Maíz Apellániz et al., 2016). For the B stars we used the criteria and spectra published by Gray & Corbally (2009).

We looked for the presence of [O i] 6300. This forbidden line is associated with bi-polar outflows/jets, disk winds, and the disk surface in Herbig Ae/Be stars. The [O i] line is thought to originate from the region where the UV radiation from the star impinges on the disk surface (Finkenzeller, 1985; van der Plas et al., 2008; Mendigutía et al., 2011). We also paid special attention to the luminosity subclass Vz, which refers to objects with a substantially stronger He ii  4686 line. These objects are hypothesised to be the youngest optically observable O type stars (the ’z’ stands for ZAMS; Sabín-Sanjulián et al., 2014; Arias et al., 2016). Given the young age of M17, one might expect such stars to be present.

Following this classification scheme we are able to determine the spectral types with an accuracy of one subtype, unless stated otherwise in Sect. 3.1. We constrained the luminosity class of the stars using the surface gravity derived in Sect. 5. An overview of the spectral features described in this section is presented in Table 2 and the results of the spectral classification are listed in Table 3. In most cases the spectral classification agrees well with the one provided in the literature by H97 and Hoffmeister et al. (2008).

3.1 Individual sources

In this section we discuss the X-shooter spectra of the individual targets. In parentheses we list the updated spectral types corresponding to this work. The spectra are displayed in full in Appendix A.

B111 (O4.5 V): This clearly is an early O-type spectrum: the He ii lines at 4200, 4541, and 4686 Å are strong, consistent with an O4.5 V type star. The spectrum does not show He i 4387, and He i 4471 is weak. The hydrogen Balmer (including H$\alpha$) and Paschen series lines are strongly in absorption. Using its NIR spectrum, H97 classified B111 as a kO3-O4 based on the He ii, N ii lines and the lack of C iv. According to our $K$-band spectrum B111 is a kO5-O6 type star, based on the presence of He ii lines and weak, if any, C iv emission. The spectrum does not include any (nebular) emission lines.

B163 (kA5): Due to severe interstellar extinction it is not possible to identify any photospheric features in the UVB arm of the X-shooter spectrum. Longward of $\sim$8500 Å  the first Paschen line becomes detectable. The NIR spectrum also shows the Brackett series in absorption, no elements other than H are seen. Based on this spectrum we assign a mid-A spectral type. Two of the CO first overtone transitions are detected, indicating the presence of a circumstellar disk. The spectrum includes nebular emission lines in the hydrogen series, He i lines (e.g., 10830 Å), and [N iii] 9069,9532 Å. This points to an energetic ionising radiation field as expected given the location of B163 near the core of NGC 6618 (Fig. 1). According to our criteria B163 is a mYSO.

B164 (O6 Vz): The presence of He ii 4541, 4686, and the fact that these lines are stronger than He i 4387 and 4471 demonstrate that B164 is an O-type star. The ratio He ii 4200 / He i (+ii) 4026 Å is near unity which indicates a spectral type O6 V. Based on the ratio EW(He ii 4686) / EW(He ii 4542) $>1.10$ we conclude that it belongs to the Vz class. From the $K$-band classification we conclude that B164 is a kO5-O6 type star. H97 classified B164 as kO7-O8 based on the N iii and He i lines. No emission lines are detected.

B215 (IRS15; B0-B1 V): No He ii lines are detected, which indicates that B215 must be of later spectral type than B0. The hydrogen lines of the Balmer, Paschen and Brackett series are in absorption, and the He i absorption lines are weak (e.g., 4922). The absence of Mg ii 4481 suggests that the spectral type is earlier than B1; therefore, we conclude that this star is of spectral type B0-B1. The $K$-band classification is consistent with kO9-B1. Nielbock et al. (2001) classified this object as an extreme class I source, and Chini et al. (2006) pointed out that it is a $\sim 26$ $M_{\odot}$ star surrounded by an extensive remnant disk. The spectrum contains nebular emission lines in the hydrogen series, He i lines and forbidden emission lines [N ii] 6548,6584, [S ii] 6716,6731 and [S iii] 9069,9532. We also identify [O i] 6300 emission: a broad component ($\sim$100 km/s) on top of a nebular component. By inspecting the 2D frames we conclude that the emission seen in H$\alpha$ has a nebular origin. From the X-shooter spectrum alone we cannot conclude that this is a mYSO.

B243 (B8 V): The hydrogen lines of the Balmer series are in absorption. The strongest Balmer lines show a central emission component; H$\alpha$ is strongly in emission. The Paschen series absorption lines are filled in with double-peaked emission lines. The Brackett series lines mainly exhibit a double-peaked emission component. Furthermore, the ‘auroral’ O i 7774 and 8446 lines show prominent double-peaked emission. The Ca ii triplet is not present. The spectrum includes He i absorption lines (e.g. 4471, 5876, 10830); the ratio He i 4026 / He i 4009 is around two. The He i 4471 and Mg ii 4481 line ratio is close to unity which points towards spectral type B8. The $JHK$ CCD indicates NIR excess, weak CO-bandhead emission is detected, and it presents [O i] 6300 emission. Therefore we conclude that B243 is a mYSO.

B253 (B3-B5): The spectrum displays strong and broad hydrogen absorption lines, with a central (nebular) emission component. He i absorption lines are present (e.g., 4026, 4471, 4922, 10830 Å), also with a central (nebular) emission component; the He i lines are most prominent in the blue part of the spectrum. No He ii lines are detected, nor the C ii 4267 Å line; therefore, we conclude that B253 is a B3-B5 type star. The NIR CMD and CCD show no evidence for a NIR excess. Its $K$-band spectrum is kB5 or later. Many forbidden emission lines are present (e.g., [O iii], [N ii], [S ii]) that, together with the He i emission, indicate a high degree of ionisation of the surrounding H ii region.

B268 (B9-A0): The hydrogen Balmer and Paschen series are prominently in absorption and include a central nebular emission component. The H$\alpha$ and H$\beta$ profiles are filled in by a circumstellar emission component. The He i lines are very weakly present, the nebular emission component dominates the line profiles. The fact that He i 4471 and Mg ii 4481 have almost the same strength indicates that the spectral type is late-B to early-A. The Ca triplet lines show double-peaked emission and a red-shifted absorption component that reaches below the continuum and could indicate the presence of an accretion flow (Hartmann et al., 1994). We confirm the findings of H97 who classified B268 as a mYSO, based on the presence of CO bandhead emission and the Pa$\delta$ line.

B275 (B7 III): Ochsendorf et al. (2011) performed a detailed spectral classification, resulting in spectral type B7 III. They concluded that B275 is a pre-main-sequence star contracting towards the main sequence. The He i 4009 Å and C ii 4267 Å lines are very weak and when considering the Si ii 4128 Å to Mg ii 4481 Å ratio the spectral type becomes B6-B7. The O i 8446 Å and Ca triplet lines show pronounced double-peaked emission, as do several of the hydrogen lines on top of a photospheric absorption profile. Both first- and second-overtone CO emission is detected. The [O i] 6300 emission line is present. Together with double-peaked emission features and NIR excess, this points to the presence of a rotating circumstellar disk. B275 is a massive YSO.

B289 (O9.7 V): The spectrum includes He ii 4686 and He ii 5411 absorption lines, but the presence of He ii 4200 is hard to confirm given the low signal-to-noise ratio. The He i (+ii) 4026 line is present, He i 4144 and 4387 are weak and He i 4471 and 5876 are strong. In addition, it is possible to identify the C iii 4647/50/51 line complex with a similar strength as the He ii 4686 line, so that the spectral type is O9.7 V. From the $K$-band we obtain a kO9-B1 spectral type. According to H97, B289 might have a NIR excess and be a late O star. The $JHK$ CCD does not display evidence for a NIR excess and there are no emission lines, therefore we cannot conclude that it is a mYSO. The 1.5 and 2.0 $\mu$m photometry by Nielbock et al. (2001) using MANIAC mounted at the ESO La Silla 2.2 m telescope shows IR excess. It is one of the sources surrounded by an IR-bright dusty disk as reported by Chini et al. (2005). The spectrum shows some weak nebular emission lines.

B311 (O8.5 Vz): The He i (+ii) 4026 and He i 4471 lines are almost equally strong and sharp. He i $\lambda\lambda$ 4144, 4387, and He ii $\lambda\lambda$ 4541, 4686 also show up in the spectrum. The C iii 4068/69/70 complex is in absorption together with Si iv 4089 and 4116. We classify this star as O8.5 Vz because the ratio EW(He ii 4686) /EW(He i 4471) $>$ 1.10. From the NIR spectrum, H97 classified it as later than O9-B2 because it has He i in absorption and lacks N iii and He ii. Nielbock et al. (2001) and Chini et al. (2005) detected an IR-bright dusty disk in the $N$ and $Q$-bands. However, Povich et al. (2008) resolved a bow shock associated with it, which is responsible for the NIR emission. B311 is a main-sequence star according to our criteria. The result of the $K$-band classification agrees well with the visual spectral type (kO9-B1). Nebular emission is weakly present.

B331 (late-B): The UVB part of the spectrum is not detected due to the severe interstellar extinction towards this source. The red and near-infrared spectrum is dominated by strong emission lines. In the case of the hydrogen series lines these are superposed on a broad and shallow photospheric profile. According to Hoffmeister et al. (2008), B331 is a B2 V type star based on its visual spectrum. No helium lines are detected in the NIR part of the spectrum, indicating that B331 is a late-B or early-A-type star. It exhibits Br$\gamma$, O i, CO bandhead emission, and several double-peaked emission lines indicating the presence of a rotating circumstellar disk. The SED includes a NIR excess, making it a bonafide mYSO. The spectrum contains some weak nebular emission lines.

B337 (late-B): This object is deeply embedded so that it is not possible to detect the blue part of the spectrum. It does show the Paschen series in absorption with a narrow central nebular emission component. The Ca triplet lines show pronounced double-peaked emission and a blue absorption component remnant of the sky subtraction. The SED presented by H97 is consistent with a B5 V type star when correcting for the distance. No He i lines can be identified in the X-shooter spectrum pointing towards a late-B or early-A-type star. The NIR CCD indicates a NIR excess and the Ca ii triplet lines are in emission and double-peaked; we do not detect CO bandhead emission. B337 is a mYSO.

4 Spectral energy distribution

We construct the spectral energy distribution (SED) of the star by dereddening the X-shooter spectrum, as well as the available photometric data (extending into the mid-infrared wavelength range). We fit the slope of the SED in the photospheric domain (400 – 820 nm) to Castelli & Kurucz models (Kurucz, 1993; Castelli & Kurucz, 2004). As an initial guess, we took the associated $T_{\rm eff}$ and $\log\,g$ from the Castelli & Kurucz666Table 2 in http://www.stsci.edu/hst/observatory/crds/castelli_kurucz_atlas.html model corresponding to the spectral type reported in Sect. 3. We then cross-check our choice of $T_{\rm eff}$ and $\log\,g$ with the results obtained in Sect. 5 and perform the necessary iterations. The final values of $T_{\rm eff}$ and $\log\,g$, corresponding to the best fit Castelli & Kurucz models used, are displayed in Figure 7. We thereby constrain the visual extinction $A_{V}$ and the total to selective extinction $R_{V}=A_{V}/$E$(B-V)$ towards the objects, adopting the parametrisation of Cardelli et al. (1989). From comparing these SEDs to observed spectra a near-infrared excess should become apparent (Figure 7).

The flux calibrated X-shooter spectra and the photometric data points are shown by the black lines and squares, respectively. The photometric data in the $UBVRI$ bands were taken from Chini et al. (1980) who used the 123 cm telescope on Calar Alto, and the 50 cm telescope on the Gamsberg in South West Africa. These measurements have an accuracy of $\pm$0.5 mag. We highlight that the apertures used for this study are quite large, therefore the measurements might be contaminated by close neighbours. For B215 and B331 we obtained the $V$-band magnitude from Hoffmeister et al. (2008), and for B111 from the AAVSO Photometric all sky survey (APASS)777https://www.aavso.org/apass. The $JHK$ magnitudes were taken from the 2MASS survey (Skrutskie et al., 2006) which provides photometry for our objects with an accuracy better than 0.3 mag, with the exception of B215 whose $JHK$ photometry we obtained from Chini & Wargau (1998). The mid-IR magnitudes were measured with the Infrared Array Camera (IRAC) on the Spitzer Space Telescope (Fazio et al., 2004) as part of the Spitzer Legacy Science Program GLIMPSE (Benjamin et al., 2003); the photometric accuracy achieved for our targets is better than 0.2 mag. The $N$ and $Q$-band magnitudes (10.5 and 20 $\mu$m) were obtained (when available) from Nielbock et al. (2001) who observed M17 with the ground-based infrared camera MANIAC. The aperture sizes were 2” for B337, 3” for B331 and B275, 5” for B311, 7” for B215, and 10” for B289; the accuracy in their measurements varies from 10 to 30%. The dereddened spectra and photometry are shown in blue. The dashed grey line represents the Castelli & Kurucz model corresponding to the spectral type of the star. The temperature and $\log{g}$ corresponding to the best Castelli & Kurucz model is indicated in the bottom-right corner; note that this is not identical to the temperature and $\log{g}$ of the star, but it corresponds to the adopted Castelli & Kurucz model. The resulting extinction and stellar parameters are labeled as well.

4.1 Extinction parameters

In order to obtain independent values of the total extinction $A_{V}$ and the total to selective extinction $R_{V}$, we implemented a $\chi^{2}$ fitting algorithm. We first dereddened the X-shooter spectrum varying $R_{V}$ from 2 to 5.5 in steps of 0.1 and then fitted the slope of the Castelli & Kurucz model to that of the dereddened X-shooter spectrum in the photospheric domain to determine the total V-band extinction, $A_{V}$. This allowed us to constrain $R_{V}$ for all sources, except for B163. For this star we lack spectral coverage at $\lambda<850$ nm. As the infrared excess starts at 1000 nm, a too limited spectral range is available for constraining $R_{V}$.

The obtained values for $R_{V}$ range from 3.3 to 4.7 and $A_{V}$ varies from $\sim$6 to $\sim$14 mag (Table 3). H97 observe a similar range, while Hoffmeister et al. (2008) found $R_{V}=3.9\pm 0.2$ for their sample. The latter authors argue that the extinction to M17 is best described by a contribution of foreground extinction ($A_{V}=2$ mag with $R_{V}=3.1$, the average Galactic value) plus an additional contribution produced by local ISM dust. We derive $R_{V}$ and $A_{V}$ for each individual sight-line. The $A_{V}$ values that we obtained agree within the errors with those calculated by Povich et al. (2009) from the total hydrogen column density $N$(H).

The X-shooter spectra include strong diffuse interstellar bands (DIBs). H97 reported that their strength did not vary with $A_{V}$ (or E($B-V$)) as seen in other Galactic sightlines. Massive star forming regions are known to exhibit anomalous extinction properties (Dahlstrom et al., 2013; Oka et al., 2013; Ellerbroek et al., 2013a). We will report the DIB behaviour towards sightlines in M17 in a separate paper.

4.2 Stellar radius

To estimate the radius of the stars we scaled the observed flux to the flux produced by the stellar surface given by the Castelli & Kurucz model, using the distance to M17 ($d=1.98$ kpc). The difference between the measured flux $\rm F_{\lambda}$ and the flux from the model $\rm F_{Kur}$ can be corrected for by multiplying by a factor $(R_{\star}/{d})^{2}$, where $R_{\star}$ is the stellar radius, and $d$ is the distance to the Sun. We calculated the luminosity $\log{L/L_{\sun}}$ using the $V$-band magnitude from Chini et al. (1980), the distance to M17, the $A_{V}$ values from our fit (Sect. 4.1), and the bolometric correction corresponding to the spectral type from Cox (2000).

To assess the uncertainties in the parameters, we calculate the probability $P$ that its $\chi^{2}$ value differs from the best-fit $\chi^{2}$ due to random fluctuations: $P=1-\Gamma(\chi^{2}/2,\nu/2)$, where $\Gamma$ is the incomplete gamma function and $\nu$ the number of degrees of freedom. We normalise the $\chi^{2}$ such that the reduced $\chi^{2}$ corresponding to the best fitting model has a value of unity and we select all models with $P\geq 0.32$ as acceptable fits representing the 68% confidence interval. The finite exploration of the parameter space may result in an underestimation of the confidence interval near the borders of $P(\chi^{2},\nu)$ = 0.32. To avoid an underestimation of the errors we select as boundaries of the 68% confidence interval the first combination of parameters that do not satisfy $P(\chi^{2},\nu)\geq 0.32$ (see Tramper et al., 2014; Ramírez-Agudelo et al., 2017). Like in the case of B275 (Ochsendorf et al., 2011) the obtained stellar radius of the mYSOs does not correspond to the value expected for a main-sequence star. This corroborates our earlier finding that these mYSOs indeed are massive pre-main-sequence stars that are still contracting towards the main sequence. In Sect. 5 we will show that this is consistent with the $\log{g}$ values independently measured from the broadening of the diagnostic lines.

5 Modeling the photospheric spectrum

To further constrain the temperature, luminosity, projected rotational velocity, and surface gravity of the stars we used an automatic fitting algorithm developed by Mokiem et al. (2005) and Ramírez-Agudelo et al. (2017), which compares the H i, He i, He ii, and N absorption lines with model profiles produced by the non-LTE stellar atmosphere model FASTWIND (Puls et al., 2005; Rivero González et al., 2012). For our PMS stars we did not fit the He and N abundances but fixed them to be consistent with the solar values. This is in agreement with what we would expect for such a young stellar population. This method applies the genetic fitting algorithm PIKAIA (Charbonneau, 1995) which allows us to explore the parameter space in an extensive way (Mokiem et al., 2005; Tramper et al., 2014). FASTWIND calculates non-LTE line-blanketed stellar atmospheres and accounts for a spherically symmetric stellar wind. It can be used to examine the dependence of H, He i, and He ii photospheric lines on $T_{\rm eff}$ and $\log{g}$. Its application to B stars has been successfully tested by Lefever et al. (2010).

Inputs to the fit are the absolute $V$-band magnitude ($M_{V}$) and the radial velocity (RV) of the star. We calculated $M_{V}$ using the apparent magnitude reported by Chini et al. (1980) and the extinction coefficient $A_{V}$ found in Sect. 4 adopting a distance of 1.98 kpc. RV was measured following the procedure described by Sana et al. (2013) where we simultaneously adjust the spectral lines for a given object, taking into account variations in the signal-to-noise ratio. The RV values vary from -11 to 20 $\rm km$ $\rm s^{-1}$ and are listed in Tab. 3; for B243, B268, and B289 we list the RV values obtained from the 2013, 2013, and 2012 observations, respectively. We calculated the RV dispersion using the standard deviation of a gaussian distribution ($\sigma_{1D}$) and we find the strikingly low value of $\sim$5 km $\rm s^{-1}$. This is in contrast with the expectation that the massive star population in M17 includes many (close) binaries (e.g. Kouwenhoven et al., 2007; Sana et al., 2012; Dunstall et al., 2015). This points to a lack of short period binaries or a low binary fraction. A quantitative investigation of the low $\sigma_{1D}$ is presented in Sana et al. (2017).

Several atmospheric parameters are obtained from the fitting procedure. The effective temperature ($T_{\rm eff}$) can be constrained using the relative strength of the H i, He i, and He ii lines. The surface gravity ($\log{g}$) is obtained from fitting the wings of the Stark-broadened H i lines. This parameter allows us to constrain the luminosity class (see Mokiem et al., 2005). The projected rotational velocity ($v\sin{i}$) is a natural outcome of the fitting procedure as the models are convolved with a rotational profile to reproduce the observed spectrum. The mass-loss rate is mainly determined by fitting the H$\alpha$ line. Only for B111 ($\log{\dot{M}}=-6.00\pm 0.1$) and B164 ($\log{\dot{M}}=-6.35\pm 3.65$), the two hottest stars in our sample, we obtain a reliable measurement. The parameter describing the rate of acceleration of the outflow ($\beta$) cannot be constrained by our data and was therefore fixed to the theoretical value predicted by Muijres et al. (2012), i.e. $\beta$ = 0.8 for main-sequence stars. Given the low mass-loss rates, the value adopted for $\beta$ does not affect the model results. Another parameter in the fitting procedure is the microturbulent velocity: we allowed this parameter to vary from 5 to 50 $\rm km$ $\rm s^{-1}$, but due to the lack of Si lines in our spectra this parameter remains poorly constrained. This does not affect the determination of the other stellar parameters.

The bolometric luminosity ($\log{L/L_{\odot}}$) is obtained by applying the bolometric correction to the absolute magnitude used as input. The luminosity and effective temperature are used to calculate the radius of the star $R_{\star}$. Using $\log{g}$ we can calculate the spectroscopic mass of the star $M_{\rm spec}=gR^{2}/G$, where $G$ is the gravitational constant. Given the uncertainty in $\log{g}$ it is very difficult to constrain $M_{\rm spec}$ as is evident from Tab. 3.

To calculate the errors in the parameters we follow the procedure described by Tramper et al. (2014) and Ramírez-Agudelo et al. (2017). The best-fitting model is selected based on the $\chi^{2}$ value in the same way as described in Sect. 4.2. The normalisation of the $\chi^{2}$ in this approach is only valid if the best fit is a good representation of the data, which we checked visually for each case (see Appendix B). We select all models with $P(\chi^{2},\nu)\geq 0.05$ as acceptable fits, representing the 95% confidence interval to the fitted parameters. The results from the fitting procedure are listed in Tab. 3; the temperatures obtained with FASTWIND are consistent with the spectral types derived in Sect. 3. The best fitting model, and other acceptable fits (5% significance level or higher models) are shown in Appendix B. The lines used for each fit are displayed in figures 17 to LABEL:P1:fig:FASTWIND_fit_B311; we have given the same weight to all the lines in the fitting procedure.

It is important to point out that the confidence intervals cited in this paper represent the validity of the models as well as the errors on the fit. They do not include the contribution from systematic errors due to the model assumptions, continuum placement biases, etc. For a detailed analysis of the systematic errors and their impact on the parameter determination the reader is referred to Ramírez-Agudelo et al. (2017). As can be seen in Tab. 3 the error bars differ significantly from star to star. Large error bars are obtained in two cases: (i) when the parameter space is poorly explored near the border of the confidence interval and, therefore, the first model that does not satisfy $P(\chi^{2},\nu)=0.05$ is considerably outside the confidence interval; (ii) when the parameter is poorly constrained due to e.g., a low signal-to-noise ratio.

For B275 a good fit of the spectrum could not be obtained while leaving all parameters free. To mitigate this, we first constrained $v\sin\,{i}$ using a fit to only the helium lines. The obtained range of valid values was then used in a fit including the full set of diagnostic lines. H$\beta$ was not included in this fit because of the presence of a very strong 4880 Å DIB blending the red wing of the line. This approach results in an acceptable fit to the spectrum, consistent with the results presented in Ochsendorf et al. (2011); we note that the red wings of He i 6678 and 5875 are not well represented. The morphology of these lines and other He i lines could be an indication that this star is in a binary system; follow-up observations of this source are required to confirm or reject this hypothesis.

5.1 Comparison of the radius estimated by two methods

The radius of the star is estimated in two different ways following the procedures described in the previous sections. In Fig. 8 we show $R_{\star}$ obtained by fitting the SEDs to Castelli & Kurucz models versus the one obtained via the genetic algorithm fitting. The diagonal line represents the one to one correlation and each symbol corresponds to one of our targets. B163, B331, and B337 were left out because for these stars it was not possible to identify any absorption lines, and therefore we did not include them in the GA fitting. The conclusion is that the two methods yield radii that are consistent within the errors.

5.2 Hertzsprung-Russell diagram

Figure 9 shows the theoretical Hertzsprung-Russell diagram (HRD) based on the values of $T_{\rm eff}$ and $L$ obtained in the previous sections. We plotted the PMS tracks from Hosokawa & Omukai (2009) with the ZAMS mass labeled and open symbols indicating lifetimes. We also present the birthline for accretion rates of $10^{-3}$, $10^{-4}$, and $10^{-5}$ $M_{\odot}yr^{-1}$. Assuming that the accretion is constant and that the Hosokawa tracks provide a correct description of the PMS evolution we conclude that 80% of our sample must have experienced an on average high accretion rate (10*-3* - 10*-4* $M_{\odot}$ $\rm yr^{-1}$) and the remaining stars an accretion rate of at least 10*-5* $M_{\odot}$ $\rm yr^{-1}$.

All our confirmed PMS stars (blue circles) are located far away from the ZAMS and their positions can be compared with PMS tracks. The O stars (B111, B164, and B311; grey squares) and the stars without disk signatures (B289 and B215; green triangles) in the X-shooter spectra are located at or near the ZAMS. The position of the B-type star B253 is consistent with it being a PMS star, but its spectrum shows no signatures for the presence of a circumstellar disk.

6 Evidence for the presence of circumstellar disks

In this section we present an overview of the emission-line features thought to be produced by the circumstellar disks. We discuss the nature of the disks of the PMS objects identified in the previous sections and measure the disk rotational velocity, $V_{\rm disk}$, using several hydrogen lines. In some of the objects we identify CO bandhead emission. The infrared excess observed in the SEDs (Sect. 4) provides information on the dust component of the disk.

We observe several double-peaked emission lines and/or CO bandhead emission in the visual to near-infrared spectrum in six of our targets: B163, B243, B268, B275, B331, and B337. A selection of double-peaked lines along the X-shooter wavelength range is shown in Figs. 10 to 14. The fact that we see double-peaked emission is indicative of a rotating circumstellar disk. For all of these objects we observe an infrared excess in their SEDs (Sect. 4). The atypical morphologies observed in the Balmer and Brackett series of B243 and B268 might be an indication of active accretion; nevertheless, further higher resolution observations are needed to confirm this scenario.

6.1 Velocity structure of the gaseous disks

To measure the characteristic projected rotational velocity properties of the gaseous component of the circumstellar disks, we selected a sample of double-peaked emission lines for each of the sources and fitted two Gaussian functions to measure the peak to peak separation. The projected rotational velocity measured in this way corresponds to half of the peak to peak separation. The Gaussian functions fitted to the lines are shown with the solid red lines in Figs. 10 to 14 and the results from these measurements are listed in Tab. 4. The Balmer, Paschen, and Brackett series and Ca ii triplet lines are plotted in each of the columns (from left to right) and lines are labeled in the left part of the plot. The measured peak to peak separation is indicated to the right.

Figure 15 shows the disk projected velocities measured from each of the hydrogen lines against $\rm\log{\lambda f_{lu}}$, where $f_{lu}$ is the oscillator strength. We use $\rm\log{\lambda f_{lu}}$ as a measure of the relative strength of a given line within a line series. The oscillator strengths for the hydrogen lines were obtained from Goldwire (1968). For B163 and B337 we do not have sufficient velocity measurements to draw any conclusions. For four out of the six gaseous disks detected (B243, B268, B275, and B331) we see a clear trend of the projected disk rotational velocity from the hydrogen recombination lines with line strength. This suggests a structured velocity profile of the gaseous disks, in qualitative agreement with the prediction that the high excitation lines form in the inner region (dense and hot) of the disk while the low excitation lines form over a larger area (and on average more slowly rotating part) of the disk.

B268 and B275 also exhibit Ca ii triplet and O i double-peaked emission lines. The velocities measured from these lines are higher than those of the hydrogen lines. This indicates that these lines are formed closer to the star than the hydrogen lines.

Assuming that the disks have a Keplerian velocity structure and adopting the mass of the central object, it is possible to roughly calculate the distance from the star at which the lines are formed. This results in a range from a few hundred to one thousand $R_{\odot}$ (tens to a few hundreds $R_{\star}$).

6.2 CO bandhead emission

CO overtone emission is produced in high-density (1010 – 1011 cm*-3*) and high-temperature (2500 – 5000 K) environments (Charbonneau, 1995). CO is easily dissociated, and therefore must be shielded from the strong UV radiation coming from the star. These conditions are expected in the inner regions of (accretion) disks, which makes the CO bandheads an ideal tool to trace the disk structure around mYSOs (e.g. Bik et al., 2006; Ilee et al., 2013). The shape of the CO lines can be modelled by a circumstellar disk in Keplerian rotation (Bik & Thi, 2004; Blum et al., 2004; Davies et al., 2010; Wheelwright et al., 2010; Ilee et al., 2013). The blue shoulder in the bandheads is a measure of the inclination of the disk: an extended blue shoulder indicates a high inclination angle (i.e., near ”edge-on” view).

We detect CO-bandhead emission in B163, B243, B268, B275, and B331 (see Fig. 6). Although CO-bandhead emission is rare, it is also seen in some B[e] supergiants (Lamers et al., 1998). However, these evolved B[e] stars show numerous forbidden emission lines in their optical spectra, while we only observe the well known [O i] 6300 in some of the sources. Hence, the observed spectra support our hypothesis that the observed disks are related with the accretion process in contrast to a possible origin such as in (more evolved) B[e] stars (see Bik & Thi, 2004). Modelling of the bandheads and the double-peaked hydrogen emission lines will be the subject of a forthcoming paper.

7 Discussion

7.1 Age distribution of the massive PMS population in M17

Age estimates of M17 show that its stellar population is not older than 1 Myr. There is no evidence for a supernova explosion in that area (Broos et al., 2007), consistent with such a young age. We estimated the age of the main-sequence stars by comparing their position in the HRD with Milky Way evolutionary tracks from Brott et al. (2011). We compare the position of our stars with evolutionary tracks that account for the effects of rotation considering a flat distribution of spin rates starting at the measured $v\sin{i}$. To do so, we used the bayesian tool BONNSAI111111The BONNSAI web-service is available at www.astro.uni-bonn.de/stars/bonnsai. from Schneider et al. (2014). Within the uncertainties of the models, the ages obtained for B111, B164, and B311 are 1.60, 0.98, and 0.82 Myr, in fair agreement with the age estimates for the region (Tab. 5). As a way of testing our classification of the PMS stars, we estimated their ages using main sequence tracks through BONNSAI. We obtained main sequence ages ranging from 10 to 90 Myr, which would be inconsistent with M17 being a young region. Assuming that our classified PMS stars are part of M17, we conclude that these stars cannot be post-, but must be pre-main sequence objects.

We estimated the age of the PMS stars by comparing their effective temperatures and luminosities with PMS tracks (Hosokawa & Omukai, 2009). For these stars we find an age span from tens of thousands of years to a few hundred thousand years.

If we compare the position of the (presumable) main-sequence stars B215, B253, and B289 to isochrones, we obtain age estimates of $\sim$9, $\sim$50, and $\sim$4 Myr, respectively. B253 does not present emission lines nor IR excess in its SED, and B215 and B289 do not present emission lines in the X-shooter spectrum but have IR excess longward of 2.3 $\mu$m. If these sources are PMS stars this suggests that their circumstellar disks have already (at least partially) disappeared (see also Sect. 7.2).

Given the size of our sample and the uncertainties in the models we cannot draw a firm conclusion about the age of M17 nor about the possible presence of two distinct populations in this region (Povich et al., 2009). Nevertheless we observe a trend of age with luminosity: the less luminous objects are further away from the ZAMS than their more luminous counterparts. This is in line with the idea (McKee & Tan, 2003; Zinnecker & Yorke, 2007; Davies et al., 2011) that more massive stars form faster and therefore spend less time on the PMS tracks than lower-mass stars.

7.2 Extinction towards the PMS stars in M17

Table 3 lists the extinction properties of our targets. We find that the extinction in the $V$-band varies from $\sim$5 to $\sim$14 mag and that the total to selective extinction is $3.3<R_{V}<5$. The two sources with the highest extinction are situated near or in the irradiated molecular cloud (see Fig.1), but there is no general trend with location in the H  ii region. The sources with higher $A_{V}$ tend to have larger values of $R_{V}$, although the correlation is weak. The overall conclusion is that the extinction is quite patchy, with substantial variation on a spatial scale of 50 arcsec (corresponding to a geometrical scale of 0.5 pc at the distance of M17), similar as to the findings of H97. A dust disk local to the star, for example, may dominate the line-of-sight extinction towards the sources (see Ellerbroek et al., 2013a). Studies of individual sources therefore should not rely on average properties of the region, but should be based on a detailed investigation of the extinction in the line of sight.

7.3 Presence of circumstellar disks around massive PMS stars

We detect signatures of circumstellar disks in six of our sources: B163, B243, B268, B275, B331, and B337. The full sample includes two O-type stars, B111 and B164, that are not classified as potential YSO sources by H97, but as main sequence objects, and do not reveal disk signatures. Two stars, B215 (IRS15) and B289, do not show evidence for gaseous disk material in their spectral lines but do feature excess infrared continuum emission indicative of dust in the circumstellar environment. We thus find that 60% of our massive YSO candidates show clear evidence for circumstellar disks in their spectrum and another 20% likely also feature disks based on NIR excess emission.

In Fig. 9 we have labeled the sources with a gaseous disk detectable in the X-shooter spectrum with a blue dot, the ones with only NIR excess in the SED with a green triangle and the O and B stars with a grey square. The stars with gaseous disks are located further away from the ZAMS than the other sources. However, we cannot link this in a straightforward way to an evolutionary effect. Having noted this, and taking into account that our sources with non-detectable disks are younger than $10^{5}$ yrs, we can conclude that massive stars up to $\sim$20 $M_{\odot}$ retain disks up to less than $10^{5}$ years upon arrival on the ZAMS. We identify the presence of [O i] 6330 emission in three sources (B215, B243, and B275). In Herbig Ae/Be this line is mostly formed in stars with a flaring disk (van der Plas et al., 2008). The spectral profiles of some lines in B268 and B243 might be an indication that this star is actively accreting. For the other sources, as we do not find strong evidence for infall or the presence of jets, we do not know whether these sources are actively accreting or that the observed disks are structures remnant of the formation process.

7.4 Spin properties of the massive PMS stars in M17

The projected rotational velocity ($v\sin{i}$) distribution for 216 O stars in 30 Dor has been published by Ramírez-Agudelo et al. (2013). Their distribution shows a two-component structure consisting of a peak at 80 $\rm km$ $\rm s^{-1}$ with 80% of the stars having $0<v\sin{i}<300$ $\rm km$ $\rm s^{-1}$ and a high velocity tail (containing 20% of the stars) extending up to 600 $\rm km$ $\rm s^{-1}$. Dufton et al. (2013) studied 300 stars spanning spectral types from O9.5 to B3 in 30 Dor; they find a bimodal distribution with 25% of the stars with $0<v\sin{i}<100$ $\rm km$ $\rm s^{-1}$ and a high velocity tail between $200<v\sin{i}<350$ $\rm km$ $\rm s^{-1}$. They estimated $v\sin{i}$ using a Fourier transform method, which allows them to separate the rotational broadening from other broadening mechanisms. Penny & Gies (2009) measured $v\sin{i}$ for 97 OB stars in the Milky Way. They find that 80% of their sample rotate slower than 200 $\rm km$ $\rm s^{-1}$, and that the remaining 20% has $200<v\sin{i}<400$ $\rm km$ $\rm s^{-1}$. Tab. 3 shows the $v\sin{i}$ measured for the stars in our sample. We find that around 30% of our sample is rotating relatively fast (around 200 $\rm km$ $\rm s^{-1}$).

We show the cumulative $v\sin{i}$ distribution functions of these works together with our results in Fig. 16. To quantitatively compare these distributions we performed a Kuiper test, which allows to test the null hypothesis that two observed distributions are drawn from the same parent distribution. The significance level of the Kuiper statistic, $p_{K}$, is a percentage that indicates how similar the compared distributions are. Small values of $p_{K}$ show that our cumulative distribution is significantly different from the one it is compared to. The $p_{K}$ values obtained from comparing our distribution with the ones of Penny & Gies (2009), Ramírez-Agudelo et al. (2013), and Dufton et al. (2013) are $p_{K}=13\%$, 99%, and 18%, respectively. As these values are not lower than 10%, they do not allow us to reject the null hypothesis.

Upon arrival on the ZAMS our sample will span a roughly similar range of masses as the B-star sample in 30 Dor (Dufton et al., 2013). Therefore, it is most appropriate to compare our findings to this sample. According to the models by Brott et al. (2011) and Ekström et al. (2012) the $v\sin{i}$ distribution of stars in our mass range will not change significantly during the first few Myrs of evolution. Dufton et al. (2013) point out that macroturbulent motions only need to be taken into account in cases of relatively slow spinning stars ($v\sin{i}<80$ km $\rm s^{-1}$). We note that Dufton et al. (2013) discuss the possibility of their presumably single star sample to be polluted by (relatively rapidly spinning) post-interaction binaries. We estimated the $v\sin{i}$ that our PMS stars will have upon arrival on the ZAMS using the ZAMS radii corresponding to the end of the PMS tracks from Hosokawa & Omukai (2009). We applied angular momentum conservation and assumed that the stars are rigidly rotating (which is analogous to homologous contraction). When comparing this $v\sin{i}$ distribution (dashed line in Fig. 16) with the B-star sample in 30 Dor we find that $p_{k}=87\%$. Assuming that our stars represent a progenitor population of the LMC B-star sample studied by Dufton et al. (2013) (neglecting metallicity effects and pollution from post-interaction systems in the B stars in 30 Dor sample) we find that the contraction of the PMS stars during the main sequence is consistent with being homologous. Of course, it is premature to firmly state this given the small sample size and the caveats mentioned above.

8 Conclusions

We performed VLT/X-shooter observations of young OB stars in M17. We classified and modelled the photospheric spectra using FASTWIND in order to derive their stellar parameters. We identified the presence of gaseous and dusty disks in some of them based on the emission lines in the spectrum and on the IR excess observed in the spectral energy distribution. We confirm the PMS nature of six objects in this region and conclude that they are on their way of becoming B main sequence stars with masses ranging from 6 to 20 $M_{\odot}$. This constitutes a unique sample of PMS stars that allows us to test theoretical star formation models. Our findings can be summarised as follows:

• •

We confirm the PMS nature of six of the mYSO candidates presented by H97. We conclude that most of our objects must have experienced, on average, high accretion rates.

• •

The age of the PMS objects has been obtained by comparing their position in the HRD with pre-main sequence tracks of Hosokawa & Omukai (2009), whereas the age of the OB stars was estimated using the tool BONNSAI (Schneider et al., 2014) and comparing to tracks of Brott et al. (2011). For the O stars we obtained ages younger than 2 Myr. The pre-main sequence stars have estimated PMS lifetimes of a few hundred thousand years. Given the uncertainty in the age and the fact that we have only a few stars, we cannot conclude that the PMS population corresponds to a second generation of star formation in M17 but found no indication in favour either.

• •

We measured the visual and total to selective extinction towards our objects by fitting Castelli & Kurucz models corresponding to the spectral types. We confirm that the extinction towards M17 is highly variable, as usually observed in star forming regions. We point out that a dust disk local to the star may dominate the line-of-sight extinction towards the sources. Based on the NIR excess (¿ 2 $\mu$m) observed in the spectral energy distributions we found dusty disks in eight of our targets.

• •

Via (double-peaked) emission lines we found evidence for gaseous circumstellar disks in six of our targets, all of which also include a dust component. We measured the projected rotational velocity of the disks from each of the double-peaked lines and found a structured velocity profile among the hydrogen recombination lines for four out of the six disks. For two out of the six disks we were able to identify Ca ii triplet and O i double-peaked emission lines. The velocities measured from these lines are larger than from the ones measured from hydrogen suggesting that they are formed closer to the star than the hydrogen lines.

• •

We measured the projected rotational velocities, $v\sin{i}$, of our stars. About 30% of our sample rotate at around 200 $\rm km$ $\rm s^{-1}$ or faster. The PMS objects are expected to contract and therefore have spun up upon arrival on the main-sequence. Assuming homologous contraction and the absence of processes causing angular momentum loss in their remaining PMS evolution, the $v\sin{i}$ distribution of our sample appears consistent with that of the B stars in 30 Dor once they ignite hydrogen. We note that the contraction in PMS stars is not well understood. Two of the objects will have a $v\sin{i}>300$ $\rm km\,s^{-1}$ upon arrival on the MS.

With this unique, though still small sample we show the potential for constraining models of star formation. A larger sample is needed in order to robustly assess the validity of different theories.

Acknowledgements.

We thank the anonymous referee for for carefully reading the manuscript and many helpful and insightful comments and suggestions. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration, 2013). MCRT is funded by the Nederlandse Onderzoekschool Voor Astronomie (NOVA).

Appendix A Full X-shooter spectra.

Normalised VLT/X-shooter spectra of the objects studied in this paper (Table 1), combining the UVB, VIS and NIR arm. Line identifications are indicated above the spectrum. Note the presence of telluric absorption features centred at 690, 720, 820, 920, 1400, 2010 and 2060 nm. The attenuated region between the H- and the K band (1830 - 1980 nm) has been omitted. The blue part of the spectrum for B163, B331, and B337 is not visible due to severe interstellar extinction. We have included the spectra of B111, B275, and B311 as an example, the version of the appendix for the full sample will be available on A&A.

Appendix B FASTWIND fitting results

In this appendix we show in detail the model fits for our targets and the spectral lines used in each of the cases. We have included the models for B111 and B164 as an example, the version of the appendix for the full sample will be available on A&A.