© The Authors 2024

1. Introduction

Massive stars (with initial mass ≳9 M_⊙; Poelarends et al. 2008) undergo a core-collapse supernova at the end of their life. The formation mechanism(s) of massive stars and the impact of early feedback on their natal environment are generally not well understood. The youngest massive stars are hard to study from an observational perspective. While massive stars are incredibly luminous, they are relatively scarce. Most massive stars are born in young massive clusters that result from the collapse of a giant molecular cloud. The sample of young massive clusters within 2–4 kpc with ages ≲1–2 Myr is small, creating further limitations (see, e.g., Kuhn et al. 2019; Maíz Apellániz et al. 2022a). In the first few million years, the massive stars may not have had sufficient time to disperse the natal molecular cloud through their ionising photons and stellar winds (Geen & de Koter 2022), which causes significant extinction that hampers the observations. As the formation timescale for massive stars is relatively short (several ∼100 kyr), observing the formation process is a challenge.

The outcome of massive star formation is more clear; most of the massive stars end up in binaries or higher order systems (Chini et al. 2012; Kiminki & Kobulnicky 2012; Sana et al. 2012, 2014). The binary orbital period (P_orb) distribution of massive stars of ages 2–4 Myr reveals that ≳70% transfer mass to their companion at some point in their life, with ∼50% of the massive binaries having P_orb shorter than one month (Sana et al. 2012). Ramírez-Tannus et al. (2021) suggest that the P_orb of massive binaries harden over the first 1–2 Myr of their life. This provides potential evidence that massive binaries are not formed with a relatively short P_orb, but rather harden their orbit after star formation. While this trend of massive binaries hardening is relatively clear, the mechanism with which this occurs has so far not been identified.

An important phenomenon in the early life of a young massive cluster is dynamical interaction, where single-binary, binary-binary (or higher-order) systems can exchange angular momentum (Poveda et al. 1967; Leonard & Duncan 1988). In the case of a significantly close passage, two binaries could exchange companions (see e.g. Gualandris et al. 2004) or eject stars at a high velocity, so called runaway stars (Blaauw 1961). This is thought to be most important early on as the young massive cluster should be at its densest during or right after the collapse of the star-forming natal cloud (Clarke & Pringle 1992; Oh & Kroupa 2016). Dynamical interactions are thus closely related to the formation mechanism of stellar clusters.

With the new data releases of Gaia, accurate astrometry and photometry are available for a large fraction of the brighter stars within 2–3 kpc (G ≲ 18 − 21 mag; Gaia Collaboration 2016, 2023; Babusiaux et al. 2023). Young massive clusters and their member stars can now be separated from interloper field stars in greater detail than ever before (see, e.g., Cantat-Gaudin et al. 2018; Castro-Ginard et al. 2022). For an increasing number of nearby young massive clusters in our Galaxy, runaways have been found that are produced by dynamical interactions; these runaways can be traced back to their birth sites (see, e.g., Drew et al. 2021; Maíz Apellániz et al. 2022b; Stoop et al. 2023). These runaways may convey information about the initial conditions of young massive clusters, such as the cluster radius, stellar density, initial mass function and age (Clarke & Pringle 1992; Oh & Kroupa 2016).

M 17 is a nearby giant H II region located in the Sagittarius spiral arm of our Galaxy. The central young massive cluster NGC 6618 is one of the youngest known in the Galaxy, with most age estimates ≲1–2 Myr (Hanson et al. 1997; Hoffmeister et al. 2008; Povich et al. 2009; Ramírez-Tannus et al. 2017). The distance to NGC 6618 and M 17 has been debated in the literature, ranging from ∼1.3 to 2.1 kpc (Hanson et al. 1997; Povich et al. 2007; Hoffmeister et al. 2008). The distance estimate to a spatially co-existing maser source G15.03-0.68 results in $d=1.{98}_{-0.12}^{+0.14}$ kpc (Xu et al. 2011). More recently, Gaia astrometric distance estimates are consistent with a distance of d ∼ 1.6 − 1.7 kpc, implying that the cluster is significantly closer than the maser source (Kuhn et al. 2019, 2021; Maíz Apellániz et al. 2022a). Povich et al. (2009) hypothesise the presence of an older progenitor cluster north-east of NGC 6618, which may explain the coinciding emission bubble and several older OB stars. We follow their nomenclature and refer to this progenitor cluster as NGC 6618 PG, and to the emission bubble as M 17 EB.

NGC 6618 is heavily obscured by dust, with A_V ranging from ∼5 to 15 mag for sources visible in the optical and/or near-infrared (Povich et al. 2009; Ramírez-Tannus et al. 2017, 2018). The central stellar population in NGC 6618 contains (at least) 15 O stars and more than ∼100 B stars (Chini et al. 1980; Hoffmeister et al. 2008), not accounting for possible runaways located outside the cluster. Massive young stellar objects have also been identified which are still on the pre-main-sequence (Ramírez-Tannus et al. 2017). The radial velocities of the OB stars show a relatively low dispersion compared to other young massive clusters (Sana et al. 2017). A low binary fraction is not what one expects for massive stars, which is why the massive binaries in NGC 6618 are hypothesised to have large separations (Sana et al. 2017). One of the implications of this is that dynamical interactions could be even more frequent than initially thought. Wide binaries increase the cross section for dynamical interactions, facilitate transfer of angular momentum, and possibly produce runaways. The same dynamical interactions would cause the wide binaries to harden.

We have studied NGC 6618 and M 17 with Gaia DR3. Section 2 describes the Gaia data processing, membership selection, and the spectral classification of massive stars in NGC 6618 (PG). Section 3 describes the results of the membership selection and key parameters of NGC 6618. In Sects. 4 and 5 we explain how we searched for runaways coming from NGC 6618, and determined some of their physical parameters. We provide a discussion in Sect. 6, and a conclusion and outlook in Sect. 7.

2. Gaia data release 3

We have searched for candidate members of NGC 6618 in Gaia data release 3 (DR3). First, we gathered all Gaia sources in a cone-region with a radius of 0.15 deg. This radius corresponds to ∼4.5 pc at a distance of ∼1.7 kpc and is sufficiently large to include members up to several half-light radii (see Sect. 3). A larger search radius increases the odds of mistakenly including field stars present around NGC 6618. This cone-search region was centred on the brightest and most massive system B189 (also known as CEN 1; Chini et al. 1980; Bumgardner 1992) located in the core of NGC 6618, with (l, b) = (15.0565 deg, –0.6884 deg). This yields at first 3108 sources.

To determine the members of NGC 6618, we applied a set of corrections and filters¹. In Gaia DR3, the correction to the G-band flux and magnitude for sources with 2 and 6-parameter astrometric solutions is already incorporated in the catalogue. We have applied the correction to the parallax to account for the zero-point offset estimated from quasars (Lindegren et al. 2021). To prevent contamination of spurious astrometric solutions, which are typically caused by over-crowding or binaries, we have applied the following filters: first, the renormalised unit weight error (ruwe) should be less than 1.4. If ruwe is larger than 1.4, it likely indicates that the astrometric solution is unreliable (Lindegren 2018). This may exclude binaries and higher order multiples, possibly biasing our membership against multiplicity. Second, the visibility periods used in the astrometric solution (visibility_periods_used) should be 10 or more, which could otherwise indicate astrometric or photometric biases (Gaia Collaboration 2021). Third, we discard sources for which the image parameter determination goodness of fit amplitude (ipd_gof_harmonic_amplitude) is larger than 0.15. Fourth, we exclude sources for which more than one peak (ipd_frac_multi_peak) was identified in more than 10% of the windows used by Gaia. The latter two statistics give indications for crowding or binarity if they exceed their respective threshold (Gaia Collaboration 2021). Last, we remove sources for which more than one source identifier was used in the data processing (duplicated_source), possibly indicating issues in the astrometric solution. After applying these filters, we are left with 2230 sources.

Next, the sources must have had 5 or 6-parameter astrometric solutions, as the proper motion and parallax play a key role in membership selection. Since the distance to NGC 6618 is estimated to be in the range of 1.3–2.1 kpc, we considered stars with parallax (ϖ) smaller than 1.0 mas to be foreground stars. We also required parallax_error (σ_ϖ) to be smaller than 0.12 mas (ϖ/σ_ϖ ≳ 5 at a distance of 1.7 kpc), which allowed to more easily identify whether a source is a member or a field star. This left a total of 498 sources.

While Gaia DR3 now provides non-single-star astrometric models for ∼800 000 sources, the sample of non-single-stars beyond ∼1 kpc is limited (Gaia Collaboration 2023). We have not found multiple star systems consistent with being member of NGC 6618.

2.1. Membership selection using Gaia

We separated the member stars which are part of NGC 6618 from field stars. We have applied PYUPMASK²: the PYTHON port that builds upon the unsupervised membership algorithm UPMASK. We refer to Krone-Martins & Moitinho (2014) and Pera et al. (2021) for a detailed description of the UPMASK algorithm. The ICRS astrometric data and errors are converted to Galactic coordinates with the transformation given in the Gaia documentation³.

We have applied PYUPMASK to the 5D astrometric space consisting of the two Galactic coordinates l and b, the parallax ϖ, and proper motion ${\mu }_{{\mathrm{l}}^{\ast }}$ (${\mu }_{{\mathrm{l}}^{\ast }}\equiv {\mu }_{\mathrm{l}}\cos b$) and μ_b. While radial velocities are now available for nearly 34 million stars in Gaia DR3, for only ∼20% of our 498 sources the radial velocity has been determined (Katz et al. 2023). We have not included the radial velocity as a parameter in PYUPMASK to avoid introducing biases between stars with and without radial velocity.

In the inner loop of PYUPMASK, we used an average of 25 stars per so-called candidate cluster. Candidate clusters are created with the Gaussian mixture model, which was determined by Pera et al. (2021) to be the best performing clustering method. In the outer loop of PYUPMASK, we took uncertainties into account by re-sampling the astrometric parameters using their respective Gaussian uncertainties (${\sigma }_{{\mu }_{{\mathrm{l}}^{\ast }}}$, σ_μb, σ_ϖ), ignoring correlations between these parameters and stars. Uncertainties on l and b are assumed to be negligible compared to those of the proper motion and parallax. We used 10 000 iterations for the outer loop, and applied the Gaussian-uniform mixture model to clean out possible false positives.

The membership probability p for each star is given by the number of times it is assigned member relative to the total number of iterations. We show the membership p distribution in Fig. 1, binned in steps of 0.01. We note that the first two bins between p ∈ [0.00, 0.01] and [0.01, 0.02] contain 250 and 36 members respectively, but have been cut-off for visual clarity. The two distinct peaks in the distribution towards 0.0 and 1.0 show that there is a clear difference between field stars and members. Nevertheless, the field stars heavily outnumber the member stars. Between p of ∼0.2 and 0.8 there are 59 stars for which their membership is less certain. This is a low level of “noise” of on average one star per bin which should be taken into account when deciding a cut-off membership p, to separate the member stars from the field stars. We have decided to constrain this cut-off by having a minimum of five times this “noise level” at five stars per bin (black dashed line in Fig. 1). This resulted in a minimum membership p of 0.96 and yields 47 member stars. The field and member stars are shown in red and blue, respectively.

Fig. 1. Membership probability p distribution for the 498 stars selected in the direction of NGC 6618 with parallax < 1.0 mas. The field stars have p < 0.96, shown in red, while member stars have p ≥ 0.96, shown in blue. The dashed black line (N = 5) shows the threshold used for the membership p cut-off between field and member stars. The first and second bin from the left contain 250 and 36 stars, respectively, and are cut off for better visualisation of the other bins.

As a last step, we removed five stars which deviate by more than 3 standard deviations from the mean for ${\mu }_{{\mathrm{l}}^{\ast }}$, μ_b, ϖ, ending up with 42 members. We list the members in Table 1, sorted by their K_s magnitude and separate them based on whether the stars have been classified on the basis of spectroscopy or not. We have cross-matched the members with the object identifiers in Bumgardner (1992) and Chini et al. (1980), and list their 2MASS object identifier otherwise (Skrutskie et al. 2006). We give a complete table of the members with the Gaia astrometry and photometry in Appendix B to avoid cluttering. Only eight of the 42 members have been spectroscopically classified in the literature or here. One of these, B358, is likely a background post-asymptotic-giant-branch star (Chen et al. 2013; Ramírez-Tannus et al. 2017) not associated with NGC 6618. This star is located in the centre, with a parallax consistent within 1σ with that of NGC 6618 (see Sect. 3). The proper motion of B358 is marginally consistent within 3σ. We have not excluded B358 to maintain consistency in our method. Other than B358, we found no evidence for field-star contamination. The (tentative) spectral types found in the literature all indicate the presence of O and B-type stars.

2.2. Spectroscopic observations

Only for the bright end of the stellar population in the M 17 region spectroscopy is available. Spectra are important to asses the spectral type and determine the radial velocity of the stars. We have obtained spectra for 10 bright candidate OB stars in the M 17 region with either the High Resolution Spectrograph (HRS) mounted on the Southern African Large Telescope (SALT) in Sutherland, the Intermediate-dispersion Spectrograph and Imaging System (ISIS) mounted on the William Herschel Telescope (WHT) at La Palma, and the medium resolution spectrograph X-shooter on the Very Large Telescope (VLT) at Paranal. We observed OB stars with relatively low extinction (WHT), three were part of a monitoring campaign of massive young stellar objects in M 17 (SALT, programme 2022-1-SCI-002), and a few that are part of a multiplicity study in M 17 (VLT). We describe and classify the stars in Appendix A; see Table 2.

2.3. O stars in and around NGC 6618

O stars are the brightest objects in a young stellar population and therefore make up a relatively complete sample. Also, they are potential runaways as the fraction of runaways is claimed to be a strong function of spectral type (de Wit et al. 2005). We have compiled information on all 22 confirmed and candidate O stars in and around NGC 6618 and NGC 6618 PG in Table 2. For 15 of these sources spectra are available that allow to confirm an O-star nature; for six out of these 15 this identification is based on new data presented in Appendix A. For seven of the sources no spectra are available and they are considered O-star candidates on the basis of their K_s magnitude. Dwarf O stars should have M_Ks ≲ −3.0 mag (Martins & Plez 2006; Pecaut & Mamajek 2013). For a distance modulus of 10.6–11.6 (d ∼ 1.3 − 2.1 kpc), this yields K_s ≲ 8.6 mag for O stars in and around NGC 6618. However, the central cluster in NGC 6618 suffers heavily from extinction by dust. Allowing for ∼10 mag of extinction in the V-band, stars would have ∼1 mag of extinction in the K_s-band. We therefore required candidate O stars to have K_s < 10.0 mag. This does not exclude missing out on O stars with A_V ≳ 15 mag which could create a bias, but such stars are too faint for Gaia anyway. Table 2 provides for each O star or O-star candidate their separation from the centre of NGC 6618, distance determined from their individual parallax, ruwe and their member classification (true, false, or candidate). In the case that membership could not be determined due to ruwe > 1.4, we have denoted their member classification as “candidate” (C). Several O stars have ruwe > 1.4 so that we can not determine their distance.

Of the 22 confirmed and candidate O stars, 15 are located near the centre of NGC 6618 (i.e., they are inside the solid circle in Fig. 2). The seven that are located further away are not included in the membership search region. Four sources near the cluster centre have ruwe > 1.4 and are also excluded from the membership search. Of the remaining 11 confirmed and candidate O stars, eight are identified as member. The remaining three (with names B0, B189b, and B311) have deviating parallax or proper motion. Still, these three stars are likely part of NGC 6618 considering their extinction properties and proximity to the cluster centre.

Fig. 2. Position of all identified O stars or O star candidates in the search field centred on NGC 6618 and NGC 6618 PG, save one. The O star candidate 2MASS J18182392-1721517 is located outside of the field shown. Identified member stars are shown with white plusses; those that are O star (candidates) as magenta circles. O star (candidates) that are not identified as member of NGC 6618 are shown as blue circles. For context, we show the DSS2 B, R and I colour image and the Very Large Telescope Survey Telescope OmegaCAM image (ESO/INAF-VST/OmegaCAM). The inset zooms in on the centre of NGC 6618. The green cross in this inset indicates the determined cluster centre. With dashed circles the location of the emission nebulae and open cluster(s) are shown. The bars denoting physical size correspond to ∼0.2 deg (6 pc) and ∼0.03 (1 pc).

The brightest and most massive system B189 is separated into the a and b component since Gaia resolves both. Both B189a and B189b are likely spectroscopic binaries (Hoffmeister et al. 2008), and interferometry also shows that both components could be binaries or triples (Bordier et al. 2022). We have not determined the distance to B189a and B189b as their astrometric solutions are likely spurious.

3. Members of NGC 6618

Studies of membership in Galactic young massive clusters using Gaia typically identify ∼50–500 members (Cantat-Gaudin et al. 2018; Kuhn et al. 2019; Maíz Apellániz et al. 2022a). The exact number depends on several factors, including cluster distance, the amount of line-of-sight and local extinction, the intrinsic number and spatial dispersion of stars in the cluster. We identify 42 members in NGC 6618, which is at the lower end of what is typical. The detectability of stars in NGC 6618 suffers heavily from local extinction by dust in the natal molecular cloud. Previous membership analysis of NGC 6618 from optical and infrared surveys report hundred(s) to thousands of stars belonging to M 17 (Chini et al. 1980; Lada et al. 1991; Hanson et al. 1997; Hoffmeister et al. 2008). Although we identify far fewer members, the quality of Gaia astrometry is such that it allows us to better constrain properties such as the cluster half-light radius, distance, and proper motion. Accurate infrared astrometry that is less affected by extinction, such as the proposed future mission GaiaNIR (Hobbs et al. 2016), could help identify more cluster members.

We show the spatial distribution of the members and all identified O stars in the field in Fig. 2, with the 0.3 deg diameter cone-search region marked by the solid white circle. We have found eight O star members, which we show in magenta. The other 14 O stars, shown in blue, can not be assigned membership due to several reasons. Eleven O stars are located far away from the centre of NGC 6618, several O stars have ruwe > 1.4 making their astrometry spurious, and several O stars have a significantly deviating proper motion from the cluster average.

The 42 members, shown with the black and white plus symbols, display significant spatial clustering. We show the central 0.1 × 0.1 deg of NGC 6618 in the inset in the top-right panel in Fig. 2. The members are all concentrated in this inner panel of the cone-search region. Several factors could be responsible for this. First, since NGC 6618 is young, the members could still be highly concentrated. Second, extreme extinction in the outer regions of NGC 6618 could obscure our view, causing us to only see the central region. Third, the membership algorithm favours stars in higher stellar density regions, which could also create a bias.

We summarise the astrometric, kinematic and physical properties of NGC 6618 in Table 3. To determine the centre of the cluster, we perform Monte Carlo simulations and bootstrapping. We randomly select 80% of the members (without replacement) and determine for this bootstrapped sample the mean l and b. This is repeated 10⁶ times to obtain a distribution of means in l and b. The centre is determined as the 50th percentile, and the 1σ uncertainty as the 16th and 84th percentile of the distribution of means. This results in (l_NGC 6618, b_NGC 6618) = (15.058 ± 0.001 deg, –0.686 ± 0.001 deg). We show the position of the centre with the green cross in the inset in the top-right panel in Fig. 2, close to the O4 V + O4 V central multiple system B189. That B189 (CEN 1) is exactly in the centre of NGC 6618 is therefore a good assumption. Yanza et al. (2022) determine the centre from compact radio sources in NGC 6618. They find (l, b) = (15.050 deg, –0.689 deg), which is in reasonable agreement considering the aforementioned biases.

The radius of NGC 6618 can be defined in different ways. Here we determine the 2D projected half-light radius r_hl, 2D. To do this, we weigh each member by its relative K_s-flux, which should be least affected by extinction. This results in r_hl, 2D ∼ 0.006 deg. Another way to calculate a radius is to take the standard deviation of the 2D projected radii to the centre of the members, which gives ∼0.011 deg. At the distance of NGC 6618 (1674 pc, see below), these two radii translate to ∼0.18 and 0.32 pc respectively. If we assume a Plummer density distribution, the 3D r_hl will be given by 1.3r_hl, 2D (Plummer 1911), which is 0.23 and 0.42 pc for our two radius definitions, respectively.

Yanza et al. (2022) determine a 2D radius of 0.014–0.018 deg (0.4–0.7 pc) from the two standard deviations on a 2D ellipsoidal Gaussian. This results in a 3D radius of 0.52–0.91 pc assuming again a Plummer density distribution. Our r_hl is consistent within a factor of two, a result that could be subject to observational biases due to extinction and membership selection.

We mention that the compact radio sources in Yanza et al. (2022) are likely low mass stars, while our members have higher masses, are brighter, and have higher quality Gaia astrometry and photometry. The difference between our determined radius and that of Yanza et al. (2022) could therefore be due to an intrinsically mass dependent spatial distribution of stars. Massive stars could have preferentially formed in or near the centre while lower mass stars could originate from further out regions. A conservative range for the radius of NGC 6618 r_NGC 6618 is therefore between 0.2–1.0 pc.

We show the proper motion of the members and O(B) stars in the left panel in Fig. 3, adopting similar markers and colours as before. The grey dots show the 456 determined field stars. A large percentage of the O stars with ruwe < 1.4 can be seen to have a deviating proper motion with respect to the members. At the distance of NGC 6618 (see below), this would correspond to a transverse velocity difference of ≳10 km s⁻¹. We investigate whether this is due to their runaway nature in Sect. 4. The proper motion of the members is similar compared to other young massive clusters (1–3 km s⁻¹; Kuhn et al. 2019). A precise determination of the cluster proper motion is referred to Sect. 4.

Fig. 3. Proper motion and parallax for members, all identified O-type stars, and field stars in the search field centred on NGC 6618. The members and O-type stars are coloured and marked similarly as in Fig. 2. O-type stars with ruwe > 1.4 are not shown in both panels. Left: Proper motion distribution in the Galactic coordinate frame. We show for context the transverse velocity difference equal to ∼1 mas yr−1 proper motion difference at the distance of NGC 6618 (1674 pc). The cluster proper motion determined with the mean of the members and from the runaways (see Sect. 4) are shown with the green and red cross respectively. Right: Parallax distribution as a function of the G-magnitude. The best-fit parallax and 1σ uncertainty of NGC 6618 determined from the members are shown with the grey dashed line and bar respectively.

The proper motion dispersion of NGC 6618 in the l and b direction is σ_{μl^*, NGC 6618} = 0.26 mas yr⁻¹ and σ_{μb, NGC 6618} = 0.44 mas yr⁻¹, respectively. At the distance of NGC 6618, this translates into σ_{l, NGC 6618} = 2.1 km s⁻¹ and σ_{b, NGC 6618} = 3.5 km s⁻¹, respectively. The radial velocity dispersion σ_R of NGC 6618 is 5.5 ± 0.5 km s⁻¹ (Sana et al. 2017; Ramírez-Tannus et al. 2021) and differs by a factor of about two with σ_l and σ_b. A larger σ_R may indicate the presence of binaries, which contribute to the radial velocity dispersion due to orbital motion. Considering that the three velocity dispersions differ by a factor of about two, it is unclear whether these are truly different or are introduced by small number statistics (42 members considered here; 12 stars in Ramírez-Tannus et al. 2017).

We show the parallax of the members and O stars against their G-magnitude in the right panel of Fig. 3, again adopting similar markers and colours as before. We determine the distance to NGC 6618 by setting up a log-likelihood function (without priors) similar to Cantat-Gaudin et al. (2018) and Bailer-Jones et al. (2021) with

$$\begin{array}{cc}\hfill P(d|\varpi ,{\sigma }_{\varpi })& ={\displaystyle \prod _{i=1}P({\varpi }_i|d,{\sigma }_{{\varpi }_i})}\hfill \\ \hfill & ={\displaystyle \prod _{i=1}\frac{1}{\sqrt{2\pi {\sigma }_{{\varpi }_i}^2}}\exp (-\frac{{({\varpi }_i-\frac{1}{d})}^2}{2{\sigma }_{{\varpi }_i}^2}),}\hfill \end{array}$$(1)

where P is the unnormalised probability distribution for a given ϖ and assumed Gaussian uncertainty σ_ϖ. The best-fit distance is given by the mode and the 1σ uncertainties by the 16th and 84th percentile of the probability distribution, respectively. We find a distance d_NGC 6618 = 1674${}_{-18}^{+19}$ pc, with a parallax ϖ_M 17 = 0.5974 ± 0.0065 mas. The 1σ uncertainty on the parallax range is depicted by the grey bar in Fig. 3. Most of the members are consistent within 2σ with this distance. The north-east component of the visual binary B189 (G ∼ 12 mag) shows a deviating parallax with ϖ = 0.472 ± 0.049 mas, which could be attributed to a spurious astrometric solution due to its multiple nature (not including this star yields a distance of 1667 pc).

The parallax of most of the non-member O stars is consistent within 2σ with the determined parallax of NGC 6618. The clear two outlier O stars, which deviate more than 5σ are BD–16 4831 (∼2.4 kpc) and B311 (∼1.4 kpc). Since BD–16 4831 is in this paper determined to be an O9.7 Ia, this star is significantly more evolved than any age estimate for NGC 6618. BD–16 4831 is not located in NGC 6618 and instead is supposedly one of the main ionising sources in NGC 6618 PG (Povich et al. 2009). It is unlikely that BD–16 4831 is associated with NGC 6618 (PG) and we assume this star to be a background star. The O8.5 Vz star B311 (d ∼ 1.4 kpc) is positioned in the centre of NGC 6618 and has a radial velocity of 4.2 ± 0.4 km s⁻¹ (Ramírez-Tannus et al. 2017), such that it is unlikely to be a radial runaway. B311 suffers from an extinction A_V ∼ 6 mag. With only 2 magnitudes of foreground A_V, B311 is likely located inside NGC 6618 (Hoffmeister et al. 2008). One might attribute its deviating distance to a spurious astrometric solution, but we find no evidence for this in the Gaia astrometry and goodness-of-fit indicators. More investigation is needed to determine whether B311 is a member or foreground star.

The colour-absolute magnitude diagram (CAMD) gives us more insight in the age and extinction properties of the members and O stars. We show the CAMD from the Gaia photometry in Fig. 4, propagating the individual uncertainties on the distance in the M_G value, and adopting similar colouring and marking as before. The reddening line of an O9-9.5 V star (R_V = 3.1) is shown with the black dashed line and red crosses superimposed on this line denote A_V from 2.0 to 14.0 mag in steps of 2.0 mag. The O stars can be seen to follow a similar pattern as the reddening line. The least extincted O stars located in the top-left thus have around A_V ∼ 2.0 mag, while the most extincted O stars have A_V ≳ 10 mag. Figure 4 shows that (variable) extinction is a major issue in the nebulous region M 17. The R_V also ranges from 3 to 5 within the H II region (Ramírez-Tannus et al. 2018).

Fig. 4. Colour – absolute magnitude diagram for members, O-type stars, and field stars in the search field centred on NGC 6618 (PG). The members and O stars are coloured and marked similarly as in Fig. 2. Absolute magnitudes are determined from individual parallaxes. Shown are three PARSEC isochrones with AV = 2.0 mag and age of 0.5, 1.0 and 2.0 Myr with the solid, dashed and dot-dashed tracks, respectively. We give the reddening line (RV = 3.1) for an O9–9.5 V star with the dashed black line and superimpose on this line the position of 2–14 mag of AV, in steps of 2 mag, with the red crosses.

We include isochrones in Fig. 4 from the Padova and Trieste Stellar Evolution Code PARSEC (v1.2S) + COLIBRI (S_37 + S_35 + PR16) models (Bressan et al. 2012; Chen et al. 2014, 2015; Tang et al. 2014; Marigo et al. 2017; Pastorelli et al. 2019, 2020). Shown are three isochrones, each with A_V = 2.0 mag, but with different ages. Current age estimates for NGC 6618 are typically ∼1.0 Myr or younger, with upper limits around 2.0 Myr (Hanson et al. 1997; Hoffmeister et al. 2008; Ramírez-Tannus et al. 2017). The displayed isochrones have ages of 0.5, 1.0, and 2.0 Myr shown with the solid, dashed, and dash-dotted yellow lines respectively. These isochrones display a vertical main-sequence in the top-left above M_G ∼ 0 mag. Towards fainter and redder magnitudes, they transition to the pre-main-sequence.

Most of the members are to the red of three isochrones. This suggests that these stars have A_V ≳ 2 mag, similar to the O stars. Around seven members are located below G_Bp – G_Rp ≲ 2 mag and are fainter than the isochrones. We deem these members more likely to be false positive field stars, rather than members with A_V below 2 mag or ages older than 2.0 Myr. As described in Sect. 2.1, we can still expect several false positive field stars.

The members and available isochrones allow us in theory to estimate the age of NGC 6618 by finding the best-fit isochrone (see e.g. Jørgensen & Lindegren 2005). However, for NGC 6618 we can clearly see that the spread of the members, caused by strong variable extinction, complicates this method greatly. To solve this, we would need to correct for the variable extinction by obtaining their intrinsic magnitudes and colours through their spectral types. Unfortunately, the spectral types are only available for the brighter and thus more massive stars in NGC 6618. De-reddening these stars places them on the main-sequence, which gives no information on the age (the main-sequences of the isochrones nearly perfectly overlap in Fig. 4).

4. Dynamically ejected runaways

We have searched for runaways originating from NGC 6618. As NGC 6618 is estimated to be younger than 2 Myr, runaways are most likely produced by dynamical interactions rather than the supernova mechanism (e.g. van der Meij et al. 2021). The dynamically ejected runaways are expected to be produced preferentially in the densest region in a cluster (Fujii & Portegies Zwart 2011).

We have searched for runaways within a 5 deg circle centred on NGC 6618. At a distance of ∼1.7 kpc, we are able to find a runaway with v_T = 100 km s⁻¹ up to 1.5 Myr ago for this 5 deg radius. These stars are subject to the same filters and corrections as described in Sect. 2. ensure that the found runaways come from NGC 6618 and are not foreground or background stars, the parallax of the runaways should be consistent within 3σ with that of NGC 6618. On top of this, runaways should have a fractional parallax uncertainty ϖ/σ_ϖ > 10. This stricter requirement on the accuracy in distance than the cut on the parallax to determine membership (see Sect. 2.1) is to avoid false positives. It also indirectly ensures that the proper motion and thus trace-back of the runaways is accurately known, since these astrometric parameters are correlated.

We have adopted a cut-off magnitude at K_s = 11 mag, where we have cross-matched the Gaia catalogue with the 2MASS catalogue. This allows us to find not only O-type runaways, but also B-type runaways. The absolute magnitude for each runaway is calculated with their individual parallaxes. For fainter magnitudes, it becomes increasingly more difficult to distinguish between runaways and interloper field stars. The effect of extinction is less by a factor of ∼10 in the K_s-band than in the G-band, facilitating the discovery of heavily extincted runaways (provided that they have accurate and reliable astrometry).

The runaways should be escaping or have escaped from NGC 6618. Runaways should therefore be moving away from NGC 6618 and not towards NGC 6618, that is, the relative proper motion vector is pointing away from the centre. The velocity dispersion ranges from 2 to 4 km s⁻¹, so we require runaways to have Δv_T > 5 km s⁻¹ to not pick up on members of NGC 6618. We note that typical 3D escape velocities are 3–5 km s⁻¹ for young massive clusters with masses of ∼10⁴ M_⊙, while we have only considered 2D velocities. The aforementioned constraints result in a total of 2929 stars.

The majority of these stars will move in a “random” direction and only few are actually runaways coming from the cluster. To identify the latter, we have traced back in time the stars to investigate whether they were in NGC 6618 in the past, using

$$\begin{array}{cc}\hfill l_{\mathrm{sep}}(t)& =(l+\frac{t·{\mu }_{l^{\ast },}}{3.6\times {10}^6·\cos (b)})\hfill \\ \hfill & -(l_{\mathrm{NGC}6618}+\frac{t·{\mu }_{l^{\ast },\mathrm{NGC}6618}}{3.6\times {10}^6·\cos (b_{\mathrm{NGC}6618})})\mathrm{deg},\hfill \\ \hfill b_{\mathrm{sep}}(t)& =(b+\frac{t·{\mu }_b}{3.6\times {10}^6})\hfill \\ \hfill & -(b_{\mathrm{NGC}6618}+\frac{t·{\mu }_{b,\mathrm{NGC}6618}}{3.6\times {10}^6})\mathrm{deg},\hfill \end{array}$$

where t is the time in years, and l_sep and b_sep are the separations between a star and the centre of NGC 6618. Since NGC 6618 is located in the Galactic plane, b will likely be between –1 to 0 deg such that the impact of the cos(b) and cos(b_NGC 6618) factors is relatively small. We also do not account for the gravitational potential of our Galaxy, which is a good approximation up until a few Myr (Hoogerwerf et al. 2001).

We expect that true runaways have come from the centre of NGC 6618. The minimum separation between the runaway candidates and the centre of NGC 6618 should therefore be relatively small, taking into account the uncertainties in proper motions. With a previously determined r_NGC 6618 between 0.2–1.0 pc, we require that the runaways come within a conservative separation of 1.0 pc from the centre of NGC 6618. The runaway candidates are traced back in time for a maximum of 1.5 Myr considering the current age estimates of M 17 (≲1 Myr; Ramírez-Tannus et al. 2017).

While searching for these runaways we found that the minimum separation of the runaway candidates with respect to the cluster centre becomes larger (specifically in b_sep) the longer ago (in time) they were ejected. We attribute this to a marginally different cluster proper motion than what can be determined with the limited sample of 42 members.

To alleviate this issue, we leave both μ_{l^*, NGC 6618} and μ_{b, NGC 6618} as free parameters and adopt an iterative approach. Every iteration, we find runaways satisfying the above conditions for initial cluster proper motions. With the found runaways, we determine new cluster proper motion which minimises the total separation of the found runaways weighted by the number of runaways found. With the new cluster proper motions, we repeat the process of finding runaways, until the cluster proper motion and separation of the found runaways reached a stable minimum. In each iteration, we filter out bright late-type stars (red giants) described below.

The procedure could in theory converge to a non-physical solution, where the separation of only one or two runaways would be minimised. In practice, the first iteration finds 11 runaways and the second (and final) iteration finds 13 runaways with the re-determined cluster proper motion. After this, the runaways, cluster proper motion, and minimum separation were optimised. The re-determined cluster proper motion is μ_{l^*, NGC 6618} = −1.33 mas yr⁻¹ and μ_{b, NGC 6618} = −0.96 mas yr⁻¹. This differs by ∼0.07 and 0.1 mas yr⁻¹, respectively, from the mean proper motion of the members. At the distance of NGC 6618, this would correspond to a velocity difference of ∼1 km s⁻¹, or 1 pc Myr⁻¹, which is significant compared to the radius of ≲1.0 pc. While we could have increased our adopted r_NGC 6618, this would cause us to wrongly determine important physical parameters of the runaways such as their impact parameter and kinematic age. A different cluster proper motion causes the (time of) minimum separation to change. We show the proper motion determined with the runaways and with the mean of the members with the green and red cross respectively in the left panel of Fig. 3. We adopt the proper motion determined here with the runaways.

The found runaways may contain false positives, that is, interloper field stars that happen to satisfy the conditions listed above. Bright late-type interlopers, such as red giants, are easily identified in colour–colour diagrams, where the effective temperature of stars determines their location. Specifically, we made use of the (J − H) – (H − K_s) and the (J − H) – (G_BP − K_s) diagrams, which we show in the left and right panels of Fig. 5, respectively.

Fig. 5. Colour–colour diagrams for the O stars (blue), late-type stars (red) and early-type runaways (green). Left: 2MASS colour–colour diagram. We indicate the main-sequence with the solid grey line and highlight the position of an M0, A0, and B2 or earlier main sequence star with a grey label and dot. Reddening lines for an M0 V and a B2 V star are plotted assuming RV = 3.1. Right: Gaia and 2MASS colour–colour diagram.

We show the main-sequence track in both panels with the solid grey line, where the early and late-type stars occupy a different position. The O and early-type B stars all have the same intrinsic colours in these diagrams, which we indicate with < B2 V. The black dashed lines show the reddening lines for a B2 V and M0 V star (R_V = 3.1). Even accounting for reddening, early and late-type stars should be separated from each other. For context, we show the O stars in Table 2 in blue, which mostly follow the B2 V reddening line. The O star B197 is located on the M0 V reddening line in the left panel, which could indicate it is a red giant, since it has only been classified based on photometry. The runaways consistent with being early and late-type stars are shown in green and red, respectively. The six identified late-type stars have colours consistent with late-type stars.

On top of this, we have investigated the T_eff and log(g) estimated with the General Stellar Parameterizer from Photometry (GSP-Phot) in Gaia DR3 (Bailer-Jones 2010, 2011; Bailer-Jones et al. 2013; Andrae et al. 2023). GSP-Phot infers stellar parameters from the G_BP and G_RP spectra. Three of the six identified late-type stars have T_eff (teff_gspphot in Gaia DR3) in the range of 4500–5500 K, while T_eff for the early-type stars are almost all not available. Similarly, the log(g) (logg_gspphot in Gaia DR3) of these three late-type stars are in the range of 1.7–2.7, while those for the early-type stars are again not available. We identify the six late-type stars as false positives and have removed them, leaving us with 13 runaway stars in our final runaway search iteration.

5. Physical properties of runaways

We list the spectral and kinematic properties of the runaways in Table 4. Seven of these runaways are either spectroscopically confirmed or photometrically consistent with being an O star and are also listed in Table 2. We visualise the position and tangential velocity of the runaways in Fig. 6. The arrows indicate the runaway’s motion relative to NGC 6618, with the length of the arrow proportional to the absolute proper motion. The runaways consistent with being an O star are shown in blue, and the likely B or later-type runaways in red. Six to seven of the thirteen found runaways move in the direction of NGC 6618 PG. Several of these runaway stars were identified by Povich et al. (2009) to be part of NGC 6618 PG. We discuss this in Sect. 6.6.

Fig. 6. Position and motion of the found runaways relative to NGC 6618, with the DSS2 B, R, and I colour image and the Very Large Telescope Survey Telescope OmegaCAM image (ESO/INAF-VST/OmegaCAM). We show the proper motion of the runaways relative to NGC 6618 with the arrows, which are proportional to their absolute proper motion. Runaways which are consistent with being O stars in Table 2 are coloured in blue, and coloured red otherwise. We show the position of NGC 6618 (PG) and M 17 (EB) using white circles. The angular size equivalent to 10 pc is given at the distance of NGC 6618 (1674 pc).

All except one runaway are separated by ∼0.1–0.4 deg from NGC 6618, corresponding to ∼3–12 pc. Only one runaway, 2MASS J18182392-1721517, is separated by ∼1.25 deg (35 pc). There are several reports in the literature of excess infrared emission at the position of this runaway, possibly indicating an H II region or “bubble” (Simpson et al. 2012; Su et al. 2018). We show in Fig. 7 the Spitzer Space Telescope Multi-Band Imaging Photometer (MIPS) field around 2MASS J18182392-1721517 (24 μm; Werner et al. 2004; Rieke et al. 2004). We can see that the infrared emission appears to have a bowshock-like shape. The blue arrow indicates the proper motion of this runaway relative to NGC 6618, from which we can see that the bowshock shape is roughly in the direction of motion. We further discuss the nature of 2MASS J18182392-1721517 and the bowshock in Sect. 6.7.

Fig. 7. Spitzer 24 μm Multi-Band Imaging Photometer grey-scale image of the bowshock around the runaway star 2MASS J18182392-1721517. The blue square shows the position of the star. We show the relative proper motion with respect to NGC 6618 with the blue arrow.

We determine the transverse velocity with respect to NGC 6618 (|Δv_T|) for each runaway from their individual proper motions and parallaxes, propagating their uncertainties. We list |Δv_T| in Table 4. Almost all runaways have |Δv_T| ranging from ∼10 to 20 km s⁻¹, with only 2MASS J18182392-1721517 having a significantly higher transverse velocity (∼65 km s⁻¹). Gaia provides radial velocities for two runaways, and no other radial velocities are known in the literature. The stars in NGC 6618 have a mean radial velocity of ∼6.4 km s⁻¹ (Ramírez-Tannus et al. 2017), indicating that both runaways have radial velocities that deviate from this by more than 3σ. This could either be due to their runaway nature or to binarity (or both).

We next focus on the kinematic age (t_kin) and the minimum separation between the runaway and cluster centre (r_imp, 2D, the projected impact parameter) of the runaways, which require having to deal with several uncertainties. The exact origin of each runaway in NGC 6618 is unknown. We therefore assume an uncertainty region around the centre of NGC 6618 with a radius of ∼1.0 pc, determined to be the upper limit on the cluster radius of NGC 6618 in Sect. 3. The proper motion of NGC 6618 has been determined in Sect. 4 and we assume 1σ uncertainty of 0.02 mas yr⁻¹ in both μ_l^* and μ_b. The uncertainty are similar to that of the runaways, which have a 1σ uncertainty in μ_l^* and μ_b in the range of 0.01–0.03 mas yr⁻¹.

We resort to MC simulations to determine both t_kin and r_imp, 2D simultaneously. In each iteration, we randomly draw the runaway μ_l^* and μ_b from a 2D Gaussian distribution with mean, standard deviation and covariance matrix equal to the observed value, 1σ uncertainty and correlation. The proper motion of NGC 6618 is also randomly drawn from a 2D Gaussian distribution similar to the runaways, with the uncertainty listed above. In each iteration, we determine the time when the runaway “enters” the uncertainty region around NGC 6618, is closest to the centre, and “exits” the uncertainty region of NGC 6618. After 10 000 iterations, we obtain three distributions of time from which we calculate the lower limit, best-fit, and upper limit on t_kin as the 16th percentile, mean, and 84th percentile, respectively. We determine the best-fit r_imp, 2D and its uncertainty as the mean and standard deviation on the distribution of closest approach to the centre.

Table 4 lists the determined t_kin and r_imp, 2D for the runaways. The uncertainty in t_kin is not symmetrical. Since the uncertainties in the motion of the runaways and cluster increase with time, the positive error will be larger than the negative error when calculating further back in time. We have taken the absolute values to avoid confusion. We note that r_imp, 2D is a projected distance and that we have minimised the total r_imp, 2D in Sect. 4 by having the cluster proper motion as a free parameter.

We show |Δv_T| as a function of t_kin in Fig. 8, adopting similar colouring and marking as in Fig. 6. The t_kin of the runaways ranges from ∼100 to 1250 kyr, concentrating between 200–700 kyr. The already identified O stars are shown in blue and have t_kin less than ∼600 kyr. The latest O star runaway could have been ejected as recently as 70–170 kyr ago. We show the adopted minimum |Δv_T| and maximum t_kin with the dashed black lines. This shows that the typical runaway |Δv_T| between 10–20 km s⁻¹ is not a result of our cut-off. We would have found runaways with |Δv_T| between 5–10 km s⁻¹ if they existed. Similarly, we do not find runaways with t_kin ≳ 1200 kyr, making it unlikely that there is a significant runaway population older than this.

Fig. 8. Relative transverse velocity (|ΔvT|) of the runaways as a function of their kinematic age (tkin). The identified O stars in Table 2 are shown in blue; red is used for B or later-type stars. The cut-off |ΔvT| and tkin used in the runaway search are indicated with the dashed black lines.

Figure 9 shows the absolute K_s magnitudes, determined from their individual parallaxes, adopting similar colours as before. The already spectroscopically confirmed O stars can be seen to be brighter than M_Ks < −1.5 mag, as expected. The spectroscopically identified B0 V + B1 V system is the brightest B or later-type system, with M_Ks ∼ −2.2 mag. We are likely missing runaways fainter than M_Ks ≳ 0 mag (K_s > 11 mag), indicated with the dashed black line. As previously mentioned, this cut-off is intentional as for fainter magnitudes it becomes increasingly more difficult to distinguish between false positives and true runaways. Figure 9 shows how the adopted magnitude cut-off significantly impacts the found runaways. If we adopted M_Ks < −1 mag, we would have only found runaways with t_kin < 600 kyr.

Fig. 9. Absolute MKs magnitude of the runaways as a function of their kinematic age (tkin), marked and coloured similarly as in Fig. 8. The cut-off MKs and tkin used in the runaway search are indicated with the dashed black lines.

Last, Fig. 10 shows r_imp, 2D as a function of t_kin. We have converted this quantity to the physical projected separation at the distance of NGC 6618 assuming a negligible radial motion. Almost all runaways are consistent with being ejected within 0.2–0.4 pc of the centre of NGC 6618. The identified O stars could all have been ejected from within ∼0.25 pc. Two outliers can be seen which are determined to have r_imp, 2D between 0.5–1.0 pc. These two runaways are also ejected relatively long ago compared to the other runaways and are also relatively faint with M_Ks between –0.5 and –1.0. It could be possible that these were ejected from a different part of the cluster or represent false positives, and we have searched too far back in time and at too faint M_Ks magnitudes. Nevertheless, we have included them for completeness.

Fig. 10. Projected impact parameter (rimp, 2D) of the runaways as a function of their kinematic age (tkin), marked and coloured similarly as in Fig. 8. The cut-off rimp, 2D and tkin used in the runaway search are indicated with the dashed black lines.

6. Discussion

6.1. Distance to M 17

The distance to M 17 has been debated in the literature for decades, with estimates ranging from 1.3 to 2.9 kpc (Ogura & Ishida 1976; Hanson et al. 1997; Russeil 2003). More recent distance estimates resulted in ∼1.6 and 2.1 kpc (Povich et al. 2007; Hoffmeister et al. 2008). The distance of 1.98${}_{-0.12}^{+0.14}$ kpc by Xu et al. (2011) to a maser associated with M 17 lies in between the latter two distance estimates and settled these discrepancies. We note, however, that Xu et al. (2011) fit three astrometric parameters (parallax and proper motions) to only four data points, leaving room for observational biases.

With Gaia, the astrometric distance to M 17 is closer to ∼1.7 kpc, accounting for the parallax zero-point offset. Kuhn et al. (2019) estimates a distance of 1.68${}_{-0.11}^{+0.13}$ kpc with DR2, Maíz Apellániz et al. (2022a) recover 1.70${}_{-0.039}^{+0.041}$ kpc with EDR3 and we find 1.674${}_{-0.018}^{+0.019}$ kpc with DR3. Kuhn et al. (2021) estimate a distance of 1.54${}_{-0.20}^{+0.16}$ kpc, from a sample of 11 member stars. The distance determined by Xu et al. (2011) is consistent within 2σ with a distance of 1.75 kpc, not including the aforementioned uncertainties.

The recent Gaia based estimates for the distance of NGC 6618 differ by a factor of ∼0.85 from the previously used maser estimate (see e.g. Ramírez-Tannus et al. 2017). This has implications for luminosity and mass estimates of young stellar objects in NGC 6618. Masses estimated from spectral-energy-distribution fitting and spectroscopy scale with d², which results in mass estimates for these objects decreasing on average by a factor of ∼0.71.

6.2. Runaway properties

We have identified 13 runaway systems consistent with being ejected from the centre of NGC 6618. BD–16 4826 is a known binary system of which the primary is an O5 V((f))z star and the secondary could be either an O9 or B0 V star. Six to seven of these systems are likely O-type stars, while six to seven others are likely to be B-type stars. We note several observational biases. We have introduced a brightness cut-off at M_Ks < 0 mag, which is why no runaways were found with A or later spectral types. Similarly, we could miss out on heavily extincted, therefore fainter runaways either due to such systems having large uncertainties in the Gaia astrometry or simply because they are not detected at all.

We can estimate the percentage of massive runaways in two ways. First, we compare the number of O-type runaways to the total number of O stars originating from NGC 6618. We have six O-type runaways and 15 O-type stars located in NGC 6618, which includes the runaways and excludes the background O9.7 Ia star. The O star runaway percentage is therefore 29%. Not only could we miss O-type runaways due to the aforementioned biases, we could also miss O-type stars in NGC 6618 because of the extreme extinction in the cluster itself. The 29% of O-type runaways is therefore still somewhat uncertain.

Second, we use a brightness cut-off to estimate the percentage of massive runaway stars. The detectability of massive stars in NGC 6618 diminishes for K_s-magnitudes fainter than ≳9.0 mag. We have seven runaways and 22 stars in NGC 6618 brighter than this. This results in a similar runaway percentage of 32% for stars brighter than K_s ≲ 9.0 mag. We have included two systems for which their origin is unclear (BD-16 4822 and LS 4972) as part of NGC 6618.

Since these runaways originate from a cluster with an age estimated to be ≲1 Myr (Hanson et al. 1997; Hoffmeister et al. 2008; Povich et al. 2009; Ramírez-Tannus et al. 2017), significantly less than the lifetime of a massive star, they should have been ejected by dynamical interactions. We compare the percentage of ejected O stars to other clusters where dynamical ejections have been observed and to the runaway percentage in numerical simulations of star clusters. Dynamically ejected O-type runaways have been observed for an increasing number of young massive clusters. In-depth investigations on these runaways with Gaia have now been done for Westerlund 2, NGC 3603, NGC 7000, the Orion Nebula Cluster and NGC 6611 (Drew et al. 2018, 2019; Maíz Apellániz et al. 2022b; Farias et al. 2020; Stoop et al. 2023). Westerlund 2 is estimated to have eight O star runaways, and considering that there are still 30 O stars in the centre of the cluster, this would translate to an O star runaway percentage of ∼21% (Vargas Álvarez et al. 2013; Drew et al. 2021). On the one hand, undiscovered O + O binaries may cause this percentage to be overestimated. On the other hand, more O star runaways might exist with velocities below the adopted threshold of |Δv_T|≳20 km s⁻¹. Our adopted runaway velocity threshold of > 5 km s⁻¹ makes direct comparisons more difficult. We have discovered four O star runaways with |Δv_T|≳20 km s⁻¹, yielding a runaway percentage of ∼17% above this velocity.

NGC 6618 is similar in O star content to NGC 6611 (∼21 compared to 19 respectively). Adopting the same |Δv_T| threshold of ≳5 km s⁻¹, we obtain 4–5 O star runaways in NGC 6611 (or ∼21–26%). Several O stars in NGC 6611 are more likely to be radial runaways (see Table 2 in Stoop et al. 2023), and Maíz Apellániz et al. (2018) also note that up to 50% of the runaways may not be found if only considering the proper motions. The analysis of Westerlund 2, NGC 6611 and NGC 6618 suggest that O star runaway percentages of 15–20% are typical for 2D velocities above 20 km s⁻¹, but that a considerable fraction of O stars can also have 2D velocities of 3–20 km s⁻¹. In the case of NGC 6611 and NGC 6618, another 15–20% may have 2D velocities in this range. In order to make more robust statements on the fraction of runaways and their velocities, more detailed investigations are needed.

Almost all runaways move with transverse velocities of ∼10–20 km s⁻¹, except for the high-velocity runaway 2MASS J18182392-1721517 with |Δv_T|∼65 km s⁻¹. It is unclear why these 12 runaways have such similar transverse velocities. There could be a common physical origin where the initial star formation, stellar density and stellar dynamics resulted in similar ejection velocities. This is speculative at best and we note that the radial motion of these runaways have not been taken into account yet and that the true velocity distribution could look different.

6.3. Constraints for the dynamical calculations of young massive clusters

The runaway properties can be compared to dynamical simulations of young massive clusters to constrain cluster properties such as the initial radius, mass segregation, binary fraction, separation, and mass ratio (Fujii & Portegies Zwart 2011; Oh & Kroupa 2016). We compare our results to the dynamical simulations of Oh & Kroupa (2016). The ejection fraction of O stars (∼0.3) agrees best with models that assume a relatively small half-mass radius r_h(0) of ∼0.3 pc and a high initial binary fraction⁴. Moreover, most of these models assume an initial mass segregation, with more massive stars located closer to the cluster centre, and mass-ratios (between the secondary and primary star) that tend to unity. The best value r_h(0) agrees well with our results that the O stars were all ejected from within the inner ∼0.1–0.2 pc. The one model with r_h(0) = 0.1 pc has a significantly higher O star ejection fraction of ∼0.52. Models with r_h(0) = 0.8 pc can not produce O star ejection fractions comparable to our results. A pairing of primary and secondary masses that is random causes the massive stars to preferentially be in binaries with lower mass stars, as these dominate the total number of stars. Such a pairing is unfavourable relative to an ordered pairing, which produces mass ratios that are close to unity.

The peak of the velocity distribution of the ejected O stars in the best fitting Oh & Kroupa (2016) models are in excellent agreement with most of the observed transverse velocities of ∼10–20 km s⁻¹. These models predict that most O stars are ejected 0–1 Myr after the start of the simulation, but some may be ejected at late as 2–3 Myr. Interestingly, Oh & Kroupa (2016) find that the binary properties of the ejected O stars differ for these five models. Their multiplicity fraction ranges from ∼0.1–0.2 for model MS3OP to ∼0.4–0.6 for MS3OP_SP. The orbital period, mass ratio, and eccentricity distribution of the runaway population also differ significantly for the five models. MS3OP_SP and MS3OP_SPC have typical mass ratios between 0.8 and 1.0 as a result of the ordered pairing of massive binaries. We have observed two runaway binaries, BD–16 4832 and BD–16 4826, the latter has an orbital period P_orb ∼ 16 days. The mass ratio of these two binaries are ∼0.6–0.8 and 0.5, respectively. Due to the small number statistics, we can not draw definitive conclusions on which specific model setup yields the best fit. At face value, our results best match model MS3UQ_SP. This model has r_h(0) = 0.3 pc, assumes primordial mass segregation for massive stars, a uniform binary mass fraction distribution, and the Sana et al. (2012) orbital period distribution for massive stars.

More investigation is needed regarding the dynamical evolution of young massive clusters to draw conclusive statements on their initial properties. The models of Oh & Kroupa (2016) are tailored for young massive clusters with ∼10 O stars initially, while NGC 6618 formed ∼21 O stars. The binary properties of both the O stars currently in the cluster and runaways could be key in uncovering the initial conditions of massive binaries. The binary properties such as the P_orb, mass ratio, and eccentricity help break degeneracies in dynamical simulations. Ramírez-Tannus et al. (2021) find that massive binaries have significantly wider orbits in the 1–2 Myr after birth compared to the massive binaries typically studied in massive clusters with ages = 2–6 Myr (Sana et al. 2012). This implies that the Sana et al. (2012)P_orb distribution adopted in N-body simulations such as in Oh & Kroupa (2016) may need to be adjusted to an initially wider configuration. This could alter the O star runaway fraction, velocity distribution, and binary runaway properties.

We list here our findings which may help detailed numerical simulations of star clusters to reproduce our findings. There is a period of ∼500 kyr during or shortly after star formation where 30% of the O stars are ejected. The binary fraction of these runaways is not 0, as two of the four spectroscopically observed runaways are found to be binaries. The O star runaways should all be ejected within 0.2–0.3 pc from the centre of the cluster. The 2D velocity of most the runaways is found to be in the range between 10–20 km s⁻¹, and a fast runaway with a 2D velocity > 30 km s⁻¹ is also found. The masses of the runaways are higher than ≳3 M_⊙ considering our cut-off at K_s < 0 mag (Pecaut & Mamajek 2013).

6.4. Age of NGC 6618

NGC 6618 is especially interesting because it is suggested to be one of the youngest clusters known in our Galaxy. We have identified runaways with a kinematic age ranging from ∼100 to 1250 kyr. The O star runaways specifically have a kinematic age between ∼100 and 600 kyr. Simulations of young star clusters show that dynamical ejections commence during and right after star formation (Bate et al. 2002; Fujii & Portegies Zwart 2011; Oh & Kroupa 2016). The kinematic ages of these runaways therefore convey information about when NGC 6618 was at its densest. The runaways that were ejected first may have been lost during or closely following the star formation process. We propose that we can use the kinematic age of the first runaways to estimate the age of NGC 6618. To refrain from using a single measurement to determine the age, we split the runaway sample into two halves based on their kinematic ages. The second sample contains the seven runaways with the largest kinematic ages, which range from ∼400 to 1250 kyr. We perform a Monte Carlo simulation to take the uncertainty on the kinematic age into account. For each runaway, a random t_kin is drawn between its lower and upper bound assuming a uniform distribution. We calculate the mean and standard deviation on these t_kin in each iteration. The age and uncertainty are given by the 50th percentile of the means and the 84th percentile of the standard deviations, respectively. This results in an age of NGC 6618 equal to 0.65 ± 0.25 Myr.

This method of using the kinematic age of the runaways to match to the age of a cluster also works well for several other young clusters. Stoop et al. (2023) show that for NGC 6611 the isochrone age = 1.3 ± 0.2 Myr. Using the 50% of runaways with the largest kinematic ages in their Table 2 and the method described above, yields an age of 1.35 ± 0.31. For the Orion Nebula Cluster (ONC), the age of the oldest population is within 1σ uncertainties 2.51–3.28 Myr (Beccari et al. 2017). The runaways μ Colombae and AE Aurigae were ejected ∼2.5 Myr ago from the ONC, consistent with the age of their natal cluster. Kroupa et al. (2018) also propose that μ Colombae and AE Aurigae were ejected during or shortly after star formation to explain the two younger populations found in the ONC. NGC 7000 also ejected several massive stars 1.5–2.0 Myr ago, consistent with the isochronal age of ∼1.8 Myr (Kuhn et al. 2020; Maíz Apellániz et al. 2022b).

The runaways in NGC 6618 could only have been ejected if the “bully” stars were already present (Fujii & Portegies Zwart 2011). The kinematic age of the first ejected runaway also gives a lower limit on when the first stars formed in the cluster. The first ejected O star is the high velocity runaway 2MASS J18182392-1721517 with t_kin between 0.54–0.58 Myr. This suggests that the first stars in NGC 6618 formed at least longer than 0.54 Myr ago, which is in agreement with the previously determined age = 0.65 ± 0.25 Myr. While there are three B-type runaways with a larger kinematic age, it is unclear if they truly came from NGC 6618. They have lower transverse velocities, making it harder to determine their kinematic ages. Two of these were also ejected further away from the centre of NGC 6618 (0.5–1.0 pc), making their origin less clear. They could have been ejected from a neighbouring star-forming sub-cluster, or they could be interloper field stars. O stars are relatively rare, and almost all originate from young massive clusters or OB associations (Gies 1987; de Wit et al. 2005). We can therefore be more confident that 2MASS J18182392-1721517 is a true runaway coming from NGC 6618.

6.5. Physical properties of NGC 6618

NGC 6618 is reported to contain at least 15 O stars (Povich et al. 2009; Ramírez-Tannus et al. 2017). If we include the O star runaways and several O stars “hidden” in the nearby field, we instead have 21 O stars total born in NGC 6618. We could still miss out on more O stars shrouded by extinction, or hidden in binaries. The O star population of NGC 6618 may therefore be at least 30–40% larger than initially thought. We can make a rough estimate of the cluster mass (M_cl) by extrapolating an initial mass function (IMF) down to lower masses. We adopt here the Kroupa (2001) IMF and assume initially 21 O stars with masses ≳18 M_⊙. Integrating the IMF down to the brown dwarf lower limit of 0.08 M_⊙, we obtain M_cl ∼ 5.1 × 10³ M_⊙. We note that there is a significant uncertainty in this calculation since we are extrapolating the small number of O stars down to the lower mass stars, which make up most of the mass of a cluster.

The runaway O stars are typically not accounted for in the mass function of young massive clusters. It may be that these resemble a Kroupa mass function before accounting for runaway O stars. The fraction of O star runaways is thought to be significantly higher than for B or later-type stars (Gies 1987). If these runaway O stars are accounted for in the mass function, this may result in a significantly more top-heavy distribution than commonly assumed. This has also been noted by Maíz Apellániz et al. (2022b), who find that NGC 7000 may also have a top-heavy IMF if the runaway O stars are accounted for. Schneider et al. (2018) also find a top-heavy IMF for the 30 Doradus region in the Large Magellanic Cloud, where the power-law exponent of the IMF was determined to be ∼1.9 for masses above 15 M_⊙. Since this study focused on the entire surrounding region and not on the R136 star cluster, there were almost no biases regarding either the inclusion or exclusion of runaways. Larger sample studies on a top-heavy Galactic IMF are needed to investigate whether indeed the IMF features a kink to a shallower slope above ∼15 M_⊙.

6.6. Previous star formation

The presence of several massive stars, excess soft X-ray emission and an extended bubble have been used as evidence for a generation of older stars making up the progenitor OB association NGC 6618 PG (Povich et al. 2007, 2009). To further investigate this, we have collected all stars which could be part of NGC 6618 PG in Table 5. These are made up of the stars listed in Povich et al. (2009) and possible OB stars in Reed (2003) located inside the M 17 EB region. We have excluded stars that we now know to have a late spectral type. We list their K_s magnitude, distance of the selected stars, and whether we have identified these stars as a runaway from NGC 6618 instead.

We have shown that four of the seven massive stars that could be part of NGC 6618 PG are instead runaways coming from NGC 6618. One of the dominant ionising sources, the now spectrally classified O9.7 Ia star BD–16 4831, is likely not associated with either NGC 6618 or NGC 6618 PG as this star is located significantly further away at a distance of ∼2.4 kpc. The remaining O star BD–16 4834 has poor astrometry and we can not draw any conclusions to which cluster this star might belong. The presence and extent of NGC 6618 PG is thus significantly less clear than initially thought. NGC 6618 PG would then consist of only B or later-type stars, which will have significant implications for its age. If the O stars do not make up NGC 6618 PG, it becomes less clear what caused M 17 EB.

Instead, we investigate whether the O star runaways ejected in the direction of the M 17 EB could have caused this faint diffuse H II region (Povich et al. 2009). We can assume that the ionising power of all runaways is dominated by the O5 V star in BD–16 4826, as the most massive star easily dominates the ionising budget. This gives log(Q₀) = 49.26 ionising photons per second (Martins et al. 2005). The ionising power of the other runaways combined, which are all O9-B1 II/V, is still an order of magnitude less than the O5 V star. We assume a Strömgren sphere with a radius R_S ∼ 7.5 pc similar to Povich et al. (2009) and a recombination coefficient α_B of 2.57 × 10⁻¹³ cm⁻³ s⁻¹ (Storey & Hummer 1995). We can calculate what the particle density n in the Strömgren sphere must be, which is given by

$$\begin{array}{c}\hfill n_{\mathrm{ISM}}=\sqrt{\frac{3Q_0}{4\pi {\alpha }_A{R}_S^3}}·\end{array}$$(2)

Given these values, we determine n_ISM of the order of 40 cm⁻³. Povich et al. (2009) determine n₀ ∼ 350 cm⁻³ from the dynamics of the CO and molecular gas, which includes two denser surrounding molecular clouds: M 17 North and part of MC G15.9-0.7. Our estimate of n lies in between typical particle densities of the interstellar medium (∼1 cm⁻³) and the n₀ determined by Povich et al. (2009). Our estimated n_ISM is therefore not surprising, given that the H II region lies in between the two aforementioned molecular clouds.

We show the spatial distribution of the OB stars on a larger scale in Fig. 11 in the direction of M 17 EB. We have included the candidate OB stars catalogued in Reed (2003) in the field around M 17 EB as well for context. The only constraint we have put on the Reed (2003) candidate OB stars is that they should be located between 1.5 and 2.0 kpc to exclude obvious fore and background stars.

Fig. 11. Large-scale spatial distribution of the O stars and runaways originating from NGC 6618. We include for context the candidate OB stars in the Reed (2003) catalogue. Their markers are coloured according to their AV with the colour-bar on the left. For the entire field, we show a number-density map coloured according to the Gaia source density with the colour-bar on the right. Each pixel is 1.2′ by 1.2′ in size. We show for context the location of M 17 and M 17 EB with the solid and dashed white circle, respectively.

Still, an over-density of OB stars is present at the location of NGC 6618 PG, even if we account for five runaways. We also show in Fig. 11 for the entire field a 2D histogram of all Gaia sources, coloured according to their source density and a colour-bar on the right. Each pixel is 1.2′ by 1.2′, with brighter colours indicating more sources than darker colours. If we compare Fig. 11 to the extinction map in Fig. 9 in Povich et al. (2009), we see that these Gaia source densities match the extinction map well. The Gaia source density can be used as an analogue to an extinction map, where brighter colours indicate less extinction so that more sources are visible, and vice versa for darker colours.

The dark lanes around M 17 highlight the extreme extinction found in this region (A_V > 5 mag), so that few sources are visible with Gaia in the optical. The M 17 EB shows significantly more sources than the dark lanes surrounding it. If the ionising star in BD–16 4826 has carved out this H II region, it could have decreased the extinction in this line of sight compared to the neutral medium surrounding M 17 EB.

We have estimated the A_V for all runaways, OB star candidates in Reed (2003), and O stars in NGC 6618 (PG) from their 2MASS J − K colour, for which the intrinsic colour is ∼0.24 mag for OB stars assuming the Cardelli et al. (1989) extinction law with R_V = 3.1. The markers of these stars are coloured according to their A_V in Fig. 11, with bluer colours depicting less extinction compared to redder colours with more extinction. The O stars in NGC 6618 can be seen to be reddened (A_V > 5 mag), as is observed (e.g. Ramírez-Tannus et al. 2017).

The sources in NGC 6618 (PG) suffer from less extinction, depicted by their bluer colours (A_V < 5 mag). For example, we estimate that the most blue system BD–16 4832 (B0 V + B1 V) in Fig. 4, located in M 17 EB only has A_V ∼ 2 mag. With less extinction in this line of sight, the early to mid-type B stars are more easily visible and picked up by OB star catalogues such as in Reed (2003). Conversely, fewer OB star candidates are detected by both Gaia and other catalogues in the dark lanes surrounding NGC 6618, which is why we may miss runaways and OB stars in these directions. This can explain why we observe more runaways in the direction of M 17 EB., since we expect runaways to be ejected isotropically. The total number of runaways with stellar masses > 3 M_⊙ could therefore be intrinsically larger and our observed 13 runaways may be a lower limit. We observe six to seven runaways in the direction of M 17 EB. Only an additional four to five runaways were ejected in this quadrant (east) compared to other quadrants, which could also be explained by low number statistics. This does not imply that NGC 6618 PG does not exist and that all the O and B stars make up an “imposter” cluster. We list this as a possibility and note that NGC 6618 PG is likely less massive than initially thought, but may still exist.

6.7. Bowshock

Not all runaways apparently produce a wind bow shock, which is also the case here (Huthoff & Kaper 2002). The runaway star 2MASS J18182392-1721517 can be seen to produce a bowshock in Fig. 7. We investigate whether 2MASS J18182392-1721517 is consistent with being a late-type O star. Baranov et al. (1971) show that momentum balances between the ram pressure exerted by the stellar wind and the ISM. For a given stellar wind velocity v_w, a particle density of the ISM n_ISM, a velocity of the star relative to the ISM v_* and the distance between the star and the apex of the bowshock R₀, the stellar wind mass-loss rate should be

$$\begin{array}{c}\hfill {\dot{M}}_{\ast }=\frac{4\pi m_{\mathrm{H}}{n}_{\mathrm{ISM}}{v}_{\ast }^2{R}_0^2}{v_{\mathrm{wind}}}·\end{array}$$(3)

We assume n_ISM ∼ 1 cm⁻³, typical for what is observed in the ISM. While we have determined n_ISM ∼ 30 cm⁻³ in Sect. 6.6, 2MASS J18182392-1721517 is located in an entirely different region than M 17 EB. Particle densities derived from bowshocks in our Galaxy are typically in the range of 0.1–10 cm⁻³ and our assumed n_ISM seems a reasonable mean (Peri et al. 2012, 2015). We will take v_* to be equal to the transverse velocity |Δv_T| = 65 km s⁻¹ and determine R₀ to be ∼30–60″ which translates to 0.25–0.5 pc at the distance of 2MASS J18182392-1721517 (∼1.75 kpc). Finally, we assume v_wind ∼ 2000 km s⁻¹, which is typical for a late-type O star (Prinja et al. 1990).

We then estimate the mass-loss rate of 2MASS J18182392-1721517 and obtain Ṁ_∗ ∼ 1 × 10⁻⁷ M_⊙ yr⁻¹, or log(Ṁ_∗) = –7.0. This could easily vary by an order of magnitude given the uncertainties and assumptions. Nevertheless, this Ṁ_∗ is typical of mid-to late-type O stars (see e.g. Vink et al. 2000; Björklund et al. 2023). The tentative spectral classification based on its distance and K_s-magnitude is consistent with that needed to produce the bowshock. Optical spectra of this star could conform our findings.

7. Summary and conclusion

We have studied NGC 6618, located in M 17, and it is one of the youngest massive clusters known in our Galaxy with Gaia DR3. We summarise our findings here:

• We have spectroscopically classified six new O stars, of which five are likely associated with NGC 6618, in addition to the 17 known O stars in the cluster. In addition to this we have classified four B-type stars, of which two have been identified as SB2 B-type binaries.

• We have identified 42 members of NGC 6618. M 17 and NGC 6618 are significantly closer than initially thought at a distance of 1.65–1.70 kpc. This has implications for studies investigating the rich pre-main-sequence population that is present.

• We have found 13 runaways coming from NGC 6618, consistent with the dynamical ejection scenario. Out of the 13 runaways, six to seven are likely to be O stars and the remaining are likely to be B stars. The percentage of O stars ejected from NGC 6618 is 27–35%. These runaways are typically not accounted for in studies of the IMF of young clusters, which raises questions on whether the IMF could be more top-heavy than initially thought.

• The accurate Gaia astrometry has allowed us to determine the transverse velocity, kinematic age, and impact parameter with great precision. The runaways have transverse velocities ranging from ∼10 up to 65 km s⁻¹. The kinematic ages range from ∼0.1 to 1.2 Myr, and all runaway O stars have kinematic ages between 0.1 to 0.6 Myr ago. Most of the runaways are ejected from inside 0.2 to 0.3 pc from the centre of NGC 6618.

• We find that the kinematic ages of the runaways can be used to estimate the ages of young clusters. Using the 50% of runaways with the largest kinematic ages (ejected longest ago), we obtain an (kinematic) age = 0.65 ± 0.25 Myr.

• We have discovered a bowshock in the Spitzer MIPS 24 μm image around one of the runaways 2MASS J18182392-1721517, which likely is a late-type O star.

• We find that the progenitor cluster north-east of NGC 6618, NGC 6618 PG is likely to be less massive than initially thought. The main ionising sources in NGC 6618 PG are either runaways coming from NGC 6618 or background stars. Instead, the emission bubble present at the location of NGC 6618 PG could have been caused by the ejected O stars.

M 17 provides a unique opportunity to study a population of massive stars that has just arrived on the main sequence, to constrain the end of star formation properties, and to constrain the early history of dynamical interactions. Increasingly more accurate Gaia astrometry will reveal more details about the young population present, including multiplicity properties of both the runaway and cluster stars.

Acknowledgments

The data and data products presented in this paper are available at the following DOI: https://doi.org/10.5281/zenodo.8120575. M.S. acknowledges support from NOVA. Some of the observations reported in this paper were obtained with the Southern African Large Telescope (SALT) under programme 2022-1-SCI-002 (PI: Paul Groot). Based on observations collected at the European Southern Observatory under ESO programme 0103.D-0099. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

References

Appendix A: Spectroscopy, data reduction, and spectral classification

We have obtained spectra for 10 systems in the O star sample in either with the High Resolution Spectrograph (HRS) mounted on the Southern African Large Telescope (SALT) in Sutherland, the Intermediate-dispersion Spectrograph and Imaging System (ISIS) mounted on the William Herschel Telescope (WHT) at La Palma, or the medium resolution spectrograph X-shooter on the Very Large Telescope (VLT) at Paranal. We describe these observations and classify these stars based on the spectra below. The O stars spectrally classified here are listed first in Table 2. Stars spectrally classified here as early B-type are stars at the bottom of the table.

SALT/HRS

High-resolution spectroscopic observations have been taken of two stars, BD-16 4832 and BD-16 4834, with the HRS, a dual-beam fibre-fed, white-pupil, échelle spectrograph, in low resolution mode. The wavelength range is from 370-890 nm, divided in a Blue and Red arm, at a resolving power of R ∼ 14,000. The spectra were taken as part of a monitoring campaign of three MYSOs in M 17 (programme 2022-1-SCI-002). The data is reduced with the HRS MIDAS pipeline developed by Kniazev et al. (2016, 2017) and includes a flat field reduction, wavelength calibration and object subtraction.

WHT/ISIS

Eight stars have been observed with the ISIS optical spectrograph that includes a blue arm (388-478 nm) and red arm (815-904 nm). The slit width used for the observations is 0.79″ resulting in a resolution of R ∼ 6000 for the blue arm and R ∼ 13,000 for the red arm. An ISIS reduction pipeline was developed in PYTHON, and allows for bias subtraction, flat fielding, overscan and cosmic ray removal, spectrum, and sky extraction and wavelength calibration, based on scripts from Craig et al. (2017) and Crawford et al. (2021).

VLT/X-shooter

B181 has been observed with X-shooter, VLT, a multi-wavelength (300-2480 nm) spectrograph with three separate spectroscopic arms. The resolving power of the arms is R ∼ 5400 in UVB arm (slit width 1.0″), R ∼ 8900 in the VIS arm (slit width 0.9″) and R ∼ 8100 in the NIR arm (slit width 0.6″). The data are taken as part of a multiplicity study in M 17 (programme 0103.D-0099). The pipeline for X-shooter (esorex X-shooter 3.3.5) provided by ESO (Modigliani et al. 2010) is used for flat field correction, bias subtraction and wavelength calibration.

BD-16 4832

The HRS spectrum of BD-164832 reveals a SB2 nature of the source. One of the two stars in the spectrum shows weak He II 4686 Å. Both stars have a strong He I 4471 Å absorption compared to the Mg II 4481 Å, but for the star with the weak He II line, this ratio is slightly larger. Therefore the stars are classified as B0 + B1. The spectrum does not reveal signs of a luminosity class brighter than V (because of only weak O II 4350 Å absorption, close to Hγ), however, close-by absorption features caused by an extended atmosphere are easier obscured in a SB2 spectrum.

BD-16 4834

BD-164834 has an HRS spectrum that shows He I and He II lines and where the relative strength between the He I 4471 Å and He II 4542 Å points to a late O-type star. Since He I 4144 Å is stronger than He II 4200 Å and the latter is slightly stronger than C III 4187 Å, the spectral type is about O9. This is narrowed down to O9.5 by the similar strength of He II 4686 Å and C III 4647/50/51 Å. The luminosity class is I-III by the strength of the N II 4379 Å compared to He I 4387 Å. This is narrowed down to II by the relative strength of Si II 4116 Å to He I 4121 Å and the lack of emission features. Therefore, BD-164834 is of spectral type O9.5 II.

BD-16 4831

The ISIS spectrum of BD-164831 show weak He II lines pointing to a late O or early B spectral type. The classification based on photometry of Povich et al. (2009) is O4 V or O7 III. Since O II 4349 Å is not as strong as is expected for B-type stars, this is a late O-type star, more specifically O9.7 due to the similar strength of He II 4200 Å and C III 4187 Å. Based on the intensity of the Si IV 4089 Å and Si IV 4116 Å lines compared to hydrogen the star classifies as a supergiant. Since the emission lines of N III 4634/42 Å are absent, this is a luminous supergiant; therefore its classification is O9.7 Ia.

BD-15 4928A

This spectrum, taken with ISIS, shows strong He I lines and one He II line at 4686 Å. This points to an early B type classification, different than the O8 V classification from Povich et al. (2009). He I 4713 Å is stronger than He II 4686 Å and there is no Ni II 3995 Å which corresponds to a spectral classification of B0.5. The luminosity class is dwarf since O I 4070/76 Å and Si II 4552/68/75 Å are not particularly strong. The classification of this star is B0.5 V.

BD-15 4928B

BD-154928B is nearby BD-154928A on the sky as reported by Povich et al. (2009) with a spectral type B0.5 V. This star has strong He I lines but no He II, indicating an early B spectral type, but not B0. Since He I 4144 Å is slightly stronger than He I 4121 Å, the star is around spectral type B1. This is B1.5 because of a stronger He I 4713 Å compared to the blended O II + C III 4650 Å which also points to a dwarf luminosity class. The spectral type of BD-154928B is B1.5 V.

B260

B260 has an ISIS spectrum and is previously classified as O6-8 V (Hanson et al. 1997; Hoffmeister et al. 2008; Povich et al. 2009). The spectrum reveals strong He I and He II lines pointing to a late O or B0-type star. The relative stronger He II 4686 Å compared to He I 4713 Å suggests late O-type and the stronger He II 4542 Å compared to Si III 4553 Å, the latter being weaker than He I 4387 Å suggests O9.5. The luminosity class is V based on the strength of C III 4647/50/51 Å, and the strength of the Si II 4116 Å line compared to He I 4121 Å. The classification of this star is O9.5 V.

LS4941

The ISIS spectrum of LS4941 shows He I and He II lines and a strong He II 4686 Å compared to He I 4713 Å, ruling out the classification of Povich et al. (2009) of B1 V and instead suggesting a late O-type star. The stronger He II 4542 Å compared to Si III 4553 Å and stronger C III 4647/50/51 Å compared to He II 4686 Å are associated with spectral type O9.7. Due to the strength of Si II 4116 Å compared to He I 4121 Å, this star is a dwarf star with spectral type O9.7 V.

LS4943

Similarly to LS4941, the ISIS spectrum of LS4943 shows He I and He II lines and a strong He II 4686 Å compared to He I 4713 Å, which is suggesting a late O-type classification, in line with the O9 V photometry classification of Povich et al. (2009). By the same reasoning as for LS4941, this star is classified as O9.7 V.

LS4972

Povich et al. (2009) labelled this star as B2 V, but the ISIS spectrum of LS4972 reveals a SB2 nature with similarly strong He I lines in both components. The lack of He II lines and strong He I 4471 Å compared to the Mg II 4481 Å, suggests B1-B4 for both of the SB2 components. Since He I 4121 Å is similar in strength to He I 4144 Å, the stars are of B1-B2 spectral type. As the stars show a slight difference in strength of O I 4070/76 Å, He I 4144 Å and He I 4387 Å, the classification is B1 + B2. Both stars show a relative lack of metal lines, yielding B1 V + B2 V.

B181

Povich et al. (2009) classified B181 (or CEN61) as O9V based on photometry and Hanson et al. (1997) as O9-B2 based on its K-band spectrum. The X-Shooter spectrum of B181 reveals several He I lines and a few He II lines. The relative strengths of these lines point to a late O or early B-type star. The relative strength of the He II 4686 Å, C III 4647/50/51 Å and He I 4713 Å and the weak He II 4200 Å, of similar strength to C III 4187 Å, point to a spectral type O9.7. The strength of the Si IV 4116 Å compared to its neighbouring He I 4121 Å points to luminosity class III. Therefore, the spectral type of B181 is O9.7 III.

BD-16 4822

The star has been tentatively classified as O6.5 V by Povich et al. (2009) based on photometry, but a lack of He II lines in the ISIS spectrum points towards an early B classification rather than an O-type star. The relative strength between He I 4471 Å and Mg II 4481 Å is indicative for spectral type B1-B3. Because of the presence of C III 4068/70 Å and the strength of He I 4121 Å compared to Mg II 4481 Å this is specified to B2.5. The luminosity class of BD-164822 is I-III due to the relatively narrow shape of the hydrogen lines and the presence of N II 3995 Å. More specifically, the relative strength of Si II 4128-4130 Å points to luminosity class II, such that we settle on B2.5 II.

Appendix B: Members

All Tables

All Figures

Fig. 1. Membership probability p distribution for the 498 stars selected in the direction of NGC 6618 with parallax < 1.0 mas. The field stars have p < 0.96, shown in red, while member stars have p ≥ 0.96, shown in blue. The dashed black line (N = 5) shows the threshold used for the membership p cut-off between field and member stars. The first and second bin from the left contain 250 and 36 stars, respectively, and are cut off for better visualisation of the other bins.

In the text

Fig. 2. Position of all identified O stars or O star candidates in the search field centred on NGC 6618 and NGC 6618 PG, save one. The O star candidate 2MASS J18182392-1721517 is located outside of the field shown. Identified member stars are shown with white plusses; those that are O star (candidates) as magenta circles. O star (candidates) that are not identified as member of NGC 6618 are shown as blue circles. For context, we show the DSS2 B, R and I colour image and the Very Large Telescope Survey Telescope OmegaCAM image (ESO/INAF-VST/OmegaCAM). The inset zooms in on the centre of NGC 6618. The green cross in this inset indicates the determined cluster centre. With dashed circles the location of the emission nebulae and open cluster(s) are shown. The bars denoting physical size correspond to ∼0.2 deg (6 pc) and ∼0.03 (1 pc).

In the text

Fig. 3. Proper motion and parallax for members, all identified O-type stars, and field stars in the search field centred on NGC 6618. The members and O-type stars are coloured and marked similarly as in Fig. 2. O-type stars with ruwe > 1.4 are not shown in both panels. Left: Proper motion distribution in the Galactic coordinate frame. We show for context the transverse velocity difference equal to ∼1 mas yr−1 proper motion difference at the distance of NGC 6618 (1674 pc). The cluster proper motion determined with the mean of the members and from the runaways (see Sect. 4) are shown with the green and red cross respectively. Right: Parallax distribution as a function of the G-magnitude. The best-fit parallax and 1σ uncertainty of NGC 6618 determined from the members are shown with the grey dashed line and bar respectively.

In the text

Fig. 4. Colour – absolute magnitude diagram for members, O-type stars, and field stars in the search field centred on NGC 6618 (PG). The members and O stars are coloured and marked similarly as in Fig. 2. Absolute magnitudes are determined from individual parallaxes. Shown are three PARSEC isochrones with AV = 2.0 mag and age of 0.5, 1.0 and 2.0 Myr with the solid, dashed and dot-dashed tracks, respectively. We give the reddening line (RV = 3.1) for an O9–9.5 V star with the dashed black line and superimpose on this line the position of 2–14 mag of AV, in steps of 2 mag, with the red crosses.

In the text

Fig. 5. Colour–colour diagrams for the O stars (blue), late-type stars (red) and early-type runaways (green). Left: 2MASS colour–colour diagram. We indicate the main-sequence with the solid grey line and highlight the position of an M0, A0, and B2 or earlier main sequence star with a grey label and dot. Reddening lines for an M0 V and a B2 V star are plotted assuming RV = 3.1. Right: Gaia and 2MASS colour–colour diagram.

In the text

Fig. 6. Position and motion of the found runaways relative to NGC 6618, with the DSS2 B, R, and I colour image and the Very Large Telescope Survey Telescope OmegaCAM image (ESO/INAF-VST/OmegaCAM). We show the proper motion of the runaways relative to NGC 6618 with the arrows, which are proportional to their absolute proper motion. Runaways which are consistent with being O stars in Table 2 are coloured in blue, and coloured red otherwise. We show the position of NGC 6618 (PG) and M 17 (EB) using white circles. The angular size equivalent to 10 pc is given at the distance of NGC 6618 (1674 pc).

In the text

	Fig. 7. Spitzer 24 μm Multi-Band Imaging Photometer grey-scale image of the bowshock around the runaway star 2MASS J18182392-1721517. The blue square shows the position of the star. We show the relative proper motion with respect to NGC 6618 with the blue arrow.
In the text

	Fig. 8. Relative transverse velocity (|ΔvT|) of the runaways as a function of their kinematic age (tkin). The identified O stars in Table 2 are shown in blue; red is used for B or later-type stars. The cut-off |ΔvT| and tkin used in the runaway search are indicated with the dashed black lines.
In the text

	Fig. 9. Absolute MKs magnitude of the runaways as a function of their kinematic age (tkin), marked and coloured similarly as in Fig. 8. The cut-off MKs and tkin used in the runaway search are indicated with the dashed black lines.
In the text

	Fig. 10. Projected impact parameter (rimp, 2D) of the runaways as a function of their kinematic age (tkin), marked and coloured similarly as in Fig. 8. The cut-off rimp, 2D and tkin used in the runaway search are indicated with the dashed black lines.
In the text

Fig. 11. Large-scale spatial distribution of the O stars and runaways originating from NGC 6618. We include for context the candidate OB stars in the Reed (2003) catalogue. Their markers are coloured according to their AV with the colour-bar on the left. For the entire field, we show a number-density map coloured according to the Gaia source density with the colour-bar on the right. Each pixel is 1.2′ by 1.2′ in size. We show for context the location of M 17 and M 17 EB with the solid and dashed white circle, respectively.

In the text