© The Authors 2023

1. Introduction

Stars that end their lives in core-collapse supernovae are referred to as massive stars. They have a mass ≳8 M_⊙ at the end of the formation process and begin their lives as O or early B-type stars. Most of the O stars are in young massive clusters and OB associations. It has become increasingly clear that most, if not all, O stars in the field originate from these young massive clusters or OB associations, but have been ejected as runaway stars (de Wit et al. 2005; Gvaramadze et al. 2012).

The two most important mechanisms for these ejections are the dynamical ejection and binary supernova scenarios, or a combination of both (Zwicky 1957; Blaauw 1961; Hoogerwerf et al. 2001). In the first scenario, two massive stars are initially in a binary orbit. After the primary star explodes as a supernova, the secondary star receives a kick velocity comparable to its original orbital velocity. If the binary remains bound, such a system can for example be observed as a high mass X-ray binary (e.g. van der Meij et al. 2021). If the binary is disrupted, the secondary can be observed as a single runaway and the compact object possibly as a fast moving radio pulsar. In several Galactic supernova remnants, runaway stars are identified which may have been ejected after the supernova (Dinçel et al. 2015; Boubert et al. 2017).

In the dynamical ejection scenario, the gravitational interaction between single-binary or binary-binary stars can produce a runaway star (Leonard & Duncan 1988; Fujii & Portegies Zwart 2011). Due to their nature, these interactions are most efficient at high stellar density. Young massive star clusters are ideal environments, where a high stellar density is present before possible expansion sets in. For example, runaway O-type stars have been identified coming from Westerlund 2 and NGC 3603 (Drew et al. 2018, 2019). Both clusters are estimated to be younger than ∼2.5 Myr (Harayama et al. 2008; Pfalzner 2009), making these clusters too young to have produced runaways through the binary supernova scenario. The overall importance of dynamical interactions has remained elusive and for the intra-cluster dynamics, we often rely on simulations (see e.g. Fujii & Portegies Zwart 2011, 2013, 2014; Fujii et al. 2012; Oh et al. 2015; Oh & Kroupa 2016).

Typically, we assume that stars in young massive clusters originate in a burst of star formation. Even accounting for an age spread of approximately several 100 kyr, it has become difficult to reconcile this notion for several clusters. For example, the Orion nebula cluster is consistent with having three discrete episodes of star formation each ∼0.5 to 1.0 Myr apart (Beccari et al. 2017; Jerabkova et al. 2019). Similarly, Westerlund 1 and its stellar content are inconsistent with having a single age, but require an age spread of several million years (Beasor et al. 2021). The R136 cluster in the 30 Doradus region in the Large Magellanic Cloud also shows an age spread among the massive stars of several million years (Schneider et al. 2018).

The study of young massive and open clusters has seen great improvement with the advent of Gaia (Gaia Collaboration 2016). With the precise Gaia astrometry, the member and field stars can now be separated in greater detail than ever before without making strong assumptions on the photometric content of the cluster (see e.g. Cantat-Gaudin et al. 2018). Runaway stars can not only be directly identified based on their relative motion, but can also be traced back to their cluster of origin to accurately constrain their ejection time and kinematic age.

The Eagle Nebula (M16) is a well-known and well-studied H II region spanning over approximately 30′ by 30′ on the southern sky. Near the centre of M16 lies the young massive cluster NGC 6611, about ∼5′ in apparent diameter. Currently, NGC 6611 hosts 15 O stars and about 50 early B-type stars. These massive stars have shaped, through their ionising power, the natal molecular cloud, carving out so-called ‘elephant trunks’ (see e.g. Hester et al. 1996). A spectroscopic study of the O star population has revealed a minimal binary fraction of 0.44 (Sana et al. 2009, 2012). O stars are not only found inside the H II region, Gvaramadze & Bomans (2008) find three nearby O stars which could have been ejected from NGC 6611 in the past.

Various photometric and/or spectroscopic studies have been performed throughout the years (Walker 1961; Hillenbrand et al. 1993; Belikov et al. 1999; Evans et al. 2005; Guarcello et al. 2007a; Martayan et al. 2008; Sana et al. 2009). The mass of NGC 6611 is estimated to be about ∼2 × 10⁴ M_⊙ (Pfalzner 2009). The age of the cluster varies throughout the literature depending on the method used and ranges from 0.5 to 5 Myr, with most estimates between 1 and 3 Myr (Walker 1961; Hillenbrand et al. 1993; Belikov et al. 2000; Oliveira 2008). Indications for an older population of low-mass pre-main sequence stars has been uncovered, with a rough age estimate somewhere between ∼8 to 32 Myr (De Marchi et al. 2013). This older generation of stars is spatially more extended, has a lower X-ray luminosity, and occupies a different regime in the colour-magnitude diagram compared to the ‘younger’ population (Guarcello et al. 2012; Bonito et al. 2013).

With Gaia Early Data Release 3 (EDR3; Gaia Collaboration 2021), we have re-analysed the stellar content of NGC 6611 and M16. We present a kinematic history of the O star content of the open cluster. In Sect. 2, we describe the corrections and filters applied to the Gaia data and the clustering algorithm to obtain the members of NGC 6611. Sections 3–7 present our main results regarding its astrometric, photometric and kinematic properties. We discuss our results in Sect. 8 and provide a summary in Sect. 9.

2. Gaia EDR3 data selection

We have searched for candidate members of NGC 6611 in the Gaia Early Data Release 3 (EDR3; Gaia Collaboration 2016, 2021) catalogue by performing a cone-search centred on (α, δ) = (274.7 deg, –13.8 deg) or (l, b) ≃ (16.96 deg, +0.80 deg) with a radius of 0.25 deg, resulting in 5583 sources. We note that the astrometry and photometry for these sources in the full DR3 is exactly the same, except for the additional radial velocities. To ensure the quality of the Gaia data, we have excluded sources with two-parameter astrometric solutions (astrometric_params_solved = 3) and sources with missing photometry entries in either the G, Bp, or Rp-band, leaving 4886 sources. These filters are necessary as the proper motion, parallax and photometry play a key role in the membership selection and subsequent analysis. This dataset was next subjected to a set of corrections before applying any filters. First, we have corrected the G-band flux and magnitude for sources with 2 and 6-parameter astrometric solutions (Gaia Collaboration 2021; Riello et al. 2021). Second, we have applied the Lindegren et al. (2021a) correction to the parallax to account for the zero-point bias. Third, the Cantat-Gaudin & Brandt (2021) correction is applied, adjusting the proper motion of stars with G < 13 mag. Last, we have applied the El-Badry et al. (2021) correction to the parallax error for stars with G-magnitudes in the range of 7 < G < 21 mag.

Next, a set of filters was applied to the data. We briefly describe these and list the remaining number of sources after each step:

– We excluded sources with a renormalised unit weight error (ruwe) larger than 1.4, likely indicating poor quality astrometric solutions for these sources, possibly as a result of binarity. This leaves 4626 sources.

– We removed sources for which, in 10% of the cases, more than one peak was identified in the windows used by Gaia (ipd_frac_multi_peak). This could imply that the source may be a visually resolved double star. This leaves 4548 sources.

– We excluded sources for which the visibility periods (visibility_periods_used) used were 9 or less, possibly indicating that the astrometric parameters (such as the parallax) could be subject to larger systematic errors (see e.g. Lindegren et al. 2018). This leaves 4522 sources.

– We removed sources which had more than one source identifier during the data processing (duplicated_source) potentially indicating astrometric or photometric issues. This leaves 4520 sources.

– We excluded sources for which the statistic measuring the amplitude of the image parameter determination goodness of fit (ipd_gof_harmonic_amplitude) is larger than 0.15, implying that the source could be a double star. This leaves 4310 sources.

– We excluded sources for which the fractional parallax uncertainty (parallax / parallax_error) is smaller than 5, in order to properly distinguish stars belonging to NGC 6611 from fore or background stars. This leaves 2059 sources.

With these filters, we may create biases by excluding binaries and fainter sources. We have not filtered on the astrometric_excess_noise (ϵ_i) and astrometric_excess_noise_sig (D), which measure the disagreement between the observations and best-fitting standard astrometric model (see e.g. Lindegren et al. 2012). This was not done as this excluded nearly every OB star present in NGC 6611. The Gaia Documentation¹ also notes: “In the early data releases ϵ_i will however include instrument and attitude modelling errors that are statistically significant and could result in large values of ϵ_i and D”. We have also not filtered on the corrected G_Bp and G_Rp flux excess factor C^* (Riello et al. 2021), which is expected to be close to 0 (Riello et al. 2021), as this would exclude a major fraction of the OB stars present in NGC 6611.

2.1. Separating cluster members from field stars

With the 2059 high quality sources, we can separate the cluster members and field stars. Similar to Cantat-Gaudin et al. (2018), we have applied the unsupervised membership algorithm UPMASK to the candidate sources (Krone-Martins & Moitinho 2014). UPMASK was applied to the 5-dimensional space (α, δ, μ_α^*, μ_δ, ϖ), the observables, in the following manner:

1. For each star, the observables are randomly drawn from a multivariate normal distribution with means equal to the values of the observables and covariance matrix based on the uncertainties and correlations between the observables.

2. A grouping algorithm, k-means clustering, is applied to the randomly drawn observables in the 3D astrometric space (μ_α^*, μ_δ, ϖ), with k equal to the number of groups.

3. The groups in (2.) are tested whether they are closely distributed in sky-coordinates (α, δ). This is done by evaluating the total branch length l_obs of the minimum spanning tree (MST) of each group and comparing it to the total branch length l of the MST of randomly drawn uniform distributions. This gives for each group and thus each star a binary membership assignment (true or false).

4. Step (1.) to (3.) are repeated n times, resulting in each star having a membership ‘probability’ p; how many times a star was assigned member relative to the total number of iterations n.

In practice, instead of self-assigning the number of groups k in the k-means clustering algorithm, we can set the average number of stars per group since we know the total amount of candidate stars we start out with. Here, we set the average number of stars per group to 15 as recommended by Krone-Martins & Moitinho (2014). The total branch length l and standard deviation σ_l of randomly uniformly drawn distributions are pre-computed containing 3 to 100 stars in circular distributions with radii equal to the previously specified cone-search radius. The stars in each group with l_obs are assigned member for that iteration if l_obs deviates more than σ_l from the expected l.

We adopted a total of 10 000 iterations for n and a cut-off membership probability p of 0.9, balancing the sensitivity to the number of member stars found and the field stars excluded. After applying the UPMASK algorithm, we obtain the probable cluster members. A last step is to exclude outliers by removing stars that deviate more than three median absolute deviations from the median value for μ_α^*, μ_δ and ϖ, leaving 137 stars brighter than G ≃ 17.5 mag that we identify as members of NGC 6611.

2.2. O stars in NGC 6611

We have gathered all known O stars in and around NGC 6611 in Table 1 and list their identifiers and spectral type from the literature. This gives 19 O stars in a total of 16 systems, where we have included the three distant runaway O stars found in Sect. 6 outside of our cone-search region. Only 7 of the 16 O star systems were found to be members of NGC 6611. This is a combination of the position of O stars being outside the cone-search region (due to their runaway nature as we subsequently show), ruwe > 1.4 due to their binary nature, or proper motion directions not intersecting with the cluster proper motion. We have added the missing 9 apparent non-cluster O star systems to the sample of 137 members. This gives 130 B or later-type and 16 O star systems for a total of 146 members (Table C.1). Since several of these O stars are significantly separated from NGC 6611, and/or have different proper motions, only the O stars which were initially identified as a member (‘T’ in Table 1) are used to calculate the cluster centre, radius, motion and distance in the following section.

3. Members of NGC 6611

In order to validate the results of the UPMASK clustering algorithm, we will first evaluate the 5 astrometry parameters used (α, δ, μ_α^*, μ_δ, ϖ). We show a 0.8 by 0.8 deg Digital Sky Survey 2 B, R and I colour image of NGC 6611 and the surrounding nebula in Fig. 1. The dashed circle indicates the cone-search region, in which the members and O stars are shown. Most of the members are not concentrated towards the centre of the H II region, but several arcminutes north-west of this position. We define the centre of NGC 6611 as the position with the highest source density. The centre is determined to be (α_NGC 6611, δ_NGC 6611) = (274.67 ± 0.01°, –13.78 ± 0.01°), indicated with the purple cross. We highlight the spectroscopically confirmed O stars visible in this image and make a distinction between O stars with poor and good quality astrometry (in this case RUWE > 1.4 and < 1.4 respectively). We will motivate that these originate from NGC 6611 and we include them for visualisation.

Fig. 1. DSS2 B, R and I colour image of NGC 6611 and the surrounding nebula (M16). The member stars are shown as white plus symbols. North is up and east is to the left. We highlight the O stars in this image and show these with poor quality astrometry with the red squares, while the ones with good quality astrometry are shown with the cyan circles. Two runaway O stars are not visible in this image. The centre of NGC 6611, the position with the highest source density, is indicated with the magenta cross. At the top of the figure the angular distance equivalent to 4 pc is given for a distance of 1.7 kpc. We show the cone-search region with the dashed circle, with a radius equivalent to 0.25 deg or 7.4 pc. The open cluster Trumpler 32 (age ∼ 0.3 Gyr; Kharchenko & Schilbach 1995) can be seen in the top-right of the image.

We show the proper motion distribution in the top-left panel of Fig. 2. The members of NGC 6611 (blue) are clearly concentrated in proper motion space, relative to the field stars (grey). As in Fig. 1, we highlight the proper motion of all O stars. To determine the average cluster proper motion, we follow step 1 in Sect. 2.1 by drawing the proper motion randomly from a multivariate normal distribution with means equal to the observed values and a covariance matrix based on the uncertainties and correlations. In this step we ignore correlations between the stars and add in quadrature the systemic uncertainties of 0.0112 and 0.0107 mas yr⁻¹ from Lindegren et al. (2021b) to pmra_error and pmdec_error, respectively. Each iteration gives a mean proper motion, with a negligible difference between the mean and median proper motion. The average cluster proper motion is determined with the 50th percentile over the distribution of all mean proper motions (1000 iterations), while the 1σ uncertainty is determined with the 16th and 84th percentiles. This results in (μ_{α^*, NGC 6611}, μ_{δ, NGC 6611}) = (0.21 ± 0.01 mas yr⁻¹, −1.59 ± 0.01 mas yr⁻¹). This derived cluster proper motion is consistent with other estimates based on Gaia (see e.g. Flynn et al. 2022; Maíz Apellániz et al. 2022).

Fig. 2. Astrometric parameters of NGC 6611. Top-left: proper motion distribution of the field stars (grey) and members of NGC 6611 (blue). We show the O stars present in NGC 6611 with poor and good quality astrometry with the red squares and cyan circles, respectively. Top-right: parallax distribution of the field stars and members, coloured similarly as in the top-left panel. We also show the fractional parallax uncertainty (ϖ/σϖ) of the members as a function of their parallax. Two O star systems have RUWE = 2.544 and 8.382, and have parallaxes outside this panel, with only part of the error-bar visible. The 1σ error on the parallax of NGC 6611 is shown with the orange bar. Bottom-left: radial velocity of the members in the heliocentric frame collected from the literature. Bottom-right: radial velocity in the heliocentric frame of the B0 to B4 stars after removal of the outliers in the bottom-left panel. Outliers have been excluded as described in Sect. 3. The mean and standard deviation of the radial velocity of these stars are shown with the dashed black line and shaded region, respectively.

The next astrometry parameter is the parallax. We show the parallax distribution of the field stars and members (including the O stars) in the top-right panel of Fig. 2, adopting equal colouring and marking as before. Of the 146 members, 72 sources have ϖ/σ_ϖ ≤ 20 and 74 sources have ϖ/σ_ϖ > 20, providing (extremely) accurate parallaxes for nearly half of the members. Assuming that the members are at the same distance, we expect that the parallaxes of these stars are normally distributed around the true parallax of NGC 6611. Note that the parallaxes of several of the O stars should not be trusted due to ruwe > 1.4. We find no evidence for a skewed parallax distribution in either the full sample (p = 0.45), the ϖ/σ_ϖ ≤ 20 sub-sample (p = 0.16), or the ϖ/σ_ϖ > 20 sub-sample (p = 0.09). The parallax distribution also shows that the majority of the field stars are located closer than ϖ > 0.65 mas and steeply drops off for stars beyond ϖ < 0.6 mas. Background stars most likely go undetected because of the severe extinction by the natal molecular cloud. We find no indication of contamination of fore and background stars.

Following the advice of Bailer-Jones et al. (2021), we have set up a joint likelihood to estimate the distance directly from the parallax. The joint likelihood is similar to that used in Cantat-Gaudin et al. (2018), with the only difference being the adopted parallax zero-point offset. We have already applied the Lindegren et al. (2021a) parallax zero-point offset to the parallax in Sect. 2; it amounts to a mean offset of −0.039 mas (standard deviation of 0.008 mas) for members with 5-parameter, and −0.051 mas (standard deviation of 0.009 mas) for members with 6-parameter astrometric solutions. The likelihood is given by

$$\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)

Here, P(d | ϖ,σ_ϖ) is the unnormalised probability distribution for the distance d in kpc, given the individual parallaxes ϖ_i in mas, and adopted Gaussian parallax uncertainties σ_ϖi in mas. Also noted in Cantat-Gaudin et al. (2018), we neglect correlations between measurements of all stars, and assume that all stars are at the same distance. We determine the distance with the mode of the unnormalised probability distribution. The positive and negative 1σ error are determined from the 16th and 84th percentile, respectively. We determine the distance and parallax to be 1706 ± 7 pc and 0.587 ± 0.003 mas (symmetric distributions), respectively. The determined values for the centre, proper motion, parallax and distance of NGC 6611 agree with other determinations in the literature (see e.g. Cantat-Gaudin et al. 2018; Kuhn et al. 2019; Maíz Apellániz et al. 2022). Our distance is consistent with the distance of ${1681}_{-10}^{+9}$ pc determined in Cantat-Gaudin et al. (2018) considering the uncertainties in distance and parallax zero-point offset in Gaia DR2 and EDR3, and the distance determination of 1697${}_{-30}^{+31}$ pc in Maíz Apellániz et al. (2022) who include estimates for external uncertainties.

We have collected all known spectral types of stars in the field of NGC 6611 determined from optical spectra in the literature (Hillenbrand et al. 1993; Evans et al. 2005; Wolff et al. 2007; Martayan et al. 2008). In most cases, this also allowed the (systemic) radial velocities to be collected. This gives a total of 166 spectral types and 146 radial velocities. Out of the 146 stars with known radial velocities, 34 were excluded based on their poor astrometry. The remaining 66 stars were removed because of their deviating parallax and/or proper motion. We have cross-matched the members of NGC 6611 with these spectral types and radial velocities. Focusing on the radial velocities of the members, we show them in the bottom-left panel of Fig. 2. This results in 46 members with a known radial velocity.

In order to estimate the radial velocity of NGC 6611, we do not need to analyse the full sample. The radial velocities in the bottom-left panel of Fig. 2 display significant spread even accounting for uncertainties (partially because of binarity). To determine the radial velocity of NGC 6611 we turn to the early type B stars, the B0 to B4 stars. We do not use the O stars, which we treat in Sect. 7. We exclude known binary stars (such as W472; Walker 1961; Evans et al. 2005). We also remove stars with radial velocity uncertainties larger than 15 km s⁻¹. We use a similar approach as Evans et al. (2005) and determine for the B0-4 stars the mean of the radial velocity and exclude stars which deviate more than 2σ from the mean. While this method can be repeated multiple times, only 1 iteration was needed to converge. We show the remaining 17 B0-4 stars in the bottom-right panel of Fig. 2. The mean and standard deviation of the radial velocity is 4.7 ± 3.5 km s⁻¹ shown with the dashed black line and shaded region respectively. This is consistent with the mean radial velocity of NGC 6611 in Evans et al. (2005) of 10 ± 8 km s⁻¹ within the uncertainties.

4. Two stellar populations in NGC 6611

4.1. Gaia photometry

We have constructed the colour – absolute magnitude diagram (CAMD) of the members of NGC 6611 with the Gaia photometry using the determined cluster distance for all members. We display the Gaia CAMD in the left panel of Fig. 3, where the different spectral types (O, B or A) are indicated. We show the reddening line for the O7 V((f)) star W222 with R_V = 3.5 and A_V = 6.49 mag and for visualisation we show A_V of 2.0, 4.0, 6.0 and 8.0 mag with the black crosses (Maíz Apellániz & Barbá 2018). This R_V is similar to the average of ∼3.56 in NGC 6611 (Kumar et al. 2004). The O stars are distributed parallel to this line, likely due to variable extinction. As mentioned, we have included all O stars, including the three distant runaway O stars for completeness.

Fig. 3. Colour – absolute magnitude diagram (CAMD) for the average cluster distance for NGC 6611. Errors are smaller than the respective symbol if not shown. Left: Gaia CAMD of the members with known spectral types. The O stars are shown with the open blue squares, B stars with the open orange circles, A stars with the open red stars and unclassified members in black. The two runaway O stars discussed in Sect. 6 are included. We show a reddening line for an O7 V star with RV = 3.5, the black crosses indicate AV = 2.0, 4.0, 6.0 and 8.0 mag from left to right. Right: extinction corrected MV0 for the average cluster distance (RV = 3.56; Kumar et al. 2004) against the intrinsic colour (B − V)0 from the Guarcello et al. (2007a,b) BVI photometry. The stars are marked and coloured similarly as in the left panel. Stars with an uncertain or no spectral type are not shown. We show three PARSEC + COLIBRI isochrones with ages of 1.0, 5.0 and 10.0 Myr with the solid, dashed and dash-dotted grey lines, respectively.

The B stars are almost all fainter than M_G = –1 mag. The three B stars that are brighter constitute a B0.5 V + B1: system, a B0.5 V star (Evans et al. 2005) and a B2.5 I star (Hillenbrand et al. 1993). The B2.5 I star BD–13° 4912, located at (G_Bp − G_Rp) ∼ 2.0 mag and M_G ∼ –2 mag, stands out specifically as it is the only star classified as a supergiant among all the collected spectral types. The Gaia astrometry of BD–13° 4912 is in excellent agreement with that of NGC 6611. The B stars fainter than M_G = –1 mag show a clear over-density around (G_Bp − G_Rp) ∼ 1.0 mag and M_G ∼ 2 mag. The membership of NGC 6611 includes 4 A stars, located around (G_Bp − G_Rp) ∼ 1.0 mag and M_G ∼ 2.5 mag. Three of these are early type A stars, while the fourth has an uncertain spectral type (Ae).

The spectral types of the members are incomplete and completely missing below M_G ∼ 1 and 3 mag, respectively. Among the members with known spectral type, we identify one outlier: W611. This star is classified as a K0 V star in Evans et al. (2005) and is located around (G_Bp − G_Rp) ∼ 2.0 mag and M_G ∼ 0 mag among the O stars. The parallax of W611 is only marginally consistent at 3σ but would still place it too far away to be consistent with a K0 V star. The 3σ lower limit on the parallax places it at ∼1.4 kpc, which would still make it 5 to 6 magnitudes brighter than a typical K0 V star (Pecaut & Mamajek 2013). de Winter et al. (1997) classified W611 as a G8 III star, with the luminosity classification based on the apparent brightness. We deem W611 to be a foreground star and have therefore excluded it from the members. Other than this one outlier, we find no evidence for significant contamination with field stars. Another possibility is that W611 is a pre-main-sequence star. This would imply that this star has formed more recently than all other stars in NGC 6611.

Towards the fainter end (M_G ≲ 2 mag) of the CAMD, two groupings are visible. A ‘blue group’ is visible starting at (G_Bp − G_Rp) ∼ 1.0 mag and M_G ∼ 3.5 mag and a ‘red group’ is visible starting at (G_Bp − G_Rp) ∼ 1.5 mag and M_G ∼ 2.5 mag. They are clearly separated by a gap. As we know the distance to these stars, there are two key quantities left in explaining their positions in the CAMD; extinction and age. The clear gap between these two tracks argues against variable extinction, although it can explain the spread of the stars within each of the tracks.

4.2. BV photometry

Before drawing further conclusions from the CAMD, we take a look at the (variable) extinction in NGC 6611. To do this, we have manually cross-matched our members with the Guarcello et al. (2007b,a) Johnson BVI photometry obtained with the Wide Field Instrument on the 2.2 m telescope of the European Southern Observatory. The cross-matching provided satisfactory counterparts for all but two members: W186 and Gaia EDR3 4146405584125976576. These two stars have not been included in further analysis of the BVI photometry. This does not affect our results. With the BV photometry and known spectral types, we can determine the colour excess E(B – V) and extinction in the V-band A_V with

$$\begin{array}{cc}\hfill A_{\mathrm{V}}=R_{\mathrm{V}}E(\mathrm{B}-\mathrm{V})& =R_{\mathrm{V}}[(\mathrm{B}-\mathrm{V})-{(\mathrm{B}-\mathrm{V})}_0],\hfill \end{array}$$(2)

where R_V is the slope of the extinction curve and (B – V)₀ is the intrinsic colour of the star. We adopt (B – V)₀ for luminosity class V from Pecaut & Mamajek (2013) and use these intrinsic colours for all luminosity classes. This is an excellent assumption for O stars (Martins & Plez 2006), and only results in deviations of ∼10⁻² mag for the B and A stars. We have not estimated the extinction for the Ae star W112 and the B star W433, as the spectral types are uncertain. We assume the intrinsic colour to be that of the earliest-type star (assumed to be the brightest) in the case of a binary or triple system. These assumptions contribute to the overall uncertainties, however, the adopted spectral type and R_V are likely to be more significant sources of uncertainty. R_V is typically assumed to be 3.1 (Cardelli et al. 1989), but can be significantly higher for H II regions (for NGC 6611 see e.g. Kumar et al. 2004; Maíz Apellániz & Barbá 2018). We will adopt R_V = 3.56 from Kumar et al. (2004), which provides reasonable isochrone fits in Sect. 4.3, while an R_V of 3.1 does not.

For the stars with known spectral sub-types, we show the CAMD with the average cluster distance and the extinction corrected M_V0 against the intrinsic colour (B – V)₀ in the right panel of Fig. 3, using similar colours and markers as in the left panel. The spread in colour which is visible in the left panel has now almost completely disappeared. Almost all of the O and B stars now clearly trace out a main-sequence, which is expected as the (B – V)₀ is nearly equal for the O and B stars, with only M_V0 varying. The intrinsically brightest B2.5 I star BD–13° 4912 stands out, along with the three A stars which may not have reached the main-sequence.

For context, we show the stellar isochrones 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) for ages of 1.0, 5.0 and 10.0 Myr. At the faint-end, the isochrones show the location of stars still on the pre-main-sequence. The transition from the main-sequence to the pre-main-sequence is visible as the isochrone experiences a ‘knee’ and becomes intrinsically redder. For older isochrones, this transition moves towards fainter and redder magnitudes, as the pre-main-sequence lifetime increases for decreasing mass. This gives us a way to determine the age of a young massive cluster (Guo et al. 2021). No isochrone is able to perfectly fit all stars at first glance. Fitting the A stars, the main-sequence and the intrinsically brightest B2.5 I star is already creating problems.

So far, we have only used the members with known spectral types. Including the members without spectral types, which are mostly fainter than M_G ∼ 1 mag could give additional information to better constrain the (variable) extinction and age of NGC 6611.

4.3. Isochrone fitting

To estimate the age of NGC 6611, we have set up a log-likelihood function similar to Jørgensen & Lindegren (2005) estimating the probability of obtaining an observed set of stars given an isochrone. The log-likelihood is given by

$$\begin{array}{c}\hfill \mathcal{L}\propto {\displaystyle \sum _{i=1}\log \left(\sum _{j=1}\frac{1}{2\pi {\sigma }_{x_i}{\sigma }_{y_i}}\exp [-\frac{{(x_i-x_j)}^2}{2{\sigma }_{x_i}}-\frac{{(y_i-y_j)}^2}{2{\sigma }_{y_i}}]\right),}\end{array}$$(3)

where in our case x is equal to (B − V), and y is equal to M_V. We have set up a grid using the PARSEC + COLIBRI stellar isochrones varying the extinction A_V between 0.0 and 5.0 mag in steps of 0.1 mag (with R_V = 3.56; Kumar et al. 2004) and the age between 0.1 and 15 Myr in steps of 0.1 Myr. For each of the isochrones in the grid, we can calculate the log-likelihood with Eq. (3) in order to find the best-fitting isochrone in the entire grid. The best-fit isochrone will have the largest value for the log-likelihood as this curve deviates the least from the fit members. We are dealing with a large number of stars with ‘small’ uncertainties and a large number of points on the isochrones. The log-likelihood value will therefore quickly have large negative values, even in log-space.

For each of the isochrones, we choose an A_V. Stars with known spectral type and therefore known A_V need to be corrected to the A_V adopted in each isochrone. For example, when fitting an isochrone with A_V = 2.6 mag, we would place the stars with known spectral type in the CAMD as if they also had A_V = 2.6 mag by linear inter- or extrapolation.

It is clear from the right panel in Fig. 3 that we can not accurately constrain the age based on the early type main-sequence stars alone. The members without spectral type display, as mentioned, two clear groups (see the left panel in Fig. 3). We have divided the Gaia CAMD in three different groups, the early type stars, a blue group and a red group. The location and distribution of these three groups is shown in Fig. 4, where the red and blue groups have their respective colours and are not corrected for extinction. If the spectral type is known, we show the star as an open symbol. To constrain the average extinction and age, we will fit two sets of data comprised of the three groups in Fig. 4. The first set of data is the combination of the early type stars and the blue group (with and without spectral type). As stated above, the early type stars are assigned an A_V similar to the isochrone in order to make the main-sequence consistent with the blue group. We have excluded the O stars for this set of data, which we will show is a necessary assumption. The second set of data is the combination of the early type stars and the red group (with and without spectral type). Again, the early type stars are assigned an A_V similar to the isochrone to make the main-sequence consistent with the red group.

Fig. 4. Gaia colour – absolute magnitude diagram. We have separated the early type stars, the blue and red group shown with the black circles, blue triangles and red squares, respectively. Stars with a known spectral type are shown with an open marker.

The results for the best-fit age and mean A_V applied to the isochrones are shown in Fig. 5. For the first set of data, the early type stars and the blue group, the best fit isochrone has an average A_V = 2.0 mag and age = 7.5 Myr. This best fit isochrone is shown in the top-left panel, with the stars with known spectral type shifted to the same extinction (stars without spectral type are not shifted). The log-likelihoods of the isochrones are shown in the bottom-left panel. The cyan cross indicates the location of the best-fit isochrone. Similarly, the early type stars and the red group are shown at the right hand side of Fig. 5. The best-fit isochrone to this set of data has A_V = 3.5 mag and age = 1.3 Myr. In this method, we are primarily fitting the mean A_V and age of the pre-main-sequence stars. Since the early type stars have mostly known spectral type, these stars will always align with the isochrone. However, these early type stars are still valuable to the fit as they dictate where the main-sequence should be with respect to the pre-main-sequence. Fitting only the blue or red group allows too much freedom for the isochrones since we would only be spanning 4 to 5 mag in M_G (i.e. the early type stars give constraints to the isochrone fitting).

Fig. 5. Isochrone fitting results of NGC 6611. Top-left: best-fit isochrone with AV = 2.0 ± 0.1 mag and age of 7.5 ± 0.4 Myr shown with the solid cyan line to the early type stars and the blue group (with and without spectral type). An open symbol represents a known spectral type, while a closed symbol does not. The uncertainty on the age for AV = 2.0 mag is shown with the cyan shaded region. The stars with known spectral type are corrected for extinction according to the best-fit AV with RV = 3.56 (Kumar et al. 2004). Bottom-left: log-likelihood for each isochrone in the searched AV-age grid for the stars indicated in the top-left figure. Each point in the grid is coloured according to its (unnormalised) likelihood value. The log-likelihood gives the quality of the fit, with smaller negative values (in lighter grey) representing better fits. The cyan cross and circle indicate the best-fit AV and age with 1σ uncertainty range. Top-right: similar to the top-left figure, but now for the early type stars and the red group (with and without spectral type). The solid cyan line gives the best-fit isochrone with AV = 3.5 ± 0.1 mag and age of 1.3 ± 0.2 Myr. The uncertainty on the age for AV = 3.5 mag is shown with the cyan shaded region. Bottom-right: similar to the bottom-left figure, but now for the top-right figure.

The colour-mapped isochrone grids show more than one local maximum. Specifically, the second set of data including the red group shows seemingly ‘good fits’ with A_V between 4 and 5 mag and an age of ∼3 Myr in the bottom-right panel. These isochrones do not accurately fit the transition from pre-main-sequence to main-sequence (the ‘knee’ in the isochrones) but place this around (B – V) ∼ 1.25 mag and M_V ∼ 5 mag. The colour-map is in log-likelihood space, making fits with a log-likelihood, for example, of 100 less than the best-fit seem acceptable even though the probability difference is a factor 10¹⁰⁰ (due to the small uncertainties on the accurate Gaia photometry).

While the average extinction and age of the two sets of data are clearly different, we need to consider if this is significant. If the log-likelihood value of a fit decreases by 1, the p-value decreases by a factor 10. We would determine unrealistically small uncertainties if we were to use the log-likelihood directly. Instead, we determine the uncertainty on the average extinction and age through bootstrapping. Assuming that 10% of the members are outliers or field stars (since we initially chose a membership probability of 0.9), we randomly select new sets of data containing 90% of the original sets of data, without replacement. We determine the best-fit extinction and age and re-do the bootstrapping and fitting 1000 times. With a 1000 best-fit extinctions and ages (for each set of data), we determine the uncertainty on the extinction and age with the 16th and 84th percentiles. This gives for the set of data with the blue group A_V = 2.0 ± 0.1 mag and an age of 7.5 ± 0.4 Myr. For the set of data with the red group we have A_V = 3.5 ± 0.1 mag and an age of 1.3 ± 0.2 Myr. We show the best-fit values with the solid cyan line and cross in the top and bottom panels of Fig. 5, respectively. The uncertainty on the age (assuming the best-fit extinction) is shown with the shaded cyan region and circle in the top and bottom panels, respectively.

It becomes obvious now that the blue and red group have a significantly different mean A_V and age, so we should consider them as two distinct populations of stars. However, this distinction is not as easily made for the early type stars. If we assume an age of 7 to 8 Myr for the older population of stars, the O stars earlier than O7.5-8.5 should already have gone supernova or are considerably evolved. The presence of an O3.5 V((f+)), an O4 V((f+)) and an O6.5 V((f)) star thus argues that these early type stars are part of the younger population of stars (Sana et al. 2009). If these early type O stars are part of the younger population, we could expect many of the late type O stars and B stars to be so as well on the basis of the initial mass function. All O stars have been excluded in the fitting of the older population. One major assumption in the isochrone fitting is that the (B – V) spread within the blue and red group individually is caused by extinction, which is also observed for the O and B stars. The determined best-fit extinction including uncertainty should only be indicative of the mean extinction value in that group. This is not necessarily the case for the determined ages, as we expect the blue group to be formed in a single star formation episode and the red group to be formed in a different star formation episode, each lasting approximately several 100 kyr. If prolonged star formation took place over several Myr, the observed main-sequence(s) and the two pre-main-sequence tracks, clearly separated by a gap, would be difficult to fit within this scenario.

5. Astrometry of the two populations

Up to this point, we analysed the astrometric properties of all members of NGC 6611 together. There may be differences in the spatial distribution, proper motion, or even parallax between the young and old population (red and blue group) defined in Fig. 4.

The first thing to investigate is the spatial distribution of the two populations, which we show in the top-left panel of Fig. 6, coloured and marked similarly as in Fig. 4. We can see that the young (red) population is clearly concentrated towards the centre of NGC 6611, while the old (blue) population is spatially extended. The majority of the old population is scattered south of the centre of NGC 6611; however, it is unclear whether this is physical or an observational bias (e.g. stars with relatively low extinction are found south of the centre of NGC 6611).

Fig. 6. Distribution of the young and old population (the red and blue group, respectively), coloured and marked similarly as in Fig. 4. Top-left: B, V and R composite mosaic image of NGC 6611 (ESO, press release 0926) showing the spatial distribution. Top-right: proper motion distribution. Bottom-left: parallax distribution. Bottom-right: 2D half-light radius for each population as a function of time. The dashed line indicates the mean radius and the coloured region indicates 1σ uncertainty. The dash dotted lines give the time of minimum radius for each population. The young population shows a minimum –0.40 Myr ago, while the old population shows a minimum –0.41 Myr ago.

We also show the proper motion and parallax distribution in Fig. 6. Both the proper motion and parallax show no clear differences. A two-sided t-test and a Kolmogorov–Smirnov test between the blue and red group for the two proper motion components and the parallax gives p-values consistent within 2σ for all cases. We note that biases could possibly come into play in the membership selection, since the younger population has ∼3 times more stars than the older population (∼5 times when including the early type stars). We could therefore select only stars belonging to the older population with a proper motion and parallax consistent with the young population and exclude outliers. The older population can not be explained with the contamination of field stars in our members. While the spatial distribution of this population might be considered consistent with a field star population, both the proper motion and the parallax are in excellent agreement with that of the younger population. As a consistency check, we have selected stars associated with the older population from the 2878 candidate sources in Sect. 2.1 by their position in the CAMD which include many field stars. Still, a clear over-density of stars is present at the determined proper motion and parallax of NGC 6611. The older population is also redder (or brighter) than the field stars.

The older population may be spatially extended now, but could have been more concentrated in the past, similar to the younger population. To investigate this, we trace back in time the stars belonging to the young and old population separately. We perform 10 000 Monte-Carlo simulations randomly drawing the proper motions for each star from a multivariate normal distribution with means equal to the observed values, with the uncertainties given by the uncertainties and correlations. We trace back the young and old population –2 Myr back and forward in time by propagating their proper motions. In increments of 10 kyr we determine the 2D half-light radius for each population, given by the mean as the 50th percentile and 1σ uncertainty as the 16th and 84th percentile (without correcting for extinction). We show the result in the bottom-right panel of Fig. 6. The young population has a minimum 2D half-light radius of 1.7 pc –0.40 Myr ago, while the old population has a minimum radius of 3.0 pc –0.41 Myr ago. Both populations were smaller in the past and are currently already expanding. The position of the stars and thus the radius are likely extrapolated too far back and forward in time before and after –1 and 1 Myr, respectively.

6. Runaway stars

We have now discussed the astrometric parameters and the isochronal age(s) of NGC 6611. Dynamical interactions in young massive clusters are expected to produce runaway stars (Gvaramadze et al. 2012). The young population has an age of 1.3 ± 0.2 Myr, so core-collapse supernovae should not have happened yet. Runaways should then be produced through dynamical interactions. With the accurate Gaia EDR3 proper motions and parallaxes, we should be able to find runaways that trace back to NGC 6611.

We search for runaway stars coming from NGC 6611 by performing a cone-search in Gaia EDR3 centred on the cluster with (α_NGC 6611, δ_NGC 6611) = (274.67 deg, –13.78 deg) with a radius of 15 deg. These stars are subject to the same corrections and filters described in Sect. 2. Moreover, we require these stars to have ϖ/σ_ϖ > 10 and ϖ to be consistent within 3σ from the parallax of NGC 6611. This assumption may exclude runaways moving in the radial direction, for which the parallax could now significantly deviate. Unfortunately in Gaia EDR3 few radial velocities are available, so there is no viable method for including these. The determined radial velocity for NGC 6611 is 4.7 ± 3.5 km s⁻¹ and its radial motion has a minimal effect on the analysis. Next, we have only included stars brighter than M_G < 1 mag, where we have directly used their (accurate) parallax. We have tried to adjust this cut-off magnitude and found that for fainter potential runaways it becomes increasingly more difficult to distinguish between actual runaways and field stars aligning by chance. As we will show, the adopted cut-off magnitude allows us to separate the early type (O, B or A) stars expected from NGC 6611 from bright late type stars (e.g. red giants) in the field at this magnitude.

We identify runaways by tracing back in time their position and checking whether this position aligns with the position of NGC 6611 at the same time. Stars consistent to be runaways from NGC 6611 should satisfy at some point in time

$$\begin{array}{c}\hfill r_{\mathrm{sep}}(t)=\sqrt{{\alpha }_{\mathrm{sep}}^2(t)+{\delta }_{\mathrm{sep}}^2(t)}<r_{\mathrm{NGC}6611},\end{array}$$

where α_sep, δ_sep and r_sep are the separation between the centre of NGC 6611 and the runaway in right ascension, declination and on-sky projection respectively. r_NGC 6611 is the cluster radius and will be chosen sufficiently small to ensure that the small-angle approximation holds. α_sep and δ_sep are given by

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{sep}}(t)& =(\alpha +\frac{t·{\mu }_{{\alpha }^{\ast },}}{3.6\times {10}^6·\cos (\delta )})\hfill \\ \hfill & -({\alpha }_{\mathrm{NGC}6611}+\frac{t·{\mu }_{{\alpha }^{\ast },\mathrm{NGC}6611}}{3.6\times {10}^6·\cos ({\delta }_{\mathrm{NGC}6611})})\mathrm{deg},\hfill \\ \hfill {\delta }_{\mathrm{sep}}(t)& =(\delta +\frac{t·{\mu }_{\delta ,}}{3.6\times {10}^6})\hfill \\ \hfill & -({\delta }_{\mathrm{NGC}6611}+\frac{t·{\mu }_{\delta ,\mathrm{NGC}6611}}{3.6\times {10}^6})\mathrm{deg}.\hfill \end{array}$$

Here, t is the time in years and should be negative for runaways coming from NGC 6611 and μ_α^* ≡ μ_αcos(δ). We will take a conservative approach and use a lower limit of –3 Myr for t, considering the age of 1.3 ± 0.2 Myr. We adopt a cluster radius of r_NGC 6611 = 0.1 deg. A last filter is applied by only including runaways moving outwards with relatively high transverse velocities with Δv_T > 3 km s⁻¹ and Δv_out > 0 km s⁻¹ with respect to NGC 6611, where Δv_T and Δv_out are, respectively, the transverse velocity and the outwards velocity component relative to the centre of NGC 6611.

Runaways are typically defined to have a relative velocity > 30 km s⁻¹. However, we are limited to a 2D approach as the radial velocities are missing. In the case of the binary supernova ejection scenario, walkaways with a relative velocity < 30 km s⁻¹ are predicted (Renzo et al. 2019). But these definitions are not clear for the dynamical ejection scenario, as this is dependent on the cluster escape velocity. For example, the young massive cluster R136 in the large Magellanic Cloud with a mass of ∼5 × 10⁴ M_⊙ and a radius of ∼1 pc has an escape velocity of ∼20 km s⁻¹ (Schneider et al. 2018). A conservative estimate of the current 3D escape velocity of NGC 6611 is ∼4 km s⁻¹ at a distance of 2 pc from the centre (see Sect. 8.1). If runaways are isotropically ejected, we expect the projected Δv_T to be on average ≳3 km s⁻¹, which we adopt as our lower limit. Kuhn et al. (2019) select members from X-ray, optical and infrared studies and show that stars in NGC 6611 have a median Δv_out of 0.90 ± 0.23 km s⁻¹ and NGC 6611 is therefore on average expanding. Our lower limit on Δv_T is therefore sufficiently large to exclude members of NGC 6611. Note that dynamical interactions may contribute to the expansion of a cluster like NGC 6611. While Δv_T may not be equal to Δv_out depending on the direction of motion, our use of Δv_T instead of Δv_out does not affect our results.

We find a total of 16 candidate runaways coming from NGC 6611 in the last 3 Myr. The first thing to consider is whether these are ‘true’ runaways or field stars aligning with the position of NGC 6611 by chance. As mentioned, we assume the runaways to be early type stars (M_G < 1 mag). To investigate this, we can look at their 2MASS JHK_s photometry. We show the (J – H) – (H – K_s) colour–colour diagram in the left panel of Fig. 7. The position of the zero-age main sequence is shown with the grey track. The early type stars should be located in the bottom-left below (J – H) and (H – K_s) ≲ 0 mag. Extinction will cause stars to move towards the top-right. For example, we show the reddening line (with R_V = 3.1) of a B0 V (bottom) and M0 V star (top) with a dashed black line. Since stars with a spectral type of ∼ A0 V or earlier are more or less expected to have overlapping reddening lines, the early type runaway stars should lie roughly on the B0 V track. This is the case for 8 of the candidate runaways. The other 8 clearly lie above the M0 V extinction track making them inconsistent with the early type stars, but consistent with late type red giant stars. The late type stars also have seemingly random kinematic properties, coming from random positions within the cluster (i.e. not from the centre of NGC 6611). The late type stars also have randomly distributed kinematic ages. For this reason, we consider these 8 stars as field stars aligning by chance and have excluded them.

Fig. 7. Properties of the runaway candidates. Left: 2MASS near-infrared colour-colour diagram for the runaway candidates. The early type runaways originally from NGC 6611 in the last 3 Myr are shown in blue. Late type stars are highlighted in red and are deemed unlikely to come from NGC 6611. We show the location of the zero-age main-sequence in grey (Pecaut & Mamajek 2013). The reddening line for a B0 V (top) and M0 V star (bottom) are shown with the dashed black lines for RV = 3.1 (Cardelli et al. 1989). Right: B, V and R composite mosaic image of NGC 6611 (ESO, press release 0926). We highlight the adopted search cluster search radius, with a radius of 0.1 deg (∼3 pc at a distance of 1.7 kpc), with the white circle. The trajectories of the O star runaways are shown in blue, while the B or later-type runaways are shown in orange. The ellipse shapes give the uncertainty on the location for the given trace-back time, which places the runaways near the centre of NGC 6611. The kinematic properties of these runaways are summarised in Table 2.

We remain with 8 early type runaways, which we consider as our final runaway sample. These 8 runaways include 3 O stars and 2 B stars with a spectral classification. Based on the Gaia and 2MASS magnitude and distance from their individual parallax, we estimate the spectral type of an early O star, an early B star and a late B/A star for the 3 unclassified runaways. Two runaway O stars, BD−14° 5040 and UCAC2 27149134 were also found by Gvaramadze & Bomans (2008) and are known to produce bowshocks, which we show in Appendix A.

We show the trajectories of the 8 runaways with respect to NGC 6611 in the right panel of Fig. 7, with a trace-back time that places them near the centre of NGC 6611. We show for each of the runaways their 2σ error ellipse on the position calculated from the position, proper motion and corresponding covariance matrix. Figure 8 shows the runaways far outside NGC 6611 in the AllWISE W4 (22 μm) image in Galactic coordinates, adopting similar colouring as before. The two bowshock producing O stars BD−14° 5040 and UCAC2 27149134 are shown with the blue tracks and are consistent with coming from the same position in NGC 6611 ∼1.08 Myr ago and move in almost exactly opposite directions. Two other runaways found far outside of M16 are shown in orange. The found early type runaways have kinematic ages ≲1.6 Myr, consistent with the age of the young population. Of the 8 runaways found, 4 have Δv_T > 20 km s⁻¹ and 4 have Δv_T < 20 km s⁻¹. We summarise the kinematic properties of the runaways in Table 2 and give the identifier, spectral type and kinematic properties. We note that Δv_T for the runaways is determined with the individual parallax.

Fig. 8. AllWISE W4 (22 μm) grey-scale image of M16 and the surrounding region in Galactic coordinates. We show the current position of the runaways found outside M16 also shown in Fig. 7 adopting similar colours. We trace back the motion of these runaways relative to the cluster. UCAC2 27149134 and BD–14° 5040 have a clear bowshock which we show in Appendix A. We show for context the position of the O stars with the red crosses and indicate the position of NGC6604 and M17.

We have increased the trace-back time to find runaways consistent with coming from the older population with an age of 7.5 ± 0.4 Myr. However, we find no clear candidates coming from the older population in the last 10 Myr. There are several limitations to this, as the Galactic potential may not be ignored over such long time scales, and the centre and mean proper motion of the older population may be different.

7. The O stars in NGC 6611

For a more complete view on the importance of dynamical interactions for massive stars in NGC 6611, we turn to the O stars. The majority of the O stars in NGC 6611 have been identified and spectrally classified (Evans et al. 2005; Sana et al. 2009; Sota et al. 2011, 2014; Maíz Apellániz et al. 2016). The binary properties for part of the O stars have been well studied through radial velocity variations in optical spectra (Sana et al. 2009, 2012). Systemic radial velocities have been determined for several binary O star systems, allowing us to study the transverse and radial motion for most of the single and binary O star systems. The systemic radial velocities have been collected from the literature where possible. We list the O stars, along with the found runaways, in Table 2 and give their spectral type. We also give the absolute Δv_T and calculate the uncertainty using the individual parallax and assuming Gaussian uncertainties on the proper motion of the individual system and NGC 6611.

The relative transverse motion (Δv_T) of the O stars in the vicinity of NGC 6611 is shown in the left panel of Fig. 9. The length of each arrow is proportional to the Δv_T of the star. Two runaway O stars are missing as they are outside this image. In the right panel of Fig. 9 we show the relative radial velocity Δv_R as a function of their |Δv_T| with respect to the average of NGC 6611, coloured and marked similarly as in the left panel. The grey dashed lines and shaded regions indicate respectively the mean and standard deviation on Δv_R determined in Sect. 3. For |Δv_T|, we indicate the median outward velocity from Kuhn et al. (2019) with the grey shaded region. The runaway lower limit on |Δv_T| of 3 km s⁻¹ is shown with the green dash-dotted line and is comparable to the escape velocity for a typical runaway. The three O star systems BD–13° 4923, BD–13° 4929 and HD 168076, have poor quality astrometry for which we only show the (systemic) radial velocity on the left side. We neglect the expansion velocity in the radial direction, since the uncertainty on the average cluster radial velocity is dominated by the formal error of 3.7 km s⁻¹, compared to the average expansion velocity of ∼0.9 km s⁻¹ in the transverse direction.

Fig. 9. Kinematic properties of the O stars in the vicinity of NGC 6611. Left: position and relative proper motion of the O stars near NGC 6611 overlapped on the DSS2 B, R and I colour image. The O stars shown in blue and cyan have velocities comparable or greater than the escape velocity and move away from the centre of NGC 6611. The O stars shown in magenta do not have this clear transverse motion away from the centre of NGC 6611, but could have had dynamical interactions in the past (see Appendix B. The O stars in red have no evidence for dynamical interactions in the past. The angular distance equivalent to 3 pc is given for a distance of 1.70 kpc. We highlight the centre of NGC 6611 with the dashed white circle (r = 0.05 deg). Right: relative (systemic) radial velocity ΔvR as a function of the absolute relative transverse velocity |ΔvT| for the O stars in NGC 6611, coloured similarly to the left panel. The O stars for which we have not used the Gaia astrometry (due to ruwe > 1.4) and thus have unknown |ΔvT| are shown on the left. The uncertainty on the radial velocity of NGC 6611 determined from the B0 to B4 stars is shown with the horizontal black shaded region. The median outward velocity determined by Kuhn et al. (2019) for NGC 6611 is shown with the vertical black shaded region.

Several other O stars near the centre of NGC 6611 show a receding motion away from the centre of NGC 6611, highlighted in cyan. This includes the two runaway O stars already identified in Sect. 6, BD–13°4927 and LS IV –13 14. On top of this, two other O star systems show a similar receding motion. These are the O7 V + O8 V binary system HD 168137 with an highly eccentric orbit (e = 0.902 ± 0.058; Sana et al. 2009; Le Bouquin et al. 2017) and the O8.5 V star W161. Both O stars might only be considered to move away marginally significant (2.6 ± 0.2 km s⁻¹ and 2.8 ± 0.3 km s⁻¹) compared to the 3 km s⁻¹ adopted. Complementing the transverse velocities, the (systemic) radial velocities also indicate differential motion at the 2σ level. Specifically HD 168137 could turn out to be a runaway moving mostly in the radial direction, if the systemic radial velocity can be constrained better. This could for example be visible as a bowshock in the infrared with Spitzer, where strong emission is seen ∼25″ from HD 168137 (see Appendix A; Pilbratt et al. 1998; Flagey et al. 2011). We have therefore included both O star systems as accelerated stars.

HD 168504 and HD 168183 have RUWE > 1.4 (1.670 and 1.415 respectively) and were therefore excluded from the runaway search. These two systems have parallaxes consistent (within 2 and 1σ, respectively) with that of NGC 6611 (see Fig. 2). This might indicate that the astrometric solution and thus the proper motion may not be of poor quality (see e.g. Maíz Apellániz et al. 2021). HD 168504 is an O7.5V(n)z star; HD 168183 is an O9.5 III star in an eclipsing binary system with a mid-type B star (Sana et al. 2009; Maíz Apellániz et al. 2016). These two systems are shown with the blue squares in Fig. 9 and have |Δv_T| > 4 km s⁻¹. While we could attribute this significantly deviating Δv_T to a poor astrometric solution, the systemic radial velocity of HD 168504 is significantly deviating (> 3σ) from the average of NGC 6611. HD 168504 and HD 168183 are also located several tenths of degrees from the centre of NGC 6611 and show a clear receding motion. The combination of their location, transverse and radial velocities make these two systems clear outliers with respect to the cluster average and we have included them as runaways.

This conclusion is not easily made for BD–13° 4929. While the parallax is also consistent within 1σ from NGC 6611 with ruwe > 1.4, the transverse velocity does not clearly point away from the centre of NGC 6611 but instead from the northern part of the nebula (not shown). BD–13° 4929 is a hierarchical triple system where the observed O7 V star is in a wider orbit around two early type B stars in an inner orbit. The average radial velocity of BD–13° 4929 is 17 ± 1 km s⁻¹ (determined from the O star); this is not the systemic radial velocity and could vary over a longer time scale (Sana et al. 2009, 2012). We therefore include BD–13° 4929 as a ‘candidate dynamically interacted O star’ in Table 2. For these O stars, we can not draw conclusions on whether these stars have dynamically interacted, but this remains a possibility. See for example Sect. 8.4 on why HD 168076 could have participated in dynamical interactions.

Other O star systems show no evidence for significantly deviating transverse or (systemic) radial velocities. These systems are BD–13° 4923, BD–13° 4928, BD–13° 4930, HD 168075, HD 168076, W222 and W584, for which both the (systemic) radial velocity and |Δv_T| are consistent (or unknown) with NGC 6611. These systems are highlighted in red and magenta. Sana et al. (2009) report systemic radial velocities for BD–13° 4923 of 30.7 ± 4.3 and 3.4 ± 5.3 km s⁻¹ for the primary O4 V((f+)) star and secondary O7.5 V star respectively. A significant difference is present between these two velocities, while they should be equal. The spectral line He II λ4686 used for the orbital solution can be affected by a stellar wind for early type O stars (Sana et al. 2013). Here, we will assume the 3.4 ± 5.3 km s⁻¹ belonging to the O7.5 V star to be the systemic radial velocity, which is consistent with the average radial velocity of NGC 6611. The binary orbit for HD 168075 is solved, but the systemic radial velocity is not reported in the literature. We will adopt a conservative 17 ± 10 km s⁻¹ based on Fig. 1 in Barbá et al. (2010). W584 is reported to have a bowshock in the infrared as observed with Spitzer (Guarcello et al. 2010). However, the transverse and radial velocity are in reasonable agreement with the average of NGC 6611 and we have not included this star as a runaway.

Next, we will quantify the time scale over which the O stars and runaways have been accelerated. In order to do this, we trace back the runaways and determine how long ago these stars were near the centre of NGC 6611. Specifically, we calculate the range in time the runaways spent within a radius r_NGC 6611 = 0.05 deg from the centre coordinates (α_NGC 6611, δ_NGC 6611) = (274.67 deg, –13.78 deg). This radius is the equivalent of 1.5 pc at distance of ∼1.70 kpc, but we note that this radius has been conservatively chosen and that the closest approach of almost all runaway stars is within 1.0 pc. We show this radius with the white dashed circle in the left panel of Fig. 9. The only exception is HD 168183, for which we use a radius of 0.08 deg.

Since both the proper motion of the star and NGC 6611 have uncertainties, we resort to Monte-Carlo simulations. We neglect the radial motion and uncertainty on the initial position. We assume the proper motion of the star and NGC 6611 to follow Gaussian distributions with means equal to the observed proper motions and standard deviations equal to the observed proper motion errors (as determined in Sect. 3 for NGC 6611). In each iteration, we draw proper motions following these Gaussian distributions for the star and NGC 6611. We consequently calculate the times when the star enters and leaves the circle with r_NGC 6611 = 0.05 deg for each iteration, giving two respective distributions over all iterations. We determine the time when the star enters and leaves this circle as the 14th and 86th percentile over their respective distribution from 1000 iterations, which respectively give the lower and upper bound on the kinematic age t_kin. We list the range in t_kin for each accelerated star in Table 2.

To visualise the kinematic ages, we show them as a function of |Δv_T| for each star in Fig. 10 coloured similarly as in Fig. 9. The orange data indicate the B or later-type runaways. The kinematic ages are better constrained for stars with higher relative transverse velocities, as they spent less time in NGC 6611 when traced back, as expected. We also show the age and uncertainty of the young population, 1.3 ± 0.2 Myr, with the dashed black line and shaded region, respectively. The kinematic ages of the runaways with |Δv_T| ≳ 20 km s⁻¹ are in excellent agreement with this isochronal age. The kinematic ages of the accelerated stars with |Δv_T| ≲ 20 km s⁻¹ span a larger range between ∼ 0.0 and 2.0 Myr. This shows that runaways might not have been ejected instantaneously and that ejections could even happen one Myr after the formation of the cluster.

Fig. 10. Relative transverse velocity as a function of the determined range in kinematic age. We show the best-fit age and uncertainty determined with isochrone fitting with the black dashed line and shaded region respectively. The runaway O stars are shown in blue and cyan, while the runaway B or later-type stars are shown in orange. The O stars are coloured similarly as in Fig. 9.

8. Discussion

We have studied the young massive cluster NGC 6611 in the Eagle Nebula (M16). NGC 6611 hosts a rich early type stellar population and formed at least 19 O stars. From the CAMD, we find evidence for two populations of stars. The younger population has an average A_V = 3.5 ± 0.1 mag and an age of 1.3 ± 0.2 Myr. The older population has an average A_V = 2.0 ± 0.1 mag and an age of 7.5 ± 0.4 Myr. These two populations show a clear difference in spatial distribution, with the younger population concentrated towards the centre of NGC 6611, while the older population is spatially more extended.

We have searched for runaways with |Δv_T|> 3 km s⁻¹ and find 4 O stars and 4 B or later type stars which have been ejected in the last 1.6 Myr. Dynamical interactions between stars in the younger population are the prime mechanism to explain these runaways. We have studied the kinematics of the O stars in NGC 6611 and show how at least nine out of 19 O stars have velocities comparable to or greater than the escape velocity of the cluster (∼3 km s⁻¹). These runaway O stars can all be traced back to the centre region of NGC 6611 and have a kinematic age between 0.0 and 2.0 Myr, which is within the uncertainties less than or comparable to the age of the younger population.

8.1. Physical parameters of NGC 6611

Before further discussing our results, we will consider some of the key physical parameters of NGC 6611. These parameters include the astrometry, (initial) mass function, mass, radius, stellar multiplicity and age.

We summarise the astrometric parameters of NGC 6611 in Table 3, where we give the centre coordinates, proper motion, distance and radial velocity determined in Sect. 3. We also list σ_1D and the median v_out adopted from Kuhn et al. (2019).

We will assume a Salpeter (Salpeter 1955) or Kroupa (Kroupa 2001) distribution for the initial masses. The Kroupa initial mass function is also used to define the stellar occupation along the isochrones. With the initial mass function, we can roughly estimate the mass of NGC 6611. If we assume the sample of 19 O stars (with a conservative M > 18 M_⊙) to be mostly complete, we determine the mass of NGC 6611 to be ∼8.9 × 10³ or 5.4 × 10³ M_⊙ for a Salpeter and Kroupa initial mass function, respectively, between 0.08 and 100 M_⊙. This mass is still uncertain as we are dealing with small-number statistics, possibly excluding unclassified stars and undetected O stars. If the 7 O stars which have not been investigated for the prescence of a companion have an O star companion, we determine a mass of ∼7–12 × 10³ M_⊙. Similarly, including all spectrally classified B0-2 stars consistent with the distance of NGC 6611 yields a mass that agrees with the aforementioned values (assuming a 50% binary fraction; Banyard et al. 2021). The mass of NGC 6611 is estimated in the literature at ∼2 to 25 × 10³ M_⊙, taking into account the ∼170 early type B stars (Bonatto & Santos 2006; Wolff et al. 2007; Portegies Zwart et al. 2010). We note that if the mass is 20 × 10³ M_⊙, we would expect to find ≳67 O stars, for which we find no evidence. We will consider a mass in the range of 5–10 × 10³ M_⊙, which is mostly sensitive to the adopted initial mass function. Estimating the mass for the older population is problematic. We have found 21 stars (22 with the B2.5 I) belonging to the old population, which brings large uncertainties due to small number statistics. On top of this, the older population is spatially more extended, which could result in a bias against identification of cluster membership. In summary, we can not give a reliable estimate for the mass of the older population and leave this for future work.

Several definitions exist for the ‘radius’ of a cluster, including the half-mass radius, effective radius, virial radius and core radius (Portegies Zwart et al. 2010). We will use the half-mass radius r_hm here. Using all spectrally classified stars consistent with the distance of NGC 6611 we determine the current projected r_hm(proj) to be ∼2.0 pc from a weighted average, where the weights are taken to be the mass from Pecaut & Mamajek (2013). Using different sub-samples, such as the O stars or B stars, gives similar results in the range of 2–2.5 pc for r_hm(proj). Estimates for the half-light radius, taking into account the young population of pre-main-sequence stars, similarly give results in the range of r_hm(proj) of 2–2.5 pc. This is the projected 2D radius. Assuming for simplicity a Plummer density distribution (Plummer 1911) this corresponds to a 3D r_hm of 0.77r_hm(proj) ∼ 1.5–2.0 pc.

Most massive stars are in binaries or higher-order systems. This is also the case for NGC 6611, where 5 out of 9 systems studied in depth for the prescence for a companion are binary or higher-order systems (Sana et al. 2009, 2014). We should therefore consider binary fraction percentages between 50–100%. While Table 2 contains 10 spectrally classified ‘single’ O stars, 7 of these have not been investigated for the prescence of a companion and could still hide one. HD 168504 is an O7.5V(n)z star and is presumed to be single, although the ruwe > 1.4 gives a first indication of a possible binary system (Maíz Apellániz et al. 2016; Belokurov et al. 2020). Gaia indicates that W161 is separated from another source by only ∼3″. W161 could possibly be a binary with the radial velocity measured by Evans et al. (2005) and Martayan et al. (2008) deviating by 16 km s⁻¹, although consistent within 2σ considering the uncertainties. W161 is also relatively bright in the near-infrared for an O8.5 V star, with a 2MASS K_s magnitude brighter than the O7 V + O8 V system HD 168137, despite being more reddened. Similarly W584 located in the north-west part of the nebula also shows another source separated by ∼2″ and is only ∼2.5 mag fainter in the G-band. In summary, the binary properties of the O stars in NGC 6611 have been well studied for the apparently brightest O stars, mostly located near the centre of NGC 6611, but remains incomplete when considering the late type and runaway O stars. Thus, in order to properly determine the high-mass end of the IMF, the OB runaways have to be accounted for.

8.2. The ages of the two populations of NGC 6611

The CAMD has revealed two populations of stars, specifically for stars fainter than M_G ≳ 3 mag. Not only are the estimated ages of the two populations different with 1.3 ± 0.2 Myr compared to 7.5 ± 0.4 Myr for the young and old population, also the average extinction is different. The younger and older population require on average 3.5 ± 0.1 mag and 2.0 ± 0.1 mag extinction in the V-band. This could imply that the older population is located in front of the younger population. If so, the difference in distance would be small (≲ 20 pc), as we find no deviations in the parallax, not taking into account that we have partially selected members from the parallax.

Another difference between these two population is their spatial distribution. The young population has a half-light radius of ∼2 pc, while the old population has a half-light radius of ∼3.5 pc, corresponding to a 3D r_hm (assuming the half-light and half-mass radius are similar) of 1.5 and 2.7 pc, respectively. The fact that the size of the older population is larger than that of the younger population agrees with the expectation that the clusters are still expanding.

Assuming a cluster mass of 7.5 × 10³ M_⊙ and a mean stellar mass of 0.5 M_⊙ we derive half-mass relaxation times on the order of ∼100 Myr, the time it takes for a cluster to reach virial equilibrium. This is much longer than the age of both populations, so the cluster is not yet in hydrodynamical equilibrium (unless the initial r_hm of the cluster was much smaller than the present radius and on the order of ∼0.2 pc). From Fig. 6 the increase in the projected half light radius is about ∼0.8 pc Myr⁻¹ and 0.3 pc Myr⁻¹ for the young and old population respectively, corresponding to an expansion of ∼0.8 and 0.3 km s⁻¹, respectively. This expansion velocity for the young population is consistent with that measured by Kuhn et al. (2019) of 0.90 ± 0.23 km s⁻¹. We can also estimate the dynamical timescale t_dyn, the time it takes a typical star to cross the cluster. With a cluster mass of 7.5 × 10³ M_⊙ and 3D r_hm of 1.5 and 2.7 pc, we find t_dyn ∼ 0.4 and ∼0.7 Myr for the young and old population respectively. The t_dyn of 0.4 Myr for the young population is similar to the time of minimum radius found in Fig. 6. This indicates that we should not blindly trust the radius given in Fig. 6 at a time of –0.4 Myr or earlier. The motions of the stars could have changed direction or crossed the cluster on longer timescales. The t_dyn of 0.7 Myr for the old population differs from the time of minimum radius (0.4 Myr), but this is not surprising considering the uncertainties involved as mentioned in Sect. 8.1. The calculation for t_dyn also assumes a virialised cluster, which may not be the case here.

The age of the young population in NGC 6611 can also be derived from the astrometry. The kinematics of the runaways and O stars show that dynamical interactions happened over the last ∼2 Myr. The four runaways with |Δv_T|> 20 km s⁻¹ specifically have been kicked out 1.0 to 1.6 Myr ago. This is in excellent agreement with the isochronal age of the young population. We expect the dynamical interactions to take place right after star formation (Fujii & Portegies Zwart 2011).

For the stars not included in the blue or red group, only one star shows evidence of belonging to the old population. The B2.5 supergiant BD–13° 4912 shows a clear deviation from the young population in Fig. 5. BD–13° 4912 has been used as an argument in favour of an older age or noncoevality for NGC 6611 as a whole, with its evolutionary status implying an age of ∼6 Myr (Hillenbrand et al. 1993; Gvaramadze & Bomans 2008). Similar to Hillenbrand et al. (1993), we find no evidence for other bright stars belonging to the old population. The O stars earlier than O7.5-8.5 should have already experienced their supernova or are considerably evolved if they were to belong to the old population. The O and B stars are almost all spectrally classified as dwarfs and show a clear spatial clustering towards the centre of NGC 6611, similar to the young population. We therefore deem it justified to consider the majority of the early type stars as part of the young population.

The presence of an older population in NGC 6611 has been suggested before. Guarcello et al. (2007b); Guarcello et al. (2007a, 2009, 2010) report an over-density of blue stars in the optical which also display a near-infrared excess indicative of circumstellar discs. These ‘blue stars with excess’ (similar to our ‘blue group’ but located at M_G ∼ 7 mag) were initially attributed to stars in the young population possessing a circumstellar disc. When viewing these stars at a certain inclination, the scattering of light from the disc could make these stars appear bluer (Guarcello et al. 2010; Bonito et al. 2013). However, De Marchi et al. (2013) show that only a few percent of the stars in the young population are subject to the aforementioned effect and that most of the blue stars with excess indeed make up an older population with an uncertain age of 8 to 32 Myr. Other evidence supporting the presence of an older population comes from the X-ray luminosity, lithium absorption lines and their optical magnitudes (De Marchi et al. 2013). De Marchi et al. (2013) also find a significantly different spatial distribution for the older population, consistent with our spatially extended distribution for the blue group in Fig. 6.

Several theories exist for the formation of a cluster with multiple populations of stars. For example, the Orion nebula cluster is now observed to have three different populations of stars differing in age by ∼0.5 to 1 Myr each (Beccari et al. 2017; Jerabkova et al. 2019). If the first generation of stars (observed as the oldest) formed several OB stars, these stars could halt star formation through their ionising photons and strong stellar winds. If these OB stars are all ejected by dynamical interactions, the remaining gas in the molecular cloud can be allowed to collapse again and form a second generation, also including O stars (Kroupa et al. 2018). This is likely not the case for NGC 6611, since it would require the two populations to differ in age by ≲1 Myr, while we find ≳6 Myr. The younger population also appears to have formed significantly more (O) stars than the older population, while this is opposite to what is observed in the Orion Nebula Cluster (Jerabkova et al. 2019).

Another scenario is that a nearby older population of stars was captured by the younger population as it was experiencing cloud-collapse (Pflamm-Altenburg & Kroupa 2007). Stars present in the field, which could have been formed by earlier star formation, could be gravitationally bound to the molecular cloud that formed NGC 6611 if their paths cross. The Sagittarius Spiral Arm has seen star formation sites (see Fig. 8), with M17, M16 and NGC 6604 visible as the giant H II regions to mention a few. The capture of an older less massive cluster should therefore not be excluded. Similarly, if the older population passed by a different molecular cloud, a supernova could have triggered the molecular cloud to collapse and form the younger population (Guo et al., in prep.). The difference in ages between the young and old population of ∼6 to 7 Myr would be consistent with the core-collapse supernova of a massive star initially belonging to the old population. More detailed analysis is needed on the older population before conclusions can be drawn on the formation mechanism of the younger population.

8.3. Dynamical interactions in NGC 6611

We now turn to the kinematics of the NGC 6611 runaways and O stars. We have found a total of 9 O and 4 B or later-type stars with velocities comparable to or greater than the escape velocity for NGC 6611 (|Δv_T| ≳ 3 km s⁻¹). Since the young population is < 2 Myr old, we do not expect core-collapse supernova to have happened yet, making the dynamical ejection scenario the only feasible mechanism. The runaways coming from the centre of NGC 6611 are in line with this assumption as dynamical interactions are expected to be most prominent in the densest regions of young star clusters (Fujii & Portegies Zwart 2011). Out of the remaining 10 O stars, 5 O stars could have had dynamical interactions in the past based on their kinematic properties, but this remains inconclusive. The 5 O stars for which we find no evidence for acceleration or dynamical interactions contain three late type O stars (O9 V, O9.5 Vp and O9.5 V) and an O4 V((f+)) + O7.5 V binary system for which we have poor quality astrometry. If we assume that the runaway O stars were created through dynamical interactions, we should consider that roughly ≳50% of the O stars underwent dynamical interactions.

It is unclear if these percentages are to be expected for clusters such as NGC 6611. Simulations and models of young massive clusters have been performed to study the physical properties of runaways and the dynamical evolution of clusters (Fujii & Portegies Zwart 2011, 2014; Oh & Kroupa 2012, 2016; Oh et al. 2015). We have found two runaways with M ≳ 15 M_⊙, which is consistent with the results of Fujii & Portegies Zwart (2011). For example, their Fig. 3 predicts a runaway fraction of 0.05–0.2 for a cluster mass of 6 × 10³ M_⊙ compared to our runaway fraction of ∼0.1 (v_T ≳ 20 km s⁻¹). Our results are also in agreement with their Fig. 6, where we find that runaways with relatively low |Δv_T| are more common than runaways with relatively high |Δv_T|. This is similar to what has been found in Oh et al. (2015) and Oh & Kroupa (2016), where a large fraction of the O stars could be ejected at low velocity (≲10 km s⁻¹) under certain conditions.

This creates a consistent picture so far; however, one important parameter in these models is the initial half-mass radius r_hm. The initial r_hm typically assumed in these models is in the range of 0.1–0.5 pc. NGC 6611 is a relatively dispersed cluster for which we determine a present-day r_hm ≃ 1.5–2 pc. We visualise radii of 0.1, 0.5 and 1.5 pc in Appendix B. Even when considering the expansion of a cluster like NGC 6611 and dynamical interactions contributing to this, it is difficult to reconcile this with a small initial r_hm. The O stars can not be traced back to such a small radius, which is similar for the fainter stars. A radius of 0.1–0.5 pc is more typical of dense clusters such as Westerlund 2 and NGC 3603 and it is difficult to fit NGC 6611 among these clusters (see e.g. Pfalzner 2009; Portegies Zwart et al. 2010; Pfalzner & Kaczmarek 2013).

What so far has not been considered in models of dynamical evolution of clusters is that the massive binary system population initially had wider orbits than what is currently observed. Ramírez-Tannus et al. (2021) find that the massive binaries harden over time, implying that their orbits shrink as part of the (late) formation process. The time scale over which these binaries harden is ∼1 to 2 Myr. If the massive binaries indeed form in wider orbits, it would increase the cross-section for dynamical interactions. The orbital period distribution considered in the simulations of young massive clusters often follows the Sana et al. (2012) distribution, which has been derived for clusters that are 2–5 Myr old.

Westerlund 2, one of the clusters included in Ramírez-Tannus et al. (2021), shows a rich history of high velocity ejections (≳25 km s⁻¹) in the last ∼1 Myr, with likely 7 O stars being thrown out (Drew et al. 2018). NGC 3603 shows a similar behaviour, with 11 O star ejections at runaway velocity (Drew et al. 2019). Since both have estimated ages of ≲2 Myr, only the dynamical ejection scenario is expected to contribute to the production of runaways. For Westerlund 2, the percentage of O stars observed as runaways is ∼25% (Drew et al. 2018). Considering that we see more low velocity runaway O stars, Westerlund 2 and NGC 3603 may have also produced a significant number of low velocity runaway O stars, depending on their respective escape velocity.

8.4. Bully binary in NGC 6611

HD 168076 is the brightest visual star in NGC 6611/M16 and is spectrally classified as an O4 IV(f) star (Sota et al. 2011; Maíz Apellániz et al. 2016). Bosch et al. (1999) show a bright companion separated 0.15″ from the primary star. Sana et al. (2014) also report the prescence of a secondary at this separation. Sana et al. (2009) studied the multiplicity properties of HD 168076 with high-resolution optical spectra and find no significant radial velocity variations in line with a large binary orbital period. They also assigned a tentative spectral classification of O7.5 V / O9 III to the secondary.

This makes HD 168076 a massive binary system possibly with a wide binary orbit. Such a massive system could be responsible for a major part of the ejections of other O stars, the so-called ‘bully binary’ scenario (Fujii & Portegies Zwart 2011). Unfortunately, the astrometry of HD 168076 is of poor quality with ruwe = 8.382, and makes it impossible for us to trace its position back in time. If the lower limit on the binary orbital period of P_orb ≳ 10⁵ days reported by Sana et al. (2009) is correct, it would suggest that HD 168076 even now has a high cross-section for dynamical interactions (see e.g. Fig. 2 in Fujii & Portegies Zwart 2011). A recent role of HD 168076 in dynamical ejections could possibly be seen in the acceleration of the two O stars BD−13° 4927 and LS IV –13 14, which we show in Appendix B. If HD 168076 was the culprit of the accelerations, t_kin would be ≲0.3 Myr for both O stars. More investigation is needed on the 3D velocities of HD 168076, LS IV –13 14 and BD−13° 4927 to draw further conclusions, so it remains inconclusive as to what role HD 168076 played in accelerating the O stars.

8.5. Runaways

The two runaway O stars BD–14° 5040 and UCAC2 27149134 were also discussed in Gvaramadze & Bomans (2008), initially discovered by their bowshock. Relative to NGC 6611 the two O stars move in exactly opposite directions and are consistent with coming from the same position near the centre of NGC 6611. This could suggest that one of the binaries initially formed in NGC 6611 has been disrupted through binary-binary interactions.

This is similar to the well-known example of the two runaway stars μ Col and AE Aur which are now thought to have been ejected ∼2.5 Myr ago in or near the vicinity of the Orion Nebula Cluster (Blaauw & Morgan 1954; Hoogerwerf et al. 2001). Gualandris et al. (2004) showed that a binary-binary interaction involving μ Col, AE Aur and the binary system ι Ori could lead to the ejection of the two runaway stars.

BD–14° 5040 and UCAC2 27149134 may have been subject to a similar binary-binary interaction, though we find no evidence for other runaways consistent with coming from this event. We could therefore assume the other stars involved to remain in the vicinity of NGC 6611, similar to the case of ι Ori.

Since we know the spectral type for BD–14° 5040, an O5.5 V star, a rough estimate of the mass is ∼34 M_⊙ (Martins et al. 2005). To conserve momentum among the two stars, UCAC2 27149134 needs to have a mass ∼74 M_⊙. We note that including radial velocities in this calculation would only change the resulting mass by ∼10% assuming BD–14° 5040 moves towards us with −100 km s⁻¹. Even then, a conservative lower limit on the mass would be 60–65 M_⊙ for UCAC2 27149134.

Similar to Gvaramadze & Bomans (2008), we can roughly estimate the spectral type from the 2MASS JHK_s magnitudes and the Gaia EDR3 parallaxes. UCAC2 27149134 has K_s = 7.396 ± 0.023 mag, (J – H) = 0.41 ± 0.06 mag and (H – K_s) = 0.288 ± 0.06 mag. The distance to UCAC2 27149134 is ${1750}_{-60}^{+75}$ pc from the individual parallax. With R_V = 3.1, we estimate A_V ∼ 5.7 mag. These values agree with Gvaramadze & Bomans (2008), however Martins & Plez (2006) do not provide a photometric calibration for the 2MASS passbands, which is noticeable when comparing to Pecaut & Mamajek (2013). For UCAC2 27149134, we estimate a photometric spectral type of O3-4 V, with the dwarf luminosity classification based on the prescence of an O3.5 V((f+)) and O4 V((f+)) star in NGC 6611. A spectral type of O3 V would give a rough estimate of ∼60 M_⊙ on the mass, in reach of the 60–65 M_⊙ ‘lower limit’ mentioned earlier. Spectroscopic follow-up on this source may confirm our hypothesis.

Both UCAC2 27149134 and BD–14° 5040 produce bowshocks visible in the AllWISE data (see Appendix A). A third bowshock producing runaway O star analysed in Gvaramadze & Bomans (2008), HD 165319, is spectrally classified as an O9.7 Ib star. This star was identified as a runaway coming from NGC 6611, for which we find no evidence. HD 165319 is located at a distance of ∼1.4 kpc and would be even closer in the past considering its radial velocity of ∼25 km s⁻¹. It is thus inconsistent with being at the same distance as NGC 6611. The proper motion of HD 165319 indicates that its projected closest approach is about ∼0.5 deg away from the centre of NGC 6611 and requires a kinematic age of ∼4 to 5 Myr, again arguing against an origin in NGC 6611.

9. Summary and conclusion

We have studied the Eagle Nebula and its young massive cluster NGC 6611 in detail with Gaia EDR3. We summarise our findings below:

• We have confirmed 137 members of NGC 6611, at a distance of 1706 ± 8 pc.

• NGC 6611 includes two populations of stars. We have determined the mean age and extinction with isochrone fitting. The young population has an age of 1.3 ± 0.2 Myr and a mean extinction of A_V = 3.5 ± 0.1 mag. The old population has an age of 7.5 ± 0.4 Myr and a mean extinction of A_V = 2.0 ± 0.1 mag. Almost all OB stars are most likely belonging to the young population. Only the B2.5 I star BD–13° 4912 may be associated with the old population.

• We have found 8 runaways (|Δv_T|> 3 km s⁻¹) coming from NGC 6611, of which 4 are O stars and 4 are B or later-type stars. These runaways have kinematic ages consistent with being ejected from the young population.

• Investigating the kinematics of all O stars present in NGC 6611 in detail shows that at least 9 out of 19 O stars have velocities comparable to or greater than the escape velocity of the cluster (∼3 km s⁻¹). With the dynamical interactions as the sole ejection mechanism given the young age of the cluster, we show that the accelerated O stars and runaways have kinematic ages in agreement with the young population.

• If massive binaries are initially formed with wider separations than what is currently observed, it will increase the cross-section for dynamical interactions. The percentage of escaping O stars in NGC 6611 (≳45%) is difficult to reconcile with n-body simulations of young massive clusters. If these simulations account for an initially wide massive binary systems, it may bring their runaway fraction in line with our results.

The number of young massive clusters younger than ∼1 Myr is extremely limited. Runaway stars can give valuable insight into the early evolution of young massive clusters and massive binaries as shown here. Extending these results to other young massive clusters can give a larger sample-size to study the importance of dynamical interactions and the initial massive binary population.

Acknowledgments

The data and data products presented in this paper are available at the following DOI: http://doi.org/10.5281/zenodo.6372280. MS acknowledges support from NOVA. SR acknowledges funding from the European Research Council Horizon 2020 research and innovation programme (Grant No. 833925, project STAREX). This research has made use of NASA’s Astrophysics Data System, the SIMBAD database, operated at CDS, Strasbourg, France, the “Aladin sky atlas” developed at CDS, Strasbourg Observatory, France. 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.

References

Appendix A: Bowshocks of runaways

The runaway O stars which have visible bowshocks in the AllWISE imaging are UCAC2 27149134 shown in Figure A.1 and BD–14° 5040 in Figure A.2. W584 and HD 168137 are reported to have possible bowshocks in the Spitzer images which we show in Figure A.4 and Figure A.3.

Fig. A.1. AllWISE (W1, W2 and W4 as blue, green and red, respectively) false-colour image of the bowshock produced by UCAC2 27149134 (18′ by 18′). The blue circle gives the position of the star. The red arrow gives the absolute transverse motion of the star, while the blue arrow gives the transverse motion relative to the cluster.

Fig. A.2. AllWISE (W1, W2 and W4 as blue, green and red, respectively) false-colour image of the bowshock produced by BD–14° 5040 (18′ by 18′). The blue circle gives the position of the star. The red arrow gives the absolute transverse motion of the star, while the blue arrow gives the transverse motion relative to the cluster.

Fig. A.3. Spitzer IRAC 8.0 μm image of W584 (2.8′ by 2.8′) The blue circle gives the position of the star. The red arrow gives the absolute transverse motion of the star, while the blue arrow gives the transverse motion relative to the cluster.

Fig. A.4. Spitzer IRAC 8.0 μm image of HD 168137 (2.8′ by 2.8′) The blue circle gives the position of the star. The red arrow gives the absolute transverse motion of the star, while the blue arrow gives the transverse motion relative to the cluster.

Appendix B: The centre of NGC 6611

The runaway O stars can all be traced back to the centre part of NGC 6611. If we assume dynamical interactions to be the cause of their ejection, we might be able to identify the exact location where these interactions took place. This is best done for the O stars most recently ejected, since the uncertainties increase over the time traced back. To do this we show the Hubble Space Telescope image of NGC 6611 with the Advanced Camera for Surveys image on top of the ESO B, V and R composite mosaic image in Figure B.1. We show the location and relative proper motion traced back for the O stars in this image for t = –0.13 Myr, coloured and marked similarly as in Figure 9.

Fig. B.1. B, V and R composite mosaic image of NGC 6611 (ESO, press release 0926). We overplot the near-infrared Hubble Space Telescope image of the centre of NGC 6611 (ESA/Hubble). The purple cross indicates the centre of NGC 6611, with three white dashed circles with radii of 0.1, 0.5 and 1.5 pc at a distance of 1.70 kpc. We show the motion of the O stars and their position at t = –130 kyr relative to NGC 6611, coloured similarly as in Figure 9, with the error ellipse giving the uncertainty on the proper motion.

One O star system stands out, the bright system HD 168076 composed of an O3.5 V((f+)) + O7.5 V. HD 168076 fits the description of a bully binary (Fujii & Portegies Zwart 2011), being able to accelerate other O stars (see Section 8.4 for a discussion). While the proper motion of HD 168076 is unusable, both HD 168075 and LS IV –13 14 are roughly consistent with coming from a similar position. This is also true for the O7 V((f)) star W222, located north of BD–13° 4923 and BD–13° 4927 and not visible in this image. W222 is consistent with being at the same location as W161 for t = –0.22 Myr. These are the only cases where the positions (with uncertainty) overlap with each other in the last 0.5 Myr. This of course does not confirm that these O stars have dynamically interacted in the past, since we are ignoring the most difficult parameter to constrain; the radial depth of the cluster. We refer to HD 168076, HD 168075 and W222 in Table 2 as ’candidate dynamically interacted O stars’.

For other accelerated O star systems it is more difficult to constrain the location of dynamical interaction. BD–13° 4927 does not come across any other O star system when traced back. This is not certain as the local standard of rest of the dynamical interaction could be different and the proper motions of the O star systems HD 168076 and BD–13° 4923 are unknown. The complete 3D velocities for the O stars in NGC 6611 could give more insight on the importance of dynamical interactions in NGC 6611. We have also not included the (early type) B stars, for which (systemic) radial velocities are incomplete.

Appendix C: Members of NGC 6611

All Tables

All Figures

Fig. 1. DSS2 B, R and I colour image of NGC 6611 and the surrounding nebula (M16). The member stars are shown as white plus symbols. North is up and east is to the left. We highlight the O stars in this image and show these with poor quality astrometry with the red squares, while the ones with good quality astrometry are shown with the cyan circles. Two runaway O stars are not visible in this image. The centre of NGC 6611, the position with the highest source density, is indicated with the magenta cross. At the top of the figure the angular distance equivalent to 4 pc is given for a distance of 1.7 kpc. We show the cone-search region with the dashed circle, with a radius equivalent to 0.25 deg or 7.4 pc. The open cluster Trumpler 32 (age ∼ 0.3 Gyr; Kharchenko & Schilbach 1995) can be seen in the top-right of the image.

In the text

Fig. 2. Astrometric parameters of NGC 6611. Top-left: proper motion distribution of the field stars (grey) and members of NGC 6611 (blue). We show the O stars present in NGC 6611 with poor and good quality astrometry with the red squares and cyan circles, respectively. Top-right: parallax distribution of the field stars and members, coloured similarly as in the top-left panel. We also show the fractional parallax uncertainty (ϖ/σϖ) of the members as a function of their parallax. Two O star systems have RUWE = 2.544 and 8.382, and have parallaxes outside this panel, with only part of the error-bar visible. The 1σ error on the parallax of NGC 6611 is shown with the orange bar. Bottom-left: radial velocity of the members in the heliocentric frame collected from the literature. Bottom-right: radial velocity in the heliocentric frame of the B0 to B4 stars after removal of the outliers in the bottom-left panel. Outliers have been excluded as described in Sect. 3. The mean and standard deviation of the radial velocity of these stars are shown with the dashed black line and shaded region, respectively.

In the text

Fig. 3. Colour – absolute magnitude diagram (CAMD) for the average cluster distance for NGC 6611. Errors are smaller than the respective symbol if not shown. Left: Gaia CAMD of the members with known spectral types. The O stars are shown with the open blue squares, B stars with the open orange circles, A stars with the open red stars and unclassified members in black. The two runaway O stars discussed in Sect. 6 are included. We show a reddening line for an O7 V star with RV = 3.5, the black crosses indicate AV = 2.0, 4.0, 6.0 and 8.0 mag from left to right. Right: extinction corrected MV0 for the average cluster distance (RV = 3.56; Kumar et al. 2004) against the intrinsic colour (B − V)0 from the Guarcello et al. (2007a,b) BVI photometry. The stars are marked and coloured similarly as in the left panel. Stars with an uncertain or no spectral type are not shown. We show three PARSEC + COLIBRI isochrones with ages of 1.0, 5.0 and 10.0 Myr with the solid, dashed and dash-dotted grey lines, respectively.

In the text

	Fig. 4. Gaia colour – absolute magnitude diagram. We have separated the early type stars, the blue and red group shown with the black circles, blue triangles and red squares, respectively. Stars with a known spectral type are shown with an open marker.
In the text

Fig. 5. Isochrone fitting results of NGC 6611. Top-left: best-fit isochrone with AV = 2.0 ± 0.1 mag and age of 7.5 ± 0.4 Myr shown with the solid cyan line to the early type stars and the blue group (with and without spectral type). An open symbol represents a known spectral type, while a closed symbol does not. The uncertainty on the age for AV = 2.0 mag is shown with the cyan shaded region. The stars with known spectral type are corrected for extinction according to the best-fit AV with RV = 3.56 (Kumar et al. 2004). Bottom-left: log-likelihood for each isochrone in the searched AV-age grid for the stars indicated in the top-left figure. Each point in the grid is coloured according to its (unnormalised) likelihood value. The log-likelihood gives the quality of the fit, with smaller negative values (in lighter grey) representing better fits. The cyan cross and circle indicate the best-fit AV and age with 1σ uncertainty range. Top-right: similar to the top-left figure, but now for the early type stars and the red group (with and without spectral type). The solid cyan line gives the best-fit isochrone with AV = 3.5 ± 0.1 mag and age of 1.3 ± 0.2 Myr. The uncertainty on the age for AV = 3.5 mag is shown with the cyan shaded region. Bottom-right: similar to the bottom-left figure, but now for the top-right figure.

In the text

Fig. 6. Distribution of the young and old population (the red and blue group, respectively), coloured and marked similarly as in Fig. 4. Top-left: B, V and R composite mosaic image of NGC 6611 (ESO, press release 0926) showing the spatial distribution. Top-right: proper motion distribution. Bottom-left: parallax distribution. Bottom-right: 2D half-light radius for each population as a function of time. The dashed line indicates the mean radius and the coloured region indicates 1σ uncertainty. The dash dotted lines give the time of minimum radius for each population. The young population shows a minimum –0.40 Myr ago, while the old population shows a minimum –0.41 Myr ago.

In the text

Fig. 7. Properties of the runaway candidates. Left: 2MASS near-infrared colour-colour diagram for the runaway candidates. The early type runaways originally from NGC 6611 in the last 3 Myr are shown in blue. Late type stars are highlighted in red and are deemed unlikely to come from NGC 6611. We show the location of the zero-age main-sequence in grey (Pecaut & Mamajek 2013). The reddening line for a B0 V (top) and M0 V star (bottom) are shown with the dashed black lines for RV = 3.1 (Cardelli et al. 1989). Right: B, V and R composite mosaic image of NGC 6611 (ESO, press release 0926). We highlight the adopted search cluster search radius, with a radius of 0.1 deg (∼3 pc at a distance of 1.7 kpc), with the white circle. The trajectories of the O star runaways are shown in blue, while the B or later-type runaways are shown in orange. The ellipse shapes give the uncertainty on the location for the given trace-back time, which places the runaways near the centre of NGC 6611. The kinematic properties of these runaways are summarised in Table 2.

In the text

Fig. 8. AllWISE W4 (22 μm) grey-scale image of M16 and the surrounding region in Galactic coordinates. We show the current position of the runaways found outside M16 also shown in Fig. 7 adopting similar colours. We trace back the motion of these runaways relative to the cluster. UCAC2 27149134 and BD–14° 5040 have a clear bowshock which we show in Appendix A. We show for context the position of the O stars with the red crosses and indicate the position of NGC6604 and M17.

In the text

Fig. 9. Kinematic properties of the O stars in the vicinity of NGC 6611. Left: position and relative proper motion of the O stars near NGC 6611 overlapped on the DSS2 B, R and I colour image. The O stars shown in blue and cyan have velocities comparable or greater than the escape velocity and move away from the centre of NGC 6611. The O stars shown in magenta do not have this clear transverse motion away from the centre of NGC 6611, but could have had dynamical interactions in the past (see Appendix B. The O stars in red have no evidence for dynamical interactions in the past. The angular distance equivalent to 3 pc is given for a distance of 1.70 kpc. We highlight the centre of NGC 6611 with the dashed white circle (r = 0.05 deg). Right: relative (systemic) radial velocity ΔvR as a function of the absolute relative transverse velocity |ΔvT| for the O stars in NGC 6611, coloured similarly to the left panel. The O stars for which we have not used the Gaia astrometry (due to ruwe > 1.4) and thus have unknown |ΔvT| are shown on the left. The uncertainty on the radial velocity of NGC 6611 determined from the B0 to B4 stars is shown with the horizontal black shaded region. The median outward velocity determined by Kuhn et al. (2019) for NGC 6611 is shown with the vertical black shaded region.

In the text

	Fig. 10. Relative transverse velocity as a function of the determined range in kinematic age. We show the best-fit age and uncertainty determined with isochrone fitting with the black dashed line and shaded region respectively. The runaway O stars are shown in blue and cyan, while the runaway B or later-type stars are shown in orange. The O stars are coloured similarly as in Fig. 9.
In the text

	Fig. A.1. AllWISE (W1, W2 and W4 as blue, green and red, respectively) false-colour image of the bowshock produced by UCAC2 27149134 (18′ by 18′). The blue circle gives the position of the star. The red arrow gives the absolute transverse motion of the star, while the blue arrow gives the transverse motion relative to the cluster.
In the text

	Fig. A.2. AllWISE (W1, W2 and W4 as blue, green and red, respectively) false-colour image of the bowshock produced by BD–14° 5040 (18′ by 18′). The blue circle gives the position of the star. The red arrow gives the absolute transverse motion of the star, while the blue arrow gives the transverse motion relative to the cluster.
In the text

	Fig. A.3. Spitzer IRAC 8.0 μm image of W584 (2.8′ by 2.8′) The blue circle gives the position of the star. The red arrow gives the absolute transverse motion of the star, while the blue arrow gives the transverse motion relative to the cluster.
In the text

	Fig. A.4. Spitzer IRAC 8.0 μm image of HD 168137 (2.8′ by 2.8′) The blue circle gives the position of the star. The red arrow gives the absolute transverse motion of the star, while the blue arrow gives the transverse motion relative to the cluster.
In the text

Fig. B.1. B, V and R composite mosaic image of NGC 6611 (ESO, press release 0926). We overplot the near-infrared Hubble Space Telescope image of the centre of NGC 6611 (ESA/Hubble). The purple cross indicates the centre of NGC 6611, with three white dashed circles with radii of 0.1, 0.5 and 1.5 pc at a distance of 1.70 kpc. We show the motion of the O stars and their position at t = –130 kyr relative to NGC 6611, coloured similarly as in Figure 9, with the error ellipse giving the uncertainty on the proper motion.

In the text