© ESO 2020

1 Introduction

A complete and reliable understanding of the formation and evolution of massive stars has a direct impact on many fields of modern Astrophysics, thus serving an essential role in achieving an accurate interpretation the evolution of the Universe (see e.g., Woosley et al. 2002; Langer 2012; Abbott et al. 2016; Prantzos et al. 2018). In spite of this, important uncertainties and limits persist in our understanding of these objects, even in their early evolutionary phases. These uncertainties are related to our ability to determine accurate stellar parameters, particularly with regard to mass and effective temperatures, and to understand the implications of rotation, multiplicity, internal mixing, magnetic fields, and large scale motions for their structure and evolution. The mass and effective temperature scales of O-stars are still a matter of debate (Simón-Díaz et al. 2014; Herrero 2016; Holgado et al. 2018) and the actual role of rotation and binarity in the N and He enhancement of O stars are still controversial (Hunter et al. 2008; Rivero González et al. 2012; Martins et al. 2015). Uncertainties in mass-loss and binary interactions propagate to subsequent phases of massive star evolution, modifying them and leading to inaccuracies within the model predictions (Langer 2012). To improve this situation, large observational databases and the analysis of extensive samples of high quality OB-star spectra are required.

Galactic intense star-forming regions, due to their relatively close distance, allow us to perform a detailed observation of their massive stellar population. The Cygnus-X complex represents the most powerful star-forming region at less than 2 kpc from us (Rygl et al. 2012) that is conspicuous at all wavelengths and encompasses several rich very young OB associations. Hosting the largest number of near massive stars and an intense star-forming activity (Reipurth & Schneider 2008), this region is an ideal testbed for our theories about star formation, structure, and evolution, as well as the interplay between interstellar medium and stars and the dynamics and kinematics of OB associations and stellar groups.

Its richest association Cygnus OB2 represents the most obvious example of recent star formation. It harbors hundreds of OB stars that can be analyzed (Comerón & Pasquali 2012; Wright et al. 2015; Berlanas et al. 2018a) and although the optical extinction to Cygnus OB2 is high (A_V ~ 6 mag), it is not enough to prevent us from obtaining spectra of the bright massive members. It has been studied at all wavelengths with different spatial coverage (e.g., Massey & Thompson 1991; Hanson 2003; Herrero et al. 2002; Comerón et al. 2002; Negueruela et al. 2008; Wright et al. 2010, 2015; Sota et al. 2011; Comerón & Pasquali 2012; Rauw et al. 2015; Maíz Apellániz et al. 2016; Schneider et al. 2016; Morford et al. 2016; Berlanas et al. 2018a,b) as the actual extension of the association is not known with certainty (see, e.g., Comerón et al. 2008; Comerón & Pasquali 2012). Recent extensive studies of its massive population have been presented by Wright et al. (2015) and Berlanas et al. (2018a), who used a Hertzsprung-Russell Diagram (HRD) to discuss the evolutionary status of the association. However, their analysis was based on calibrations of stellar parameters versus spectral type, which introduces additional uncertainties, as has been shown by Simón-Díaz & Herrero (2014). Thus, an extensive spectroscopic study of the Cygnus OB2 early-type population is requisite for an accurate characterization of its massive stellar content. Moreover, since the Gaia satellite has provided accurate astrometry for the whole region (Arenou et al. 2018) and the distance to Cygnus OB2 has been re-estimated(Berlanas et al. 2019), we are poised to perform an accurate characterization of their massive stellar content, allowing us to create a much more precise HRD to interpret its evolutionary status.

This paper is organized as follows. In Sect. 2, we introduce our observing strategy. In Sect. 3, we describe the sample and selection criteria. In Sect. 4, we present the methods and tools used for the quantitative spectroscopic analysis of our sample and the results obtained. Then, in Sect. 5, we discuss and interpret all the results by creating the HRD and its spectroscopic version (sHRD). We also analyze the spatial and dynamical distribution of the sample from Gaia data. In Sect. 6, we summarize the conclusions of this work.

2 New observations and data reduction

We compiled new intermediate-high resolution and high signal-to-noise ratio (S/N) spectra to build up the most complete spectroscopic census of the currently known O-type stars in Cygnus OB2 (including the recently newly classified O-type members in the association by Berlanas et al. 2018a). These new data have been complemented with already available spectra from previous observations: our own observations and others drawn from our collaborators (see Sect. 3 for further details). We included in this sample all the stars identified as O-type stars within 1 deg radius centred on Galactic coordinates l = 79.8° and b = +0.8°. Since our goal is to perform accurate spectroscopic analysis for the whole sample, we need spectral data in the 4000–5000 Å wavelength range to access all the main diagnostic lines. Besides this, we also need to cover the Hα line in order to properly characterize the stellar wind.

To complete the database we performed four different campaigns (see Table 1 for specific dates) at the Isaac Newton Telescope (INT) in La Palma. We obtained intermediate-high resolution, high S/N (100 to 200) spectra for a total of 66 O-type stars of Cygnus OB2. Since we had already blue spectra for some of them from previous observations (see Sect. 3), we observed fewer stars in the blue than in the red, that is, 27 in the blue (3900–5100 Å) and 66 in the red (6100–6800 Å) wavelength regions. We used the Intermediate Dispersion Spectrograph (IDS) and the EEV10 detector, and chose the H1800V grating for the red wavelength region and the R1200B grating for the blue one. In the first case, the instrument configuration provides a resolving power of R ~10 000 at 6500 Å, which allows us to perform an accurate analysis of the Hα profile. Due to bad weather conditions during the 16A and 16B observing runs (dust in the air) the exposure times were increased by a factor of two with respect to the INT Exposure Time Calculator (ETC), resulting in times of around one hour to reach S∕N > 100 per pixel for a star of magnitude B = 14.0 in the red. In the second case, the instrument configuration provides a resolving power of R ~5000 at 4500 Å, which is adequate for a proper determination of the spectroscopic stellar parameters. In the blue wavelength region, the dust has a higher impact on the exposure times and, consequently, the ETC times were increased by a factor of three. Therefore, in the blue, we needed ~7000 s to reach S∕N > 100 per pixel for a star of magnitude B = 14.0. In all cases, exposure times above 30 min were divided into ranges of 20 min to reduce the effect of cosmic rays.

The observational plan was divided into two phases, extending over four semesters, where two spectra were obtained foreach star on different days (belonging to the same semester or not) to detect possible spectral or radial velocity variations. The four observing runs were held throughout the summer of 2016 and 2017, for a total of 24 nights distributed as indicated in Table 1. The spectra were reduced using IRAF¹ with standard routines for bias, flat-field subtraction and wavelength calibration. We also corrected the spectra for the radial velocity using the Doppler shift of metal and He I spectral lines.

3 The sample

We built up the most complete spectral sample for the already known O-type population in Cygnus OB2 by complementing new observationswith already available spectra (own and from collaborators) whose observational details are described as follows.

Most of the available optical–blue data (see Table 2) comes from the Galactic O-Star Spectroscopic Survey (GOSSS, Sota et al. 2011, 2014; Maíz Apellániz et al. 2016), which obtained a large number of Galactic O-type high S/N (~300) blue-violet spectra at R ~ 2500–3000. The GOSSS data used in this work were observed between 2007 and 2014 using three different instruments: the Albireo spectrograph at the 1.5 m telescope of the Observatorio de Sierra Nevada (OSN), the TWIN spectrograph at the 3.5 m telescope of Calar Alto Observatory (CAHA, Centro Astronomico Hispano-Aleman), and the ISIS spectrograph (blue arm) at the William Herschel Telescope (WHT) of the Observatorio del Roque de los Muchachos (ORM).

We also have access to spectra from The Hobby-Eberly Telescope (HET) obtained by A. Pellerin. These data were observed throughout 2012 in four different spectral ranges (3800–4700, 4700–5700, 5300–6300 and 6400–7400 Å) using the High-Resolution Spectrograph (HRS) at R ~30 000 and the 600g4739/6302 gratings. We have a total of 48 spectra for 12 stars (one in each spectral range). The S/N of these spectra is variable, ranging from 100 to 300.

Finally, we also include spectra observed by us in different runs at the WHT using the ISIS spectrograph (a total of 30 spectra, 15 in each blue and red arms). These data were acquired in two different observing runs. For 2015, we chose the H2400B grating for the blue arm, which provides a resolving power of ~13 600 at 4500 Å. As the dichroic allows us to observe simultaneously with the red arm, we used the R1200R grating (R ~ 9300 at 6500 Å). Additionally we included stars observed by A. Herrero in 2003 also using ISIS at the WHT with the same configuration (H2400B and R1200R gratings). All ISIS spectra have S∕N >100 (from 100 to 200, depending on the stellar magnitude).

A brief summary of the data used in this work is shown in Table 2. The final spectral sample is composed of 78 O-type stars (12 of them previously classified as SB2) whose spatial distribution is shown in Fig. 1. The main data and photometric information are presented in Table A.1.

3.1 Binary fraction

A large number of multiplicity studies among massive stars have been conducted during the last years (see Sana 2017, for a brief overview) where the number of detected massive O-type binaries has been continuously revised and updated. The Cygnus OB2 Radial Velocity Survey is one of the most relevant multiplicity census carried out in this association (Kiminki et al. 2007, 2008, 2009, 2012, 2015; Kobulnicky et al. 2012), where the binary properties of its massive star population have been statistically studied. Kobulnicky et al. (2014) list 48 massive OB multiple systems known in Cygnus OB2, 26 of them containing at least one O-type star. In recent years, Maíz Apellániz et al. (2016, 2019) have provided spectral classification for most of the detected companions.

In this work, we found two new SB2 systems within our sample of O-type members of Cygnus OB2: J20301838 + 4053466 and CPR2002 #A32² (see Fig. 2 for an example of the binary detection). A list of the currently known binary (or multiple) systems in the association is presented in Table 3 where we also show the binary type (single-lined SB1 or double-linedSB2), the components, spectral types, and bibliographic references. Taking into account the results of this work, 29 of 78 O-type stars in the region are part of binary/multiple systems (representing a fraction of 37%). Sana & Evans (2011) found that at least 45–55% of the O star population in clusters and OB associations is comprised of spectroscopic binaries, making it highly likely that more binaries have gone undetected in our sample.

We highlight some individual remarks for three of the most relevant multiple systems known in Cygnus OB2: Cyg OB2 #5, Cyg OB2 #22, and Cyg OB2 #8. The first is a multiple system for which Maíz Apellániz (2010) lists up to four components. Here, #5A is the star that gives its name to the system, and is an eclipsing, contact binary formed by an O6.5 – O7 and a WN9/Ofp stars (Rauw et al. 1999) that is overluminous for its mass. The #5B component, at 0.9″ from #5A, has been classified as O7 Ib(f)p var? by Maíz Apellániz et al. (2019). The GOSSS catalogue gives separated spectra for #5A and #5B. For the trapezium-like system Cyg OB2 #22 Maíz Apellániz (2010) has identified a total of ten components and the GOSSS Survey has provided spectral data for the first five O components (A, B, C, D and E). At least four O-type visual components have been detected in the Cyg OB2 #8 trapezium-like system (Maíz Apellániz et al. 2012). The first component, #8A, has been identified as a SB2 star composed by two early O-type stars (De Becker et al. 2004; Maíz Apellániz et al. 2019). Both Cyg OB2 #8 and #22 have been traditionally treated as multiple stars but as already pointed out by Bica et al. (2003), they should be considered as two clusters likely constituting the core of Cygnus OB2 and hosting the brightest stars of the area. A recent analysis by Maíz Apellániz et al. (2020) has confirmed their cluster nature, and showed that both are at a distance of 1.7 kpc (the same distance as obtained by Berlanas et al. 2019, see Sect. 4.5) and have similar proper motions. Therefore, they are a double cluster separated by 2.7 pc in the plane of the sky and the paper above proposes naming them as Villafranca O-007 and Villafranca O-008, respectively. As a final remark, we note that although A36 was classified as B0 Ib + B0 III by Kobulnicky et al. (2014), it has been revised by Kiminki et al. (2015) who indicated that both components may be slightly hotter. They classified it as O9 – O9.5 III: + O9.5 – B0 IV: thus, we include the system in our sample. Recently, it was reclassified as O9.5 II + O8.5 III by Maíz Apellaníz (in prep.) and we adopt this classification in our work.

Fig. 1 Inverse Spitzer 8 μm image showing the location of the known O-type population in Cygnus OB2 (78 stars represented by green dots). The solid line circle delimits the 1 degree radius area adopted for the association. For reference, the dashed-dotted line circle shows the area considered by Wright et al. (2015) indicating the core of the association.

Fig. 2 Example of the new SB2 detections: HeI lines in stars J20301838 + 4053466 (top) and A32 (bottom).

3.2 Selecting the final spectral sample

The detected SB2 stars (see Table. 3) are excluded from our final sample for the spectroscopic analysis. We also excluded the Cyg OB2 #5B and Cyg OB2 #22E components since they are not suitable targets for an accurate determination of the stellar parameters. The #5B component shows a peculiar variable spectra with many emission lines and the faint #22E component shows a noisy spectrum. Tentatively, we also exclude the Cyg OB2 #22D component. For this star, we have spatially resolved spectroscopy from sources 1 and 6 of Table 2, but a detailed inspection of the spectra reveals peculiar line profiles in the red spectrum that are not compatible with those in the blue. Since these spectra were obtained in different observing campaigns with a lapse of a few years, new data are required to evaluate stellar variability or even a possible target confusion. Finally, the star J20272428 + 4114458 was classified as O9.5 V by Comerón & Pasquali (2012) but it has been reclassified as B0 IV in this work. As it is the only star that changes from O to B-types, we decided to keep it in the sample to see its behavior and inspect its properties.

We therefore have a final sample of 62 O-type stars (plus the star indicated above, out of a total of 78 present in our census) suitable for a quantitative spectroscopic analysis. Since we have compiled spectral data from different sources, in many cases, we have different spectra for the same star. In these cases we adopted the following strategy. To determine the distribution of rotational velocities (henceforth, sample A), we select the diagnostic lines in each spectrum following the criteria given in Sect. 4.1. Thus, we chose the spectra with the higher resolving power since it allows us to better identify spectroscopic binaries and to disentangle the macroturbulence broadening component (v mac), obtaining more accurate projected rotational velocity (v sin i) determinations (see Simón-Díaz & Herrero 2007, 2014; Simón-Díaz et al. 2017, and references therein). However, in order to determine the main stellar parameters, we use the spectra with the best quality in terms of S/N (henceforth, sample B) and, whenever possible, the best resolution for the H and He lines.

As an example, for the Cyg OB2 #6 star, we have data from three different instruments: sources 2, 4, and 6 of Table 2. In order to determine the rotational velocity andmacroturbulence broadening component, we use the source number 4 because of the higher spectral resolution (R ~30 000) and the fact that the O III 5592 line is included in the spectrum. On the other hand, for the stellar parameters we use sources number 2 (for the optical-blue range) and 6b (for the optical-red range) because of the higher quality of the spectra (in terms of S/N). We note that we are assuming no significant variations between individual blue and red spectral data (e.g., no strong variations in Hα); however in the cases where only one spectrum is available for each wavelength range, this assumption cannot be confirmed. A final list of the whole O-type population of our sample, along with the new spectral classification and the spectral data selected for each case (A or B) is shown in Table A.2.

Fig. 3 Comparison between the stellar parameters obtained using spectra from the IACOB-OWN (R ≥ 25 000) and the GOSSS surveys (R ~ 2500) for a sample of 37 standard O-type stars in common. The same v sin i values and set of diagnostic lines were used for the spectroscopic analysis. Red arrows indicate upper limits.

3.3 Implications of using an inhomogeneous spectral sample

Since our final sample is spectrally inhomogeneous in terms of R, S/N, and wavelength coverage, we evaluated the implications of using data with different spectral characteristics for the spectroscopic analysis.

In the case of the wavelength coverage, we know that a limited set of lines affects the parameter determination (mainly for gravity and, therefore, mass) and the errors decrease as we increase the number of diagnostic lines (Berlanas et al. 2017). Nevertheless, all our spectra cover the whole range where the main diagnostic lines for O-type stars are located (4000–5000 Å) and we have compensated for the lack of the Hα and He I+II 6678 lines that existed for some stars by carrying out new spectroscopic observations (see Sect. 2).

Although all the data used for the analysis have high S/N (exceeding 100 in all cases), the resolving power varies in a wide range (from 2500 to 30 000, see Table 2). We therefore decided to evaluate the possible effect that resolution could have in the derivation of the stellar parameters using spectra from the IACOB-OWN³ (R ≥ 25 000) and GOSSS (R ~ 2500) surveys. For the test, we selected the same 37 standard O-type stars (at least one for each type and luminosity class) from both surveys. In order to avoid possible additional uncertainties we performed the spectroscopic analysis considering only the common set of lines and the same values of v sin i and v mac. In Fig. 3, we present a comparison between the derived effective temperature (T_eff) and surface gravity (log g). We can see that both parameters are only slightly affected by the spectral resolution and most of the stars are located within the usual uncertainties without a global impact. Mean and standard deviation of the difference in T_eff and log g from high resolution IACOB-OWN and low resolution GOSSS spectra is 0.4 ± 0.9 kK and −0.05 ± 0.08 dex, respectively. We note that the pair (v sin i, v mac) used in this test was obtained from the high resolution IACOB-OWN spectra since the O III 5592 line is included and it is one of the best lines for an accurate line-broadening characterization. We also know that resolution affects the v sin i determination: Simón-Díaz & Herrero (2014, henceforth, SDH14) show that a low resolution could lead to an overestimate of the derived v sin i of ~25 (±20) km s⁻¹ in stars rotating more slowly than 120 km s⁻¹. Sabín-Sanjulián (2014) shows that this could lead to errors in the gravity determination of up to 0.1 dex and a consequent change in T_eff up to 1000 K (although these are upper limits and changes would tend to be smaller in most cases).

4 Methods and results

4.1 Line-broadening characterization

Since metallic lines do not suffer from strong Stark broadening nor from nebular contamination, they are best suited for obtaining accurate v sin i values. Unfortunately, they are not present in all our spectra. Consequently, we assumed the following criteria based on the works by Herrero et al. (1992), Simón-Díaz et al. (2014) and Ramírez-Agudelo et al. (2013): (i) we first choose the O III 5592 or Si III 4552 diagnostic lines, when either of them is available. If both are present, we prioritize the O III 5592 line since it is usually stronger and more isolated than Si III 4552 for O-type stars; (ii) if none of them are present, we use nebular free or weakly contaminated He I lines (an unweighted average of the He I 4387, He I 4471, He I 4713 and He I 4922 lines); (iii) if any star suffers strong nebular contamination or the He I lines are weak or too noisy, we then rely onHe II 4542; (iv) finally, if the stellar wind is weak and the line is of photospheric origin, we also use N V 4603.

As in Ramírez-Agudelo et al. (2013), for those cases in which different diagnostic lines are available we compared the measured v sin i values. We found similar results as the aforementioned authors, obtaining a good degree of agreement. In all cases differences between lines are within the limits represented in Fig. 4 and have no impact on global or individual conclusions. A summary of the diagnostic lines used for the line-broadening characterization of our final sample of 63 O-type stars is shown in Table 4.

We used the iacob-broad tool (Simón-Díaz & Herrero 2007, 2014) to obtain the projected rotational velocities and the macroturbulent broadening component for the whole stellar sample (excluding detected SB2 and problematic stars). This is a user-friendly IDL procedure for the line-broadening characterization of OB stars, based on a combined Fourier transform (FT) plus a goodness-of-fit (GOF) methodology. It allows the user to determine easily the stellar projected rotational velocity and the amount of extra broadening (assuming a radial-tangential profile) from a specifically selected diagnostic line. The FT technique is based on the identification of the first zero in the Fourier transform of a given line profile (Gray 2008; Simón-Díaz & Herrero 2007). The GOF technique is based on a comparison between the observed and a synthetic line profile that is convolved with different values of v sin i and v mac to obtain the best-fit by means of a χ² optimization. The main advantage of the iacob-broad analysis is that we obtain two independent measurements of the v sin i (resulting from either the FT or the GOF analysis) whose comparison is used as a consistency check and to better understand problematic cases. The results, with uncertainties in the range 10–20%, are given in Table A.3.

Fig. 4 Comparison of projected rotational velocities resulting from the Fourier Transform (FT) and the Goodness-Of-Fit (GOF1) techniques obtained using the iacob-broad tool. Dashed lines represent a difference of 20 km s−1 or 20% from the 1:1 relation, whichever is the largest. Red squares indicate those stars out of these limits. Different colors indicate the diagnostic lines used for the line-broadening characterization.

4.2 Distribution of rotational velocities

Since there are no previous studies focused on the distribution of rotational velocities for O-type stars in Cygnus OB2, we explore the information contained in our data, although we have to pay attention to their inhomogeneity in resolving power when comparing results from different spectra.

We first compared the v sin i values derived from the Fourier transform technique (FT) and those derived from the goodness-of-fit technique (GOF) as an assessment of the reliability of the results from the line-broadening analysis (see Fig. 4). We see that only a few stars have discrepancies beyond the adopted uncertainties (we adopt limits of 20 km s⁻¹ or 20%, whichever is the largest since SH14 found that the agreement is always better than 20% but for low rotational velocities effect like the spectral resolution or the microturbulence may play an additional role) with no clear systematic trends due to the lines used for the analysis. For the stars for which we find good agreement (below the limit given above), we assumed v sin i and v mac provided by the GOFtechnique since it is less affected by the subjectivity in the selection of the first zero of the FT (see Simón-Díaz & Herrero 2014). For those few cases where we find larger discrepancies (red squares), we identified clear signatures of asymmetry in the line, perhaps due to a broad-line secondary component. In these cases, the GOF solution tries to fit the line by increasing v mac, hence resulting in a lower v sin i. We then used the v sin i provided bythe FT technique (which is not very much affected) and the associated v mac (obtained from the fit of the line assuming a fixed v sin i value derived fromthe FT). However, we remark that nearly all discrepancies are close to the adopted limit and, therefore, our choice has little impact on the final distribution.

To derive and correctly interpret the distribution of rotational velocities of our sample, we must take into account that our spectral sample for Cygnus OB2 is inhomogeneous in terms of resolving power (see Table 2), so we must consider the lower limiting v sin i of each of the spectral resolutions used in this work. Using a rough approximation (v_lim ~ c∕R), values of ~100, ~60, ~30, and ~10 km s⁻¹ are the lower limits that can be measured at R ~ 2500, R ~ 5000, R ~ 10 000, and R ~ 30 000, respectively. We thus show the distribution of rotational velocities for each group of stars observed at a given resolution separately (see Fig. 5, left). We see that the main peak of slow rotators (v sin i ≤ 100 km s⁻¹) is composed mainly of stars observed at R ~ 5000, with an associated limiting velocity of ~60 km s⁻¹. Therefore, we are not able to detect lower velocities and the main peak of slow rotators is shifted to higher values. The second peak of stars rotating at intermediate velocities (v sin i ~180–200 km s⁻¹) is mainly composed by stars observed at R ~ 2500, with a corresponding limiting velocity of ~100 km s⁻¹. Such a low resolution may result in some cases in v sin i values overestimated by ~20 km s⁻¹ (i.e., one bin). Thus the peaks in our distribution are probably shifted towards higher v sin i values (but the results from higher resolution spectra will not be affected). Finally, we remark two more points in the distribution: (a) the very slow rotators (v sin i ≤ 60 km s⁻¹) bins where we see the lack of extreme slow rotators (v sin i ≤ 20 km s⁻¹), indicating that some other physical process is preventing the derivation of lower values (probably microturbulence, see Simón-Díaz & Herrero 2014); and (b) the short tail of fast rotators that extends only up to v sin i ≤ 280 km s⁻¹.

To check how representative the distribution of rotational velocities that we obtained really is and to explore possible future research directions, we compared our results to other similar works, such as the study of Simón-Díaz & Herrero (2014) for a general sample of O-type stars in the Milky Way (MW) and the study done by Ramírez-Agudelo et al. (2013) in 30 Doradus in the Large Magellanic Cloud (LMC). We applied the Kolmogorov-Smirnov two-sided test (see Smirnov 1939; Hodges 1958) to check if the difference between these distributions and ours is statistically significant. The p-values returned (10⁻⁵ and 10⁻⁹ for the MW and 30 Doradus samples, respectively) reject the null hypothesis that the distributions are the same. The same test carried out for the distributions up to 280 km s⁻¹ returns values of p = 0.911 and p = 0.904, indicating that the lack of very fast rotators in the Cyg OB2 distribution is significant.

In Fig. 5 (middle, right), we compare our results with those obtained by these aforementioned studies. In the middle panel, we compare our results in Cygnus OB2 (in yellow) and the distribution obtained for a general sample of Galactic O stars in the MW (in blue). Both are Galactic samples, the first one composed by O stars belonging to the same young star-forming region and the latter composed by bright O stars belonging to many regions within the Galaxy. For the latter, the sample consists of 116 O-type stars⁴ from the IACOB database that have been observed at R = 46 000 and 25 000. Both studies have followed the same methodology for the line-broadening characterization although the wavelength coverage of the IACOB spectra allows them to use the Si III 4552 and O III 5592 lines for all their sources. The sample of Galactic O stars presents a clear bimodal distribution, with a main peak of slow rotators in the 40–80 km s⁻¹ bin (instead of the 80–120 km s⁻¹ found for Cygnus OB2) and a tail of fast rotators reaching 450 km s⁻¹ (while in Cygnus OB2 we don’t reach 300 km s⁻¹). Taking into account the lower v sin i limit that can be detected for a given resolution, the difference in the low peak can be understood as due to the impact of the lower resolution in the Cygnus OB2 sample (note also the height of the 0–40 km s⁻¹ bin in the MW sample). On the contrary, the lack of the fast rotator tail in the Cygnus OB2 sample is much more difficult to understand and demands a more detailed study.

In the right panel, we compare the v sin i distribution resulting from our study in Cygnus OB2 (in yellow) and the distribution for a sample of 216 O-type stars in 30 Doradus (in green, see Ramírez-Agudelo et al. 2013). Both are samples of stars belonging to young star-forming regions that differ in the number of targets (66 versus 216), the environment metallicity (Z_LMC = 0.5 Z_⊙) and the resolving power of the spectral data (R ~8000 for all targets in 30 Dor). Again, the same methodology was used for the line-broadening characterization, and due to similar wavelength restrictions both had to rely in many cases on He I instead of metallic lines. The peak of slow rotators in 30 Dor is found at 40–80 km s⁻¹ and there is also a more evident tail of fast rotators reaching 600 km s⁻¹, something expected for O stars located in lower metallicity regions as they lose less angular momentum through stellar winds (see, e.g., Langer 2012). The 30 Dor distribution is therefore very similar to the MW one.

In both cases we find a similar behavior, showing significant differences with respect the distribution derived for Cygnus OB2. The main peak of slow rotators is shifted to lower velocities, which can be easily explained by the different resolutions. In addition to this, the fact of having limited number of epochs for our spectra affects the detection of SB2 stars, overestimating the derived v sin i values and increasing the peak. However, the extended tail of fast rotators in the MW and 30 Dor distributions is not present in our Cygnus OB2 data, which demands further attention.

Fig. 5 Distribution of rotational velocities of our final sample of 63 O-type isolated stars in Cygnus OB2. Left: histogram color-coded on the basis of the resolving power of the spectra used for each target. Middle: v sin i histogram showing the comparison between the distribution of rotational velocities found for our sample of O-type stars in Cygnus OB2 and the distribution of O-type stars found in the MW. Vertical black lines indicate the associated Poisson noise. Right: as in middle panel but compared to 30 Dor.

4.3 Spectral classification

Since we have increased the quality of the previous available spectral data for a large number of stars in our sample, we revised the spectral classification of these stars using the following criteria.

The main classification criterion used for O-type stars is the comparison of the He II 4542 and He I 4471 lines, whose ratio is unity for an O7 type star. He I tends to increase in strength with decreasing temperature while He II decreases. Spectral types earlier than O8 can be classified using these criteria (see e.g., Walborn in Gray & Corbally 2009). For later O types the relative strengths of He II 4542/He I 4387 and He II 4200/He I 4144 may be used, representing the main diagnostic for types O8–B0 (Sota et al. 2011). We adopted the criteria of Sota et al. (2011) using the list of qualifiers for O spectral types summarized in their work. Regarding the luminosity class, the classification criteria for early O-type stars were introduced by Walborn (1971, 1973), taking into account the emission effects in the He II 4686 and N III 4634-4640-4642 lines. Sota et al. (2011) provided an updated classification system for late O-type stars, which we adopted in this work. Although the classification for our O-type stars was based on He lines as described above, we corroborated their spectral types by using the Marxist Ghost Buster tool (MGB, Maíz Apellániz et al. 2012). It compares the observed spectra with a standard library of O stars (in this work the GOSSS library, see Maíz Apellániz et al. 2016). This interactive software allows us to vary the spectral subtype, luminosity class, line broadening, and spectral resolving power of the standard spectrum until we obtain the best match. In addition, it also allows us to combine two standard spectra (with different velocities and flux fractions) to fit SB2 binaries.

The previous (from the literature) and new (from this work) spectral types are listed in Table A.2. There is quite good agreement between the previous spectral types and the new ones. In nearly all cases, the differences are not larger than half a spectral subtype and one luminosity class. In three cases, the differences reach one spectral subtype (J20291617+4057371, MT91-455 and A20). Only in one case we find a significant change in luminosity class (B18, reclassified from O7: Ib to O7.5 IV:(f)), indicating the interest of a follow-up for this target.

4.4 Determination of the main spectroscopic parameters

We determined the fundamental stellar parameters for the final sample of 63 O-type stars in Cygnus OB2 performing a quantitativespectroscopic analysis based on synthetic FASTWIND models (Santolaya-Rey et al. 1997; Puls et al. 2005; Rivero González et al. 2012) and the iacob-gbat tool (Simón-Díaz et al. 2011; Holgado et al. 2018).

The code allows us to accurately determine the basic six spectroscopic stellar parameters for OB type stars: the effective temperature (T_eff), surface gravity (log g), wind-strength parameter (Q, defined as Ṁ /(v_∞ R)^1.5), helium abundance (Y (He), defined as N(He)/N(H)), microturbulence (ξ), and the exponent of the wind velocity-law (β) from their Hand He lines⁵. Once the observed spectrum is processed, the tool compares the observed and the synthetic line profiles (from FASTWIND models, in our case) by applying a χ² algorithm. It computes the line-by-line χ² distributions, estimating the goodness-of-fit for each model within a subgrid of models selected from the global grid. Then it iteratively computes the global χ² distribution, from which the final parameter values and their associated uncertainties are estimated. The given parameters are the mean values computed from the models located within the 1-σ confidence level of the total χ² distributions (after each model has been weighed by its corresponding χ² value). Then their uncertainties are given by the standard deviation within the 1-σ level. Our grid of models covers the wide range of stellar and wind parameters considered for standard OB-type stars, from early-O to early-B types and from dwarf to supergiant luminosity classes (see Table 5 for grid details). As explained in Sect. 3.2, when different spectral data were available for the same target, we chose the spectra with the best quality in terms of S/N (and highest resolving power if possible). Derived stellar parameters are shown in Table A.3. We complement these results with a series of figures (see Appendix A) in which the best-fitting model resulting from the analysis of each star is overplotted on the observed spectrum.

4.5 Radii, luminosities, and masses from Gaia DR2

The iacob-gbat tool also allows us to obtain directly physical parameters when the absolute magnitude M_V is provided. It computes the radius (R), luminosity (L), and spectroscopic mass (M_sp) following the same methodology as above. With this aim, the tool uses the equation (Eq. (1)) introduced by Kudritzki (1980), where the $\mathcal{V}$ parameter represents the integral in wavelength of the emergent flux of the model (Matthews & Sandage 1963): $$5 \text{log}(R/R_{\odot })=29.57-M_V+\mathcal{V}$$(1)

and then we use the usual formulae to derive masses and luminosities: $$\text{log}(L/L_{\odot })=2 (R/R_{\odot })+4 \text{log}(T_{\text{eff}}/T_{\text{eff},\odot }),$$(2) $$\text{log}(M_{\text{sp}}/M_{\text{sp},\odot })=2 (R/R_{\odot })+\text{log}(g/g_{\odot }).$$(3)

By using the unprecedented parallaxes provided by Gaia in its second data release (DR2, Gaia Collaboration 2018), we derived M_V for all our stellar sample with reliable astrometry (RUWE ≤1.4, see Lindegren et al. 2018, for details on this criterion). However, we should take into account that these parallaxes could be affected by systematic uncertainties of up to 0.1 mas that are not well understood (Luri et al. 2018) and, therefore, should be taken with caution.

Since V -band photometry is not available for all stars of our sample (see Table A.1) we followed the same procedure used by Comerón et al. (2008) to derive M_V values. With M_V being the absolute magnitude we have $$M_V=5-5*\text{log}(d)+V-A_V=-\text{DM}+V-A_V,$$(4)

where DM is the distance modulus. The visual extinction A_V can be defined in terms of the (V − K_s) color excess as $A_V=E(V-K_s)+A_{K_s}$. Then, M_V can be defined as $$M_V=K_s-\text{DM}-A_{K_s}+{(V-K_s)}_o.$$(5)

We used distances provided by Bailer-Jones et al. (2018, henceforth, BJ18) to estimate individual values of the DM and adopted optical and near-IR photometry from the USNO-B and 2MASS catalogs and unreddened intrinsic colors from Martins & Plez (2006). The use of the alternative prior of Maíz-Apellániz (2001), which is optimized for early-type stars in the Galactic disk, provides similar distances within the uncertainties. We used the R_V = 3.1 extinction law from Rieke & Lebofsky (1985) to derive $A_{K_s}=0.092*E(B-K_s)$. We note that our results will be dependent both on the extinction law adopted and on the assumption of a constant average R_V. This coefficient depends on the properties of the absorbing dust grains and thus on the kind of region containing them (see, e.g., Maíz Apellániz & Barbá 2018). Nevertheless, for the whole sample we find good agreement between the derived M_V values and those obtained from the Martins & Plez (2006) calibration. The M_V uncertainties were obtained from the distance errors provided by BJ18. Uncertainties related to the apparent magnitudes are negligible compared with those obtained from distances, that are of the order of ± 0.13 magnitudes.

The final list of targets with reliable astrometry for the O-type population in Cygnus OB2 is composed of 52 stars, whose corresponding Gaia DR2 sources, estimated distances by BJ18 and derived M_V values are presented in Table A.4. This table also includes the membership group of each star from Berlanas et al. (2019, henceforth B19). Group 2 is the main Cygnus OB2 group at ~1.76 kpc and Group 1 is a closer one at ~1.35 kpc⁶. Group 3 contains foreground or background contaminants. Group 0 indicates that the object lies between Groups 1 and 2 and could not be reliably assigned to any group. Luminosities, radii and spectroscopic masses from the iacob-gbat tool are given in Table A.4 as well. Uncertainties include the iacob-gbat formal uncertainties for the stellar parameters and those related to M_V.

5 Discussion

5.1 Evolutionary status

Using the luminosities derived from M_V and the iacob-gbat tool, we placed the sample of stars that have passed the criteria for reliable astrometry (RUWE ≤1.4) in the Hertzsprung-Russell diagram (HRD). The mean uncertainties are ~0.01 dex in log T_eff and ~0.05 dex in log(L∕L_⊙). We decided to explore four different stellar models (two families with two initial rotational velocities) in order to assess possible uncertainties in the obtained distribution. We used rotating and non-rotating Geneva (Ekström et al. 2012) and Bonn (Brott et al. 2011)stellar evolutionary tracks and isochrones. In addition, we also divided the sample into the different stellar groups found by B19.

In Fig. 6 (top left hand panel), we show the HRD using Geneva non-rotating evolutionary models. The age of the O population is well constrained between 1 and 6 Myr, as estimated by other authors (Hanson 2003; Negueruela et al. 2008; Comerón & Pasquali 2012; Wright et al. 2015; Berlanas et al. 2018a). We find the most massive star at around 75 M_⊙ (Cyg OB2 #22A) and the least massive one at 18 M_⊙ (Cyg OB2 #23), as expected for O-type stars. However, different stellar age groups can be clearly distinguished within the whole population⁷. We find the older group following the ~6 Myr isochrone with initial masses between 18 and 40 M_⊙. The group is well separated in luminosity class, with the dwarfs located along the Main Sequence (MS) and the two supergiants A29 and J20283038 +4105290 close to (but before) the terminal-age main sequence (TAMS). A more numerous and younger group is found along the ~3 Myr isochrone with initial masses between 20 and 60 M_⊙. It is more tightly concentrated than the others, containing most of the population. Interestingly, all the stars belonging to Group 1 (the foreground group at ~1350 pc) and Group 0 (located between the foreground and the main group) identified by B19 seem to follow the same isochrone, except for one object in each group. A third smaller group containing the apparently youngest stars of the sample seems to follow the ~1.5 Myr isochrone with initial masses between 20 and 75 M_⊙. The hottest stars of our sample, Cyg OB2 #22A and Cyg OB2 #7 are includedin this group (both O3 supergiants, see Walborn et al. 2002). Without the presence of the most massive stars (with M >60 M_⊙) it is possible to assign the rest of the stars to the group at ~3 Myr. In this case, there would only be two star formation bursts (at ~6 and ~3 Myr) and the position of the luminous stars close to the ZAMS would be due to peculiar evolution, like merger products after binary interaction (as it was suggested e.g., in Negueruela et al. 2008), effects of the accretion rate (see Holgado et al. 2020), or chemically homogeneous evolution (see below). We make note of the lack of stars close to the ZAMS between 40–60 M_⊙ in this group, producing an apparent gap close to the ZAMS in the global diagram. This effect has been pointed out by other authors (e.g., Herrero et al. 2007; Castro et al. 2014; Sabín-Sanjulián et al. 2017) and was recently revisited by Holgado et al. (2020).

The same HRD is shown in the bottom left-hand panel of Fig. 6 but this time with a considering rotation (v∕v_crit = 0.4). As stated in previous works (e.g., Maeder & Meynet 2000; Wright et al. 2015; Berlanas et al. 2018a), Geneva rotating stellar models provide older main sequence stellar ages than their non-rotating models. The most massive members are more affected, and in average we find all the age groups ~2 Myr older than in the non-rotating case. Interestingly, the most massive stars could also belong to an older population with stars of lower mass located further away from the ZAMS as the rotating isochrones of 3–5 Myr move bluewards for the higher masses. The Bonn stellar models (right hand panels of Fig. 6) do not exhibit large differences in stellar ages when using non-rotating or rotating models with initial rotational velocities (v_ini) as high as 330 km s⁻¹. This different behavior is most likely to be ascribed to their different treatment of angular momentum transport and magnetic fields in the stellar interior. The ages derived with the Bonn models are similar to those obtained with the non-rotating Geneva models. The most relevant difference between both families of models is that Bonn models provide a more extended TAMS, and when the models include rotation, the ZAMS is slightly displaced to higher luminosities and lower temperatures. In consequence, we find more stars close to (or on) the ZAMS.

In summary, our data indicate the presence of at least two main episodes of star formation in the main group of the Cygnus OB2 association (the one at ~1.76 kpc in B19) at ages, according to Geneva non-rotating tracks, of ~6 and ~3 Myr. A third episode at ~1.5 Myr seems possible, depending on the interpretation of the most massive stars. These results are fully in agreement with the analysis by Negueruela et al. (2008). They assigned an age of 2.5 Myr to the association, with evidence of a slightly older population (much more evident in our results), and they explore the possibility that the most massive stars belong to a younger population, but ultimately conclude that this is not the most probable case. Their findings are also consistent with Wright et al. (2015), who propose a nearly continuous star formation between 1 and 7 Myr (but our results favour a more episodic history). Although the region has currently an intense stellar formation activity (seee.g., Schneider et al. 2006), there is no evidence of star formation episodes more recent than 1–2 Myr.

That situation resembles the one established by Schneider et al. (2018a,b) in 30 Dor, with a main peak of star formation 8 Myr ago, a decline in the last Myr and intense star formation in between. The ages derived for the stars in Cygnus OB2 become older by 1–2 Myr when considering Geneva models with rotation, except for the youngest population, which approaches the theoretical ZAMS when considering rotating Bonn models. For the smaller group at ~1.35 kpc and the few stars that according to B19 seem to lie between the two main groups, we find that star formation took place ~3 Myr ago, without substantial activity at other ages (but we note that the sample is small). Two stars might belong to an older population at again ~6 Myr, but the small sample size hinders us from generalizing the result and claiming that this is another episode of massive star formation. Careful research of the associated population at lower masses would be required to confirm or reject such a possibility.

We also have to consider the possibility that some of the stars were born with a different initial rotational velocity than the rest or have followed peculiar (e.g., not single) evolutionary paths. For example, the two more massive stars (Cyg OB2 #7 and Cyg OB2 #22A) are a bit puzzling. At face, we would assign them to a 1.5 Myr group, but then they leave a gap between 25–40 M_⊙ (depending on the diagram we use) and 60 M_⊙ close to the ZAMS. Alternatively, we could assume that they belong to an older population and are the result of alternative evolutionary channels. Cyg OB2 #7 and Cyg OB2 #22A could be theproduct of fast rotation evolution, while other stars may have followed an essentially non-rotating track. In Fig. 7 we plot the stars on the HRD with Geneva rotating tracks and isochrones. We see that it would be possible to assign Cyg OB2 #7 and Cyg OB2 #22A to the 3–4 Myr population, with no mass gap in the population. In this scenario we would expect these high luminosity stars to show a present-day rotational velocity of the order of the observed v sin i, but also to have a high surface Helium abundance (higher than the enhanced Helium abundance found for Cyg OB2 #7), which is not supported by our analysis. However, the distribution of fast rotators (v sin i > 200 km s⁻¹) on the HRD shown in Fig. 7 concentrates at the lowest masses and it does not offer any support for this scenario. On the other hand, the stars could be interpreted as mergers after binary interaction (de Mink et al. 2014; Wang et al. 2020), again corresponding to a 3–4 Myr age. In that case, helium abundances would be closer to the observed ones as a consequence of the fresh hydrogen supply, but the expected rotational velocities would be very large (although the observed one would be affected by the inclination of the rotational axis). Although it does neither contradict it, if the most massive stars belonged to a population of fast rotators following an isochrone between 2 and 4 Myr, we would expect them to be rotating fast and closely following a chemically homogeneous evolution. Yet, both stars have modest v sin i and only the He abundance of Cyg OB2 #7 is is slightly above the solar value.

Finally, we evaluated the possible impact of the lack of detected binary (or multiple) systems in our sample. As stated in Sect. 3.1, we may expect more binaries, as seen in other clusters (see, e.g., Sana & Evans 2011; Sana et al. 2012; Aldoretta et al. 2015). However, only the undetected SB2 binaries would have an impact on the HRD, specifically those systems with relatively separate components. Systems with very close components would have high v sin i values, thus the probability of not detecting them would be small. Therefore, we would expect that only a small number of stars could be affected and our global conclusions would not vary. Moreover, we note that this small number of undetected binaries could be related to the lack of fast rotators shown in Fig. 5 since binary interaction has been proposed to be the origin of the high rotational velocities (de Mink et al. 2014).

Fig. 6 HRD for the sample of O stars in Cygnus OB2 with reliable astrometry using Geneva (Ekström et al. 2012, left) and Bonn (Brott et al. 2011, right) evolutionary stellar tracks and isochrones. Non-rotating models are on the top and rotating models on the bottom. The sample is also subdivided into the spatial groups found by Berlanas et al. (2019). Mean uncertainties are ~0.01 dex in log Teff and ~0.05 dex in log(L∕L⊙).

Fig. 7 Distribution of fast rotating stars (v sin i > 200 km s−1, marked by red dots) in the HRD. Geneva rotating evolutionary models from Ekström et al. (2012).

5.2 HRD versus sHRD

The spectroscopic version of the HRD (henceforth, sHRD) has been widely used in the literature to interpret the evolutionary status of massive populations. It is a distance-independent diagram presented by Langer & Kudritzki (2014) based on the effective temperature and surface gravity where a spectroscopic luminosity is defined as log($\mathcal{L}$) = 4 * log (T_eff) − log g.

For simplicity, in Fig. 8 (top panels), we only show sHRDs for both non-rotating and rotating Geneva models. Comparing the HR diagrams to their respective sHR diagrams we find some notable differences. On the one hand, the whole population seems to be displaced to higher masses in the sHRD. The most massive star is found in the sHRD at around 120 M_⊙ and the less massive one at 18 M_⊙. The age groups seem to be more dispersed although still present: the older group follows the 4–5 Myr isochrones (+2 Myr for Geneva rotating models) reaching masses up to 60 M_⊙ and we find a more concentrated young group at ~3 Myr (+2) with masses up to nearly 85 M_⊙. In addition, we now find an evident gap between the predicted ZAMS and the observed stars above 32 M_⊙, similar to that found in other studies of O stars (see Sect. 5.1 and citations there) and in our discussion of the HRD (now even more conspicuous). Again, the way in which the most massive stars have evolved determines the interpretation of such gap and the possible presence of a recent episode of star formation (~1–1.5 Myr). The advantage of the sHRD is that we do not have to rely on distances. The sample of ten O stars that did not pass the selection criteria for reliable astrometry (RUWE >1.4) were placed in the sHRD (see Fig. 8, bottom panels). The only significant difference we see with the previous results is that this small sample increases the number of stars with 3–5 Myr.

In view of the above findings we conclude that we obtain a similar global information when using the HRD and the sHRD, except that we derive higher stellar masses from the sHRD than from the HRD.

The sHRD is thus useful to interpret the stellar ages and a general evolutionary status for massive stars when no reliable distances (RUWE >1.4) are available for them, although due to the significant difference in mass found here between the two types of diagrams we should be aware of a possible mass overestimation with respect to the HRD. To reconcilethe masses from the two diagrams we would need to increase the stellar luminosity by ~0.3–0.4 dex in the HRD or decrease the spectroscopic luminosity by 0.2–0.3 dex. But that would imply in the first case to increase the distance to Cygnus OB2 to values barely consistent with Gaia DR2 parallaxes (and incompatible with other measurements, such as that of Rygl et al. 2012). In the second case, we would have to increase log g by an amount incompatible with present day atmosphere models⁸. At the same time, T_eff should decrease beyond present uncertainties. We note, however, that the determination of log g depends not only on the model atmospheres, but also on a correct normalization of the observed spectrum (see Simón-Díaz et al. 2020).

Fig. 8 Top: sHR diagrams for the sample of O stars in Cygnus OB2 with reliable astrometry using Geneva (Ekström et al. 2012, considering non-rotating models on the top and rotating models on the bottom, respectively) evolutionary stellar tracks and isochrones. The sample is also subdivided into the stellar groups found by Berlanas et al. (2018b). Bottom: sHR diagrams for the sample of O stars in Cygnus OB2 that have not passed the criterion for reliable astrometry considering the same set of tracks and isochrones. The yellow triangle indicates an upper limit for log g.

Fig. 9 Spatial distribution of the samples in B19 (left) and this paper (right). Coordinates are from the Gaia DR2 catalog. Colors identify the different groups. Group 1 members are represented by blue symbols. Groups 2, 0, and 3 by red, cyan and green symbols. Stars with RUWE >1.4 (only in the right plot) are represented with grey symbols.

Fig. 10 Gaia DR2 proper motion distribution of the samples in B19 (left) and this paper (right). Colors identify the groups as in Fig. 9 and light gray is used for stars with RUWE >1.4.

5.3 Spatial and dynamical distribution from Gaia DR2

In B19, the spatial substructure of the Cygnus OB2 association using parallaxes from Gaia DR2 was explored. The sample used for the analysis comprised known OB members. Using a Bayesian methodology B19 found evidence of two different stellar groups superposed on the association but separated along the line-of-sight by ~400 pc. Taking into account the spectral analyses and the two age bursts of O-type stars found in this work, we explore further the relationship between both groups and their dynamical state. We note that B19 include objects with confirmed OB spectra and RUWE ≤1.4, whereas in the present work we only include the O stars, but independently of their RUWE value. Thus, B19 include 200 objects, whereas we only consider 78.

5.3.1 Spatial distribution

We plot in Fig. 9 the spatial distribution of both samples, namely, B19 and ours. The left panel shows the spatial distribution of the B19 groups: the main groups at 1.35 and 1.76 kpc (Groups 1 and 2, respectively), those intermediate objects that could not be clearly assigned to either group (Group 0) and those that are foreground or background contaminants (Group 3). The distribution is dominated by the large Group 2 that appears splitted in two main concentrations: one to the left (larger right ascension) around the Cyg OB2 #8 and Cyg OB2 #22 systems plus Cyg OB2 #7⁹, and the other one, with less massive members, to the right and centered around Cyg OB2 #15. Only Group 2 shows a strong stellar concentration, whereas the other groups are much more dispersed. Beyond that, it is not possible to spatially separate the groups found by B19.

In the right panel, we present the spatial distribution of our sample (the 78 O stars) where we included the group of stars with poor parallax measurements for which we have stellar parameters (and that are not in the left plot). We see the same groups as in the full sample. The spatial features remain the same as in the full sample, except that now a small concentration of O stars around Cyg OB2 #11 (to the northeast of that around the Cyg OB2 #8 complex) is more apparent. We note that there is only one representative star of the B19 Group 3 (that is composed by foreground and background objects) and it is precisely J20272428+4115458, the B0IV star that we decided to keep in our sample (see Sect. 3.2). Therefore we find that all O stars that we see belong to the Cygnus star forming complex.

5.3.2 Proper-motion distribution

Proper motions are represented in Fig. 10. Again, the B19 sample appears to the left. In that plot star J20322615+4057194 (with μ_α= −5.368 ± 0.085 and μ_δ= −14.083 ± 0.095 mas yr⁻¹) has been excluded to expand the scale. In the B19 sample we see that there is a small difference between Group 1 and Group 2: while their distributions overlap, Group 2 forms the lower envelope, whereas Group 1 defines the upper one. As a result, the mean proper motions of both groups are slightly different. However, statistically the difference is not significant. The core of Group 1 has (μ_α = −2.770 ± 0.450, μ_δ = −3.679 ± 0.732) mas yr⁻¹ and that of Group 2 shows (μ_α = −2.641 ± 0.340, μ_δ = −4.301 ± 0.380) mas yr⁻¹. We note that these numbers do not take individual errors into account and assume that the errors are independent. The difference in right ascension proper motion is thus not significant (0.129 ± 0.564) mas yr⁻¹, whereas the proper motion in declination is larger, but still within the uncertainties (0.622 ± 0.825) mas yr⁻¹. Nevertheless, the difference in the proper motion in declination is clearly deserving of a diligent inspection in the future Gaia DR3¹⁰. To the right of the main concentration, we see a small number of stars with similar proper motions. They are a mixture of early B stars from different groups and although we note that three of them belong to Group 1, they lie far apart from one another on the sky and, thus, they are not physically related. The rest of the stars have proper motions departing from the average of the central objects and can be considered as runaway candidates (except those belonging to Group 3).

On the right side of Fig. 10 we show the same plot, now for the stars in our sample. Star Cyg OB2 #22B (μ_α = −4.628 ± 0.776 and μ_δ = 8.367 ± 0.797 mas yr⁻¹) has also been excluded to expand the scale (we note that it also has a very high RUWE value). Now the difference between Groups 1 and 2 is not so clear because of the low number of targets in Group 1, although they tend towards higher μ_α values than the Group 2 members. Now we introduce the O-stars that could not be assigned to Group 1 or 2 and we see that six of them can be classified as runaway candidates. The runaways candidates are listed in Table 6. For completeness we also list the early B-stars from the original list of B19 that may be considered runaway candidates, although they are not analyzed in this work (we exclude objects included in Group 3 as they cannot be considered runaway candidates).

Also interesting is the fact that all O-type runaway candidates belong to Group 2 or have a poor RUWE value (and could not be assigned to any group), whereas B-type runaways belong mostly to Group 1, with only one member of Group 2. The three stars from Group 2 that we may consider runaway candidates are A46, A15 and A37. The last two are also fast rotators. This might be seen as supporting the relation between fast rotators and runaways (see, e.g., Blaauw 1993; Li 2020), but other five fast rotatorsfrom Group 2 share the proper motion of the majority of the stars. In addition, most runaway candidates belong to the group ofstars with poor RUWE values and none of these candidates is a fast rotator (at least, from its projected rotational velocity). Therefore, in Cygnus OB2 we cannot find a strong association between runaways and fast rotation (which is not surprising in view of the results in Sect. 4.2). A37 has been identified by Kobulnicky et al. (2016) as having a bow-shock, indicative of the star high velocity against the interstellar medium and reinforcing its status as runaway candidate. In the group of runaway candidates with poor RUWE we have Cyg OB2 #10, that was also proposed as such by Caballero-Nieves et al. (2014), although later Aldoretta et al. (2015) list it as “no runaway.” Interestingly, of the four objects with the highest over-luminosity in X-rays found by Rauw et al. (2015)¹¹, three of them (A20, A26¹² and MT91-516) appear on this list, which leads us to wonder whether the interaction of the stellar wind with the surrounding interstellar medium may play a role here. The fourth one is A11, one of the binary systems listed in Table 3 with an early supergiant of type O7.5Ib(f) as primary.

To the left of Fig. 11, we again show the spatial distribution of the O stars in Cygnus OB2, but now with arrows representing their residual proper motion after having substracted the mean proper motion of the corresponding group (for Group 0 we have averaged the means of Groups 1 and 2, and for stars without an assigned group, we used the mean of Group 2). We see that it is impossible to distinguish between the three groups in Cygnus (0, 1 and 2). The kinematical and dynamical status of Cygnus OB2 was carefully studied by Wright et al. (2016), who found the association to be highly substructured, although their spatial scales have been recently challenged by Arnold et al. (2020), who find substructure at smaller scales of 0.5 pc. A comparison of Fig. 9 in Wright et al. (2016) and on the left of our Fig. 11 makes it obvious that there is a similar pattern, although their sample is dominated by F-G stars. As these authors, we find no clear expansion nor contraction pattern in the O-star population of Cygnus OB2.

Fig. 11 Left: spatial distribution and proper motion arrows for the O stars in Cygnus OB2. Runaway candidates are represented in purple. Right: Gaia DR2 proper motion distribution of our sample divided by binary status. Cyan: SB1 systems; purple: SB2 systems; light gray: single stars. We note that the ordinate scale is different as in the previous figures to include Cyg OB2 #22B.

5.3.3 Signatures of possible dynamical ejection events

One of the puzzles in Cygnus OB2 is the lack of strong evidence of past supernova (SN) explosions in spite of the strong concentration of massive stars. The most compelling evidence to date is the presence of the pulsar PSR J2032 + 4127 (Albacete Colombo et al. 2020) that forms a binary system with Cyg OB2 #4B (a.k.a. MT91-213), a Be star (Salas et al. 2015; Ho et al. 2017). This star has a proper motion that differs from the majority of these stars and, thus, it could be considered a runaway candidate. However, a monitoring of the system by Chernyakova et al. (2020) found that the periastron passage occurred in November 2017, as predicted by Ho et al. (2017). Furthermore, Gaia proper motions could be affected by the orbit of Cyg OB2 #4B around PSR J2032 + 4127, although given the mass ratio large effects are not expected. Since we currently have no access to the Gaia individual epochs, we tentatively include Cyg OB2 #4B in Table 6. Other runaway candidates in the field could also be the result of ejection after a SN explosion or after a dynamical encounter and, in fact, all O-type runaway candidate stars listed in Table 6 have an age that is in excess of 2 Myr (as derived from the HRD and the evolutionary tracks), which makes them compatible with any scenario for a possible ejection. However, we point out that the region considered in this paper is relatively small and that evidence for runaways from Cygnus OB2 should consider a larger extension (see e.g., the case of the runaway BD + 43°3654 in Comerón & Pasquali 2007; Wright et al. 2015).

We also explored the distribution of binary systems in the proper motion diagram to see whether they are affected in some way in their proper motions. All of them but three show RUWE values below 1.4. Two of the three stars with RUWE >1.4 are just above this value (Cyg OB2 #1 and Cyg OB2 #27, with values of 1.421 and 1.521, respectively). The outlier is Cyg OB2 #22B, with a RUWE value of 5.365. Interestingly, the distribution of binaries (SB1 and SB2) in the proper motion diagram (see Fig. 11, right) shows that they all are consistent with that of the main concentration. The only exception is again Cyg OB2 #22B¹³. This implies that most probably these systems still keep their original velocities and reinforce the possibility that peculiar proper motion values are the result of a dynamical ejection or a SN kick and are not attributed to the impact of binary motion.

6 Conclusions

Our new spectroscopic census of O stars in Cygnus OB2 provides an updated view of the high-mass stellar population in one of the largest groups of young stars in our Galaxy. The main results from this study are listed as follows.

• I)

  We compiled the most complete spectroscopic census of O stars so far in the association, consisting of a sample of 78 O-type stars. We identified two new binary systems and performed quantitative spectroscopic analysis of the sample (excluding detected SB2 stars) to obtain the distribution of rotational velocities and main stellar parameters.

• II)

  The distribution of rotational velocities is similar to those obtained from analogous studies of early-type stars in other Galactic (and extragalactic) regions, except that we find a short tail of fast rotators, which requires further investigation.

• III)

  In addition, using distances from Gaia DR2 we estimated radii, luminosities, and spectroscopic masses for all isolated stars of our sample with reliable astrometry (RUWE ≤ 1.4). Using the Hertzsprung-Russell diagram the global evolutionary status of the association has been discussed. We find star formation during the last 1–6 Myr, with two main bursts centred roughly at 3 and 5 Myr. A third smaller group of stars at ~1.5 Myr containing the hottest stars of our sample might be identified, but its veracity depends on the possible evolutionary paths followed by these stars.

• IV)

  We obtained similar results when using the HRD and its spectroscopic version (sHRD), except that we derived higher stellar masses from the sHRD than from the HRD. If confirmed, this would be a cautionary remark for the use of the sHRD, although it would still be possible to use it to study the evolutionary status of a stellar population.

• V)

  We investigated the spatial and dynamical distribution of the sample from Gaia astrometry. We find the association to be quite substructured with no evidence of expansion nor contraction pattern in its O-star population. Moreover, we identified nine O stars as runaway candidates, that would support a possible dynamical ejection scenario or a past SN explosion. Proper motions of binary systems also support this possibility since the sample shares the proper motions of the global population, indicating that they have kept their initial dynamical properties.

This work has demonstrated the advantages and difficulties of performing accurate spectroscopic analysis in a very active and obscured region as Cygnus OB2. It shows the precision that can be reached at using spectroscopic parameters for the interpretation of the evolutionary status of a region, that is, revealing internal star-forming bursts that otherwise would be unnoticed. It also shows the importance of the spectral resolution at performing quantitative spectroscopic analysis (mainly for obtaining an accurate distribution of rotational velocities).

The improved view we now have of Cygnus OB2 brings us closer to a complete understanding of the origin and evolution of this association, Cygnus-X, and OB associations, in general. It provides us with a solid background to interpret the results from other studies, such as the forthcoming Cygnus survey of WEAVE, which is the next multi-object spectrograph at the 4.2 m WHT, whose first light is planned by the end of 2020. Over the next few years, the WEAVE survey will deliver high quality ground-based spectroscopy for several hundred high-mass stars in the field. This work will help to optimize and homogenize the sample, adding additional epochs for detecting multiplicity and serving as a template for the other OB associations in the complex. The procedures and results presented here will thus constrain the overall strategy for the spectroscopic analysis of a large amount of OB spectra over the whole region. Together with the upcoming improved Gaia DR3 astrometry, this will allow us to make a much more accurate interpretation of the entire Cygnus-X evolutionary status.

Acknowledgements

This paper is dedicated to the memory of Rebeca Galera Rosillo and her green fireflies. We thank the referee for useful and valuable comments that helped improve this paper. We acknowledge financial support from the Spanish Ministry of Economy and Competitiveness (MINECO) under the grants AYA2012-39364-C02-01, AYA 2015-68012-C2-01, Severo Ochoa SEV-2015-0548 and PGC2018-093741-B-C21/C22 (MICIU/AEI/FEDER, UE).

Appendix A: Best-fitting models

In this appendix, we show an example of the online material available at the Zenodo open-access repository¹⁴. We present the FASTWIND best-fitting model to the observed spectra for one star of the Cygnus OB2 O-type population analyzed in this work (Cyg OB2 #15). In this appendix, we show an example of the online material available at the Zenodo open-access repository where we present the FASTWIND best-fitting models to the observed spectra for the Cygnus OB2 O-type population analyzed in this work.

Fig. A.1 FASTWIND best-fitting model to the observed spectra for Cyg OB2 #15. The observed spectra (in black) is overplotted to the best-fitting model (in red) resulting from the iacob-gbat analysis. He I, He II, and H lines are indicated with solid, dashed, and dotted short vertical lines, respectively.

References

All Tables

All Figures

	Fig. 1 Inverse Spitzer 8 μm image showing the location of the known O-type population in Cygnus OB2 (78 stars represented by green dots). The solid line circle delimits the 1 degree radius area adopted for the association. For reference, the dashed-dotted line circle shows the area considered by Wright et al. (2015) indicating the core of the association.
In the text

	Fig. 2 Example of the new SB2 detections: HeI lines in stars J20301838 + 4053466 (top) and A32 (bottom).
In the text

	Fig. 3 Comparison between the stellar parameters obtained using spectra from the IACOB-OWN (R ≥ 25 000) and the GOSSS surveys (R ~ 2500) for a sample of 37 standard O-type stars in common. The same v sin i values and set of diagnostic lines were used for the spectroscopic analysis. Red arrows indicate upper limits.
In the text

Fig. 4 Comparison of projected rotational velocities resulting from the Fourier Transform (FT) and the Goodness-Of-Fit (GOF1) techniques obtained using the iacob-broad tool. Dashed lines represent a difference of 20 km s−1 or 20% from the 1:1 relation, whichever is the largest. Red squares indicate those stars out of these limits. Different colors indicate the diagnostic lines used for the line-broadening characterization.

In the text

Fig. 5 Distribution of rotational velocities of our final sample of 63 O-type isolated stars in Cygnus OB2. Left: histogram color-coded on the basis of the resolving power of the spectra used for each target. Middle: v sin i histogram showing the comparison between the distribution of rotational velocities found for our sample of O-type stars in Cygnus OB2 and the distribution of O-type stars found in the MW. Vertical black lines indicate the associated Poisson noise. Right: as in middle panel but compared to 30 Dor.

In the text

Fig. 6 HRD for the sample of O stars in Cygnus OB2 with reliable astrometry using Geneva (Ekström et al. 2012, left) and Bonn (Brott et al. 2011, right) evolutionary stellar tracks and isochrones. Non-rotating models are on the top and rotating models on the bottom. The sample is also subdivided into the spatial groups found by Berlanas et al. (2019). Mean uncertainties are ~0.01 dex in log Teff and ~0.05 dex in log(L∕L⊙).

In the text

	Fig. 7 Distribution of fast rotating stars (v sin i > 200 km s−1, marked by red dots) in the HRD. Geneva rotating evolutionary models from Ekström et al. (2012).
In the text

Fig. 8 Top: sHR diagrams for the sample of O stars in Cygnus OB2 with reliable astrometry using Geneva (Ekström et al. 2012, considering non-rotating models on the top and rotating models on the bottom, respectively) evolutionary stellar tracks and isochrones. The sample is also subdivided into the stellar groups found by Berlanas et al. (2018b). Bottom: sHR diagrams for the sample of O stars in Cygnus OB2 that have not passed the criterion for reliable astrometry considering the same set of tracks and isochrones. The yellow triangle indicates an upper limit for log g.

In the text

	Fig. 9 Spatial distribution of the samples in B19 (left) and this paper (right). Coordinates are from the Gaia DR2 catalog. Colors identify the different groups. Group 1 members are represented by blue symbols. Groups 2, 0, and 3 by red, cyan and green symbols. Stars with RUWE >1.4 (only in the right plot) are represented with grey symbols.
In the text

	Fig. 10 Gaia DR2 proper motion distribution of the samples in B19 (left) and this paper (right). Colors identify the groups as in Fig. 9 and light gray is used for stars with RUWE >1.4.
In the text

	Fig. 11 Left: spatial distribution and proper motion arrows for the O stars in Cygnus OB2. Runaway candidates are represented in purple. Right: Gaia DR2 proper motion distribution of our sample divided by binary status. Cyan: SB1 systems; purple: SB2 systems; light gray: single stars. We note that the ordinate scale is different as in the previous figures to include Cyg OB2 #22B.
In the text

	Fig. A.1 FASTWIND best-fitting model to the observed spectra for Cyg OB2 #15. The observed spectra (in black) is overplotted to the best-fitting model (in red) resulting from the iacob-gbat analysis. He I, He II, and H lines are indicated with solid, dashed, and dotted short vertical lines, respectively.
In the text