© ESO 2018

1. Introduction

The chemical composition of the natal molecular cloud is one of the fundamental parameters in the process of massive star formation, as it affects the opacity of the infalling matter, the chemical processes that take place during star formation or the upper mass limit, to mention some (e.g., Zinnecker & Yorke 2007; Tan et al. 2014). The previous enrichment history of the interstellar medium (ISM) through stellar evolution and processes like stellar winds, metal-rich yields from supernovae (SNe) or infall of metal-poor gas, directly affects the chemical composition of the natal molecular cloud (e.g., Langer 2012; Prantzos et al. 2018; Barrera-Ballesteros et al. 2018).

Young open clusters and OB associations have been proposed as places to explore the effects of self-enrichment by stellar winds and SNe (e.g., Spina 2017; Palla 2011; Biazzo et al. 2011a; Simón-Díaz 2010; Cunha & Lambert 1994, 1992) as many of them have a star formation history sufficiently long to have experienced the explosion of their oldest massive members as SNe, contaminating the environment through processes driven by the pollution of the ISM by subsequent generations of massive stars. We may thus expect that subsequent stellar generations show chemical differences. Inversely, the presence or absence of these differences allow us to put constrains on the processes involved in the chemical pollution of the proto-stellar cloud. For example, predictions about the oxygen yields from core-collapse supernovae (CCSN) greatly differ depending on the upper mass limit to become CCSN or on the physics of massive stars included in the models (see e.g., Suzuki & Maeda 2018; Prantzos et al. 2018; and references therein). Unfortunately, the stellar content of many OB associations is difficult to study in detail due to the combined effects of distance and extinction.

A study in Orion (Simón-Díaz 2010) resulted in a chemically homogeneous composition, contrary to previous works claiming evidence of self-enrichment (Cunha & Lambert 1994). The homogeneity has later been confirmed by other independent analysis strategies (see Nieva & Simón-Díaz 2011; Cunha et al. 2012). Moreover, this result was not only found in early-type stars, but also in late-type stars (e.g., Biazzo et al. 2011a,b). But Orion is a modest region in terms of star formation.

Because of its close distance of ∼1.4 kpc (Rygl et al. 2012) and relatively intense star formation, Cygnus OB2 provides an excellent target to study possible chemical composition inhomogeneities through self-enrichment. Massive star formation seems to have peaked around 3 Myr ago (Comerón & Pasquali 2012; Wright et al. 2015; Berlanas et al. 2018) and the distribution of stellar ages extends beyond 20 Myr (Comerón et al. 2016). Although there are no known supernova remnants clearly associated with Cygnus OB2, the current existence of a large number of B-type giants (Berlanas et al. 2018) and seven red supergiants in the region (Comerón et al. 2016) indicates an evolved population whose most massive members have already exploded as SNe. The three Wolf-Rayet stars known in the core of the region and the large number of OB stars (Wright et al. 2015), represent intense stellar wind sources susceptible to produce chemical enrichment across the association. Strong signatures from the decay of ²⁶Al in Cygnus OB2 (Knödlseder et al. 2002), other high-energy detections at γ- and X-rays (Martin et al. 2009, 2010; Bednarek 2003; Butt et al. 2006) or the presence of the young pulsar J2032+4130 in the surroundings (Abdo et al. 2009; Camilo et al. 2009), could have been caused chemical enrichment produced by stellar winds and/or SNe.

In virtue of its rich content in massive stars and an extended star formation history sufficiently long to have included the explosion of its oldest massive members, Cygnus OB2 represents a suitable region to explore possible abundance inhomogeneities. The observed spatial age gradient across its extent (Comerón & Pasquali 2012) might be associated with a chemical composition gradient produced by pollution of the ISM by successive generations of massive stars.

In Sect. 2 we describe the target selection and the observing strategy. In Sect. 3 we present the whole spectroscopic analysis, including the derived equivalent widths, stellar parameters and abundances of a sample of eight early-type slow rotators of Cygnus OB2. In Sect. 4 we discuss our results and their implications. Finally, in Sect. 5, we summarize the main conclusions of this work.

2. Target selection, observations and membership

Late O-type and early B-type stars, with their rich spectrum of silicon and oxygen lines are ideal targets to test possible chemical contamination. These stars should cover the whole association area (1 deg radius centered on Galactic coordinates l = 79.8° and b = +0.8°). In addition, to reach an optimal accuracy in the abundance analysis, we need to analyze stars with effective temperature in the range 18 000–34 000 K or, equivalently, spectral types in the range B2–O9.5 (where silicon ionization balance and O II lines are strong enough, see Simón-Díaz 2010) and low projected rotational velocities (in what follows, we call slow rotators those stars with projected rotational velocities below 80 km s⁻¹) at high resolution and signal-to-noise ratio (S/N). The fraction of early-B dwarfs and giants rotating at these velocities depends on the environment and evolutionary stage of the cluster considered, but from the results in the Milky Way (Dufton et al. 2006; Wolff et al. 2007; Huang & Gies 2006) and 30 Dor (Dufton et al. 2013) we may expect it to be 20–30%, which makes mandatory to previously observe a large sample of OB stars to identify early-type slow rotators.

In this section we describe the observations and various steps we have performed in order to build an optimal list of stars for the purposes of this paper. Briefly we started with intermediate resolution observations of a large list of OB stars in Cygnus OB2 which helped us identifying the late-O and early-B slow rotators. Then we proceeded to gather new high-resolution observations of this optimized sample. In the last step, we benefit from data about parallaxes and proper motions from the second Gaia data release (Gaia DR2, Gaia Collaboration 2018) to asses the membership of the final stellar sample.

2.1. Intermediate resolution spectroscopy and early-type slow rotators selection

In July 2014 we obtained intermediate-resolution spectroscopy of 51 stars in Cygnus OB2 classified as late-O and early-B in the large list from Comerón & Pasquali (2012) using the multi-object AF2/WYFFOS spectrograph at the William Herschel Telescope (WHT). We selected four fields in the region (see Table 1) in order to maximize the number of suitable candidates. Simulations based on a B1V spectrum broadened by different values of υ sin i and different levels of noise indicated that the best trade-off between the required exposure time and enough υ sin i accuracy for our target selection can be achieved with a combination of a resolving power of 5000 and S/N ≥ 50. We thus have chosen the H2400B grating, which provides the required resolution and wavelength coverage, setting the central wavelength to 4500 Å. The spectral region contains the Si III lines at λ4552, λ4568 and λ4575 Å, which are broadened mainly by rotation and are most prominent around the B1V type stars. Finally, in order to improve the statistics of our sample, we added six stars from previous unpublished observations in Cygnus OB2 carried out by A. Herrero in August 2012 at the WHT using the ISIS spectrograph with the R1200B and R1200R gratings (covering all the lines of interest and providing high quality spectra good enough for our purposes).

We determined the projected rotational velocities for the whole stellar sample of 57 stars using iacob_broad, an user-friendly tool for OB-type spectra which is based on a combined Fourier transform plus a goodness-of-fit methodology (Simón-Díaz & Herrero 2007, 2014). We have mainly based the analysis on the Si IIIλ4552 line or the N IIIλ4510 line in the cases in which the Si III line was weak. When both lines were weak, we have used the He Iλ4387, 4471 lines. Typical uncertainties are in the order of 10–20% (Simón-Díaz & Herrero 2014). We found a fraction of OB slow rotators consistent with the expectations, reaching a ∼23% of the targets at υ sin i ≲ 80 km s⁻¹ (see Fig. 1). For four targets, represented with dark gray color, we suggest possible spectroscopic binarity (two of them are reported as binary candidates by Berlanas et al. 2018). We can only suggest their possible binary nature due to their noisy spectra.

The whole list of 57 stars (we note that some of them are too hot for the abundance analysis) with the derived rotational velocities is given in Appendix A. From this sample, and taking into account the required conditions mentioned above, six stars were suitable for our abundance analysis.

Fig. 1. Rotational velocity distribution of a sample of 57 OB-type stars in Cygnus OB2. Dark gray color indicates possible binary candidates. The red dashed line indicates υ sin i = 80 km s−1.

2.2. High resolution and S/N spectroscopy

In July 2015 we obtained high resolution, high S/N spectra of the six selected stars. We obtained ISIS spectra at the WHT using both, the blue and red arms. We chose the H2400B grating for the blue arm, which provides the required resolving power (R ∼ 13 600 at 4500 Å). The spectral coverage of the H2400B grating requires the observation of each star in two spectral ranges centered on 4200 and 4600 Å, in order to cover all the lines of interest. But due to bad weather conditions and in order to get high S/N spectra, the H2400B grating was replaced by the R1200B one for the faintest targets (J20292449+4052599 and J20314605+4043246), which provides a lower but still sufficient resolving power (R ∼ 7500 at 4500 Å by reducing the slit width to 0.7 arcsec) and covers all lines of interest. As the dichroic allows to simultaneously observe the red arm, we used the R1200R grating (R ∼ 9300 at 6500 Å) centered at 5700 and 6500 Å for all the stars, covering the regions of O IIIλ5592, C IVλ5811-14 Å and Hα. Additionally we have included two more stars, J20314965+4128265 and J20334610+4133010, both OB slow rotators that meet all the requirements and were observed by A. Herrero in July 2003 using also ISIS at the WHT with the same configuration (H2400B and R1200R gratings). Figure 2 shows the location of the whole sample of 57 OB stars in the region, and the eight stars finally chosen to carry out the abundance analysis.

The spectra were reduced using IRAF¹ with standard routines for bias and flat-field subtraction, and also for the wavelength calibration. A summary of their basic data is presented in Table 2.

2.3. Gaia-DR2 data

We used the recent data (parallaxes and proper motions) from Gaia DR2 to check whether all the stars of our sample belong to the same region and present no peculiarities in the data that could indicate a different origin or evolutionary history.

Fig. 2. Inverse Spitzer 8 μm image showing the location of the 57 stars (dots) observed in Cygnus OB2, and the selected OB slow rotators (yellow dots) suitable for the abundance analysis. The solid line circle indicates the 1 deg. Cygnus OB2 area adopted in this work, and the four dashed line circles the observed AF2/WYFFOS fields. For reference, the red star indicates the location of the Cyg OB2 #8 Trapezium-like system.

All but one of the stars in our sample with a parallax uncertainty below 0.1 mas (see Table 3) have values consistent with the known Cygnus OB2 distance of 1.4–1.7 kpc (Massey & Thompson 1991; Hanson 2003; Rygl et al. 2012). There is one star (J20272428+4115458) which seems to be at further distance (about 1.8 sigma from the mean parallax value of our sample). However, there is a systematic (but variable) offset up to 0.1 mas that could affect all the reported Gaia DR2 parallaxes (Arenou et al. 2018; Lindegren et al. 2018; Stassun & Torres 2018). Nevertheless, the oxygen abundance obtained for this star is consistent with the mean value obtained for the sample. Other two stars have large parallax uncertainties, but still consistent with the average value.

Focusing on proper motions, they are all consistent with membership in Cygnus OB2. Our stellar sample with proper motion uncertainties below 0.1 mas yr⁻¹ have mean values of pm_RA = −2.5 ± 0.05 mas yr⁻¹ and pm_Dec = −3.85 ± 0.06 mas yr⁻¹.

We therefore conclude that our sample belongs to the normal population of Cygnus OB2 and does not include any runaway star.

3. Spectroscopic analysis

We have followed the same methodology as Simón-Díaz (2010), where the details of the FASTWIND grid used (that includes H I, He I-II, Si II-III-IV and O II) and the considered atomic models are provided. In addition, we have extended the grid to lower gravities for the supergiants analysis.

3.1. Spectral classification

We have reclassified the eight selected stars since the new spectra obtained in this work have, for many of them, a higher quality than in previous studies. We have classified the O-type stars using the MGB tool (Maíz Apellániz et al. 2012). It compares the observed spectra with a grid of O standards (in this work the GOSSS library, see Maíz Apellániz et al. 2016), allowing us to vary spectral type, luminosity class, velocity and resolution until obtaining a best match. For B-type stars we used templates gathered in the framework of the IACOBsweG spectroscopic survey (Simón-Díaz et al. 2015; Negueruela et al., in prep.) to obtain the spectral types. The new classification of the selected stellar sample is presented in Table 2 together with that from other works.

3.2. Equivalent widths

We have carried out an analysis based on equivalent widths (EW) of metal lines, similar to the method used in Orion by Simón-Díaz (2010). We used our own IDL² routines to identify metal lines in the observed spectra, and to measure the EWs and their uncertainties. For each line, a multi-Gaussian fit of the observed line profile is done in a spectral range of λ_o ± 2 max[υ sin iλ_o/c, 0.5λ_o/R] around the laboratory wavelength of the line. The uncertainty is obtained by assuming the continuum at ±(S/N)⁻¹, and comparing the value obtained by means of the Gaussian fitting with the value derived by integrating the line profile. Measured EWs are given in Appendix B.

3.3. Stellar parameters

Both parameters, effective temperature (T_eff) and surface gravity (log g), were obtained by comparing the EW ratio of Si II-III or Si III-IV (depending on the effective temperature of the star) and the wings of the H Balmer lines with FASTWIND stellar atmosphere models (Santolaya-Rey et al. 1997; Puls et al. 2005). We first used the silicon line ratios to obtain the initial values of T_eff for different log g values (e.g., Fig. 3, top). Then, the models obtained by setting each set of these parameters were compared with the wings of the observed H Balmer lines (e.g., Fig. 3, bottom) to iteratively obtain the final values for both parameters.

Fig. 3. Top: Example of the EW ratio of Si III-IV used as Teff indicator with FASTWIND models for the star J20272428+4115458. Black, gray, blue, green and purple lines correspond to log g values of 3.9, 4.0, 4.1, 4.2 and 4.3 dex respectively. The horizontal dashed line indicates the EW ratio of the star and its associated uncertainty is indicated by the dash-dotted lines. Bottom: Example of the H Balmer lines fitting to determine log g for the same star. Only two FASTWIND models are compared with the observed spectrum (black line) for clarity. In red, the Teff = 30 000 K, log g = 3.9 dex model. In green, the Teff = 31 000 K, log g = 4.2 dex model. Final fits can be seen in Figs. C.1–C.8.

We have used the ratios EW(Si IVλ4116)/EW(Si IIIλ4552) and EW(Si IIλ4130)/EW(Si IIIλ4552) as T_eff indicators, depending on the temperature of the star, because these are the strongest, unblended lines typically present in the stellar spectrum (at the spectral types of our targets and in the range of ∼4000–5000 Å). On the other hand, due to the good agreement between T_eff determined using Si III-IV and the He I-II ionization balance found by Simón-Díaz (2010), the He I-II ionization equilibrium has been used for the star J20292449+4052599 because only Si III lines were suitable for the stellar parameter determination in the observed spectrum.

Derived T_eff uncertainties from EWs of Si II-III-IV lines are, for all our stars, lower than 500 K. For the surface gravity we obtain an accuracy better than 0.1 dex. Taking into account the additional uncertainty sources (like continuum normalization) we follow the reasoning by Sabín-Sanjulián et al. (2017) and adopt in all cases errors of ±500 K and ±0.1 dex for T_eff and log g, respectively. Derived values are shown in Table 4.

3.4. Silicon and oxygen abundances

We applied the curve-of-growth method to derive silicon and oxygen abundances, considering a grid of HHeSiO FASTWIND models for a set of stellar parameters (T_eff and log g). For each star, the abundance is obtained for different values of microturbulence and for each line considered. Then, we determined the final abundance value considering the microturbulence that, for all lines, gives a similar value (see Fig. 4). More details can be found in, for example, Kilian (1992) or Simón-Díaz (2010). We have used all the lines available in our spectra, except those too faint or blended with other lines and that can not be properly separated by our IDL routines (see Sect. 3.2) at the resolving power of the analyzed spectra. This is, for example, the case for the blend of Si IV λ4089 and O II λ4089 lines. Other cases that we have also carefully checked are Si IV λ4631, and Si II λ6347, λ3856 since they could be blended with N II λ4631, Mg II λ6347, and O II λ3856, respectively.

Fig. 4. Example of the silicon (top) and oxygen (bottom) abundance versus EW diagrams for the star J20272428+4115458. All the observed Si III-IV and O II lines are indicated in the plots. The red dashed line indicates the weighed mean abundance (fit to the data), and the gray dash-dotted lines the standard weighed deviation.

In the FASTWIND oxygen grid, only lines from one ionization state were available (O II), so we could not use the oxygen ionization equilibrium. Therefore, we adopted for the oxygen abundance analysis the same stellar parameters as in the silicon analysis, including microturbulence. Although Simón-Díaz (2010) used the microturbulence obtained from the O II lines, we are more limited in the number of suitable spectral lines and decided to use the microturbulence derived from the silicon lines to avoid spurious results. In Appendix B the results of the silicon and oxygen abundance analysis obtained for each line are detailed, where the uncertainties related to the errors in EWs (Δϵ), the standard weighted deviation in the line-by-line abundances (Δϵ(σ)) and the uncertainties associated with the microturbulence (Δϵ(ξ)) are indicated.

We observed no clear dependence between the derived abundances and T_eff. However, due to the strong sensitivity of some line EWs to temperature we check the effect of varying T_eff in ±500 K (taking the corresponding gravity variation into account). In the case of silicon, the presence of two ionization stages already restricts the abundance range. In contrast, in the case of oxygen, only lines from one ionization state were available, and therefore the derived abundance is quite sensitive to small changes in temperature. Simón-Díaz (2010) tested this effect in a sample of early dwarf B stars. He showed that the minimum variation in oxygen abundance occurs for T_eff ∼ 27 000 K, increasing toward higher and lower effective temperatures. We find the same behavior in our study.

4. Results and discussion

4.1. Chemical composition

Table 5 gives the microturbulence (ξ) and final abundance values (where ϵ_x = 12 + log(X/H)) derived for our stellar sample. It also shows the uncertainties in the derived abundances associated with the uncertainties in microturbulence (Δϵ(ξ)), the standard weighted deviation in the line-by-line abundances (Δϵ(σ)) and the impact that a modification of T_eff by ±500 K and log g by ±0.1 dex have on the derived abundances (Δϵ(T_eff,  logg)).

Figure 5 represents the derived silicon and oxygen abundances as a function of Galactic longitude for our sample of stars. We find chemical homogeneity for both elements, and derive mean silicon and oxygen values of 7.53 dex and 8.65 dex, respectively, that represent the mean abundance values of Cygnus OB2. The standard deviation for the silicon abundance is 0.08 dex and for the oxygen abundance 0.12 dex. The intrinsic uncertainties (related to the microturbulence, dispersion line-by-line and stellar parameters) are on average on the order of the dispersion in abundances (around 0.1 dex for both silicon and oxygen). The larger dispersion obtained for oxygen, in comparison with silicon, is partially due to the fact that only one ionization state is available for the oxygen analysis, while two ionization states are available for silicon (Si II-III or Si III-IV, depending on the temperature).

Fig. 5. Derived silicon (left) and oxygen (right) abundance versus Galactic longitude in degrees. Dashed black line indicates the mean abundance value and dash-dotted red line the best fit to the data. Vertical lines indicate abundance errors, that include uncertainties related to the stellar parameters, microturbulence and the standard weighed deviation in the line-by-line abundances.

Thus our results indicate a homogeneous composition for our stellar sample, without evidence of a dependence on the Galactic longitude (which, according to Comerón & Pasquali (2012) or Berlanas et al. (2018) correlates with stellar age). Appendix C shows the individual spectra with the best fitting FASTWIND model including H I, He I-II, Si II-III-IV and O II.

4.2. Self-enrichment scenario

In order to see whether we can provide some constraints on the metallicity enrichment within star-forming clouds in spite of the small variations obtained for the oxygen and silicon abundances, we estimate in this section the enrichment produced by stellar winds and SNe that may have polluted the region in the last 30 Myr. This is roughly the time span during which a cluster would contain stars hotter than 20 000 K according the evolutionary tracks with and without rotation by Ekström et al. (2012) and Limongi & Chieffi (2018).

In the following we assume an Initial Mass Function between 0.1 and 120 M_⊙ with a slope α = −2.3, as in Comerón et al. (2016). We also assume three star formation histories (SFH): (1) An Instantaneous Burst (IB) having taken place 16 Myr ago (corresponding roughly to the age when the number of red supergiants (RSGs) is maximum along the IB evolution for solar metallicity models by Ekström et al. 2012). This scenario shows the enrichment due to stars more massive than 15 M_⊙ and it is used as an estimate of the maximum possible enrichment speed (see below). (2) A stair-like SFH with a constant value Φ₁ in the last 6 Myr plus a constant value ${\mathrm{\Phi }}_2=\frac{1}{6}{\mathrm{\Phi }}_1$ for ages larger than 6 Myr (see again Comerón et al. 2016); in our case we have extended Φ₂ up to 46 Myr which is the maximum age where RSGs (defined as stars with logL/L_⊙ ≳ 4 and T_eff ≲ 4000 K) would be present in a system, hence the lower limit where Comerón et al. (2016) SFH inference is applicable. And (3) a linearly increasing SFH from $\mathrm{\Phi }(t_{\mathrm{now}}-30\mathrm{Myr})=\frac{1}{6}\mathrm{\Phi }(t_{\mathrm{now}})$–Φ(t_now), where t_now is the present moment and which is a linear approximation of the previous case; this SFH extends down to 35 Myr (larger ages produce negative values). We have used yields and stellar lifetimes from Limongi & Chieffi (2018), solar metallicity evolutionary models³ (see also Prantzos et al. 2018), which include three initial rotational velocities υ_rot (0, 150 and 300 km s⁻¹), and two SNe scenarii: a first one where all stars produce SNe explosions, and a second one where stars more massive than 25 M_⊙ suffer a direct collapse to black holes. In all cases the enrichment due to stellar winds in the pre-SNe phase (hydrostatic enrichment) has been taken into account, although we have assumed that the hydrostatic enrichment is released to the ISM at the end of the evolution of the star together with the SNe contribution. However, stellar winds produce a minimal effect in the abundance, when defined as a ratio (see below). We note that the variation of the yields with the metallicity has not been considered; this choice is justified a posteriori given the small variation of the total abundance. For all SFHs we have followed the evolution of the total, H, O and Si masses released to the ISM as the stellar populations evolves. As usual in synthesis models computations, the result are normalized to total mass of ever formed stars (i.e., ${\int }_{t_{\text{SFH\hspace{0.17em}begins}}}^{t_{\text{now}}}{\Phi }(t)\text{\hspace{0.17em}d}t=1M_{\odot }$). Such choice allows to compare the different SFHs self consistently since the total normalized mass of the ejected material is equal in all cases (we would like to point out that normalized mass refers to mean mass values of the underling possible distribution (see Cerviño & Luridiana 2006) because to make choice of an absolute value of the total mass implies to provide an study of the related (sampling) distribution, which is outside the scope of this paper⁴). The only differences among the SFHs are the speed in which such material is released, being the IB case the faster case since all stars have equal ages and the only delay is due to the stellar evolution lifetimes.

In order to obtain the chemical enrichment, we assumed a close system with a global efficiency of 10%, being the efficiency the ratio between total amount of gas ever transformed into stars (the integral over time of the SFH up to present time) and the total mass of the system. We note that efficiency is enough to follow the enrichment evolution of ISM if sampling effects are neglected, since the absolute value of the total mass cancels out when the O/H is computed. As a final note, the initial metallicity assumed in computations and given by the tracks, log(O/H)+12 = 8.70 dex, is larger than the possible oxygen abundance of Cygnus OB2. As a result, rather than the nominal values of the abundance, we are interested in the abundance variation and how the slope of the chemical enrichment varies with time and, in particular, in the last 20–25 Myr when the present O-type stars were formed.

The result of the chemical evolution variation is shown in Fig. 6 for the case of O/H (results for Si/H are similar, although with a slightly narrower range of variation). We note that since we are working with the O/H ratio, the enrichment is driven by the populations where (a) the O/H of the ejected material differs significantly from the initial one, and (b) their overall contribution is large in comparison with the other stars which release material to the ISM. Such situation is only meet by the SNe contribution. Actually, there is almost no abundance variation in the IB case where stars more massive than 25 M_⊙ collapse to black holes (red and blue solid lines) up to the age when star less massive than 25 M_⊙ explode as SNe (green and black solid lines). As a secondary effect, the possible enrichment saturates with time, since the O/H of the ejected material approaches to the initial value when the initial stellar mass decreases. Actually, all models reach a steady state value of the enrichment when extended to larger (future) ages, provided the SFHs stops at t_now. Such steady state is reached around 25 Myr after the onset of the burst in the IB case, and it is shown in the figure by an horizontal arrow for the extreme case of maximum enrichment. We stress that it does not means that there is no more oxygen released to the ISM, but that such oxygen enrichment is compensated by the hydrogen enrichment. The figure shows that differences between evolutionary models with and without rotation is not larger than 0.02 dex (5%) beyond present-day accuracies in the spectral analysis. On the other hand, the scenario where all stars produce a SNe explosion (red and blue solid lines, with a maximum possible enrichment of 0.14 dex) and that where only stars less massive than 25 M_⊙ produce a SNe (green and black solid lines, with a maximum enrichment of 0.05 dex) could be distinguished with techniques like those used in the present work, provided that we extend our sample and/or increase the S/N of the spectra. The only quantity we can constrain is the star forming efficiency. With values significantly larger than the one adopted here (10%) we would see metallicity variations in our sample. Finally, if we would firmly constrain the age of the different stars, it would be much easier to discriminate between different SFHs. Gaia DR2 and future DR3 data may be unvaluable to this aim. However, given that the different scenarios produce metallicity variations which are similar to the uncertainties in abundance determination of the current analysis, we conclude that the massive star forming activity during the last 20–25 Myr in Cygnus OB2 has not caused any detectable increase in its O/H and Si/H abundance, within the limitations of the present study, in agreement with the findings in Orion by Simón-Díaz (2010) or Biazzo et al. (2011a).

Fig. 6. Time evolution of the logO/H variation obtained from solar metallicity chemical evolution models, using three different SFHs (an IB and two extended SFH, both with a difference of a factor of six between the actual SF value and the value 30 Myr ago, but with different shapes), and two SNe cases: one where all massive stars produce a SNe explosion and another where stars more massive than 25 M⊙ suffer a direct collapse to black holes (where only the IB case is shown). In all cases a star formation efficiency of 10% is assumed. The arrow in the top right shows the maximum enrichment reached by the three SFHs.

4.3. The Galactic context

Our abundance results are in agreement with the values obtained in Orion OB1 for a sample of early-B type stars (Simón-Díaz 2010), and also with those obtained for a sample of unevolved early B-type stars in the solar neighborhood (Przybilla et al. 2008). Both studies result in a small dispersion of abundances, indicating a homogeneity of oxygen and silicon among B-type stars.

Regarding oxygen, we find a mean abundance value slightly below the one obtained in both Orion OB1 and the solar vicinity. In order to check whether they are consistent with the oxygen radial abundance gradient of the Milky Way, we have compared these values with the work of Esteban et al. (2005). They presented homogeneous data of eight H II regions with Galactocentric distances (R_G) between 6.3 and 10.4 kpc using recombination lines and assuming the Sun at 8 kpc from the Galactic center. Figure 7 shows the spatial distribution of the oxygen abundance as a function of R_G for the considered Galactic H II regions and the mean values obtained for Orion OB1, the solar vicinity, and Cygnus OB2. All of them are consistent with the observed abundance gradient derived from H II regions in the Milky Way and its dispersion. Although one would expect a large abundance value for Cygnus OB2 due to its shorter Galactocentric distance, the average value is consistent within uncertainties with the Orion OB1, the solar neighborhood and the H II regions. Our conclusion is that the observed scatter associated with individual measurements of Milky Way H II regions at a given Galactocentric distance is at least of the same order than within a cluster or association.

Fig. 7. Spatial distribution of the oxygen abundance as a function of Galactocentric distances for Galactic H II regions (red dots), the mean value obtained for Orion OB1 (dark blue dot), the mean value obtained for the solar vicinity (green dot) and the mean value obtained for Cygnus OB2 (light blue dot). The red dashed line represents the least-squares linear fits to the H II regions. The Sun is assumed at 8 kpc.

5. Summary and conclusions

We have presented the spectroscopic analysis of eight late O, early B stars in Cygnus OB2 using the curve-of-growth method to derive silicon and oxygen abundances. We have used high S/N, high resolution spectra obtained with ISIS at the WHT, and FASTWIND stellar atmosphere models. Our targets are early-type slow rotators (υ sin i < 80 km s⁻¹) with temperatures below 34 000 K selected from a list of 57 spectroscopically classified OB stars in Cygnus OB2. They have late O or early B spectral type so that we can use their Si II-III-IV lines for our abundance analysis, and are distributed all across the association area. According to Gaia data there is no evidence of a peculiar behavior indicating a runaway origin of our stellar sample. Thus, the silicon and oxygen abundances of the targeted stars are considered to reflect those of the cloud out of which the stars formed.

The abundance analysis was based on the EW of silicon and oxygen lines and the Balmer line profiles. We first obtain the fundamental parameters (T_eff and log g) by comparing the EW ratio of Si II-III or Si II-IV, depending on the temperature of each star, and the observed wings of H Balmer lines with FASTWIND models. Then, using these parameters, the abundance was computed for different values of microturbulence and for each line considered. Finally, we determined the final abundance value considering the microturbulence that, for all lines, gives the same abundance value.

Both silicon and oxygen analyses suggest a homogeneous stellar chemical composition in Cygnus OB2. We find a mean abundance value of 7.53 dex for silicon and 8.65 dex for oxygen, with a standard deviation of 0.08 dex and 0.12 dex, respectively.

We roughly estimate the contribution to oxygen enrichment of SNe that may have exploded in the region, concluding that the original cloud from which stars formed in the last 20–25 Myr in Cygnus OB2 had a homogeneous metallicity and has not been significantly enriched during its lifetime as a consequence of its star forming activity in agreement with the findings in Orion OB1 and the solar vicinity. All these studies are consistent with the spatial abundance distribution for H II regions in the Milky Way. The oxygen abundance scatter based on H II regions at a given Galactocentric distance is of the same order than within a cluster or association. Therefore, abundance gradients at the local scale, like Cygnus OB2 or Orion OB1, are not expected. We estimate that the effect of self-enrichment by stellar winds and SNe is beyond the accuracy of our analyses, except for the star forming efficiency that cannot be significantly larger than the value of 10% adopted here.

Our computations illustrate the difficulties of making quantitative predictions on the level of enrichment expected in a massive association like Cygnus OB2, and show in any case that such level should be small as compared to the accuracies currently attainable.

Acknowledgments

We thank the referee for useful and valuable comments that helped improve this paper. We also thank F. Najarro, M. A. Urbaneja, and D. J. Lennon for helpful discussions and comments. We thank J. Maíz Apellániz for advice on the use of MGB and L. López-Martín for advice on spectral data reduction. We acknowledge financial support from the Spanish Ministry of Economy and Competitiveness (MINECO) under the grants AYA 2015-68012-C2-01 and Severo Ochoa SEV-2015-0548 and the Gobierno de Canarias under grant ProID-2017010115. AP acknowledges support from Sonderforschungsbereich SFB 881 “The Milky Way System” (subproject B5) of the German Research Foundation (DFG).

References

Appendix A: Rotational velocities

Table A.1 shows the whole list of stars observed in July 2014 and August 2012 at the WHT. Names, spectral types and the derived rotational velocities are indicated.

Appendix B: Individual targets

Tables B.1–B.8 show the detailed results of the silicon and oxygen abundance analysis line-by-line for each star of the sample. Abundances (ξ(x)) are given in units of 12 + log(X/H). We also indicate the final weighed mean value ($\overline{ϵ}$) and the uncertainties related to the errors in EW measurements (Δϵ), the standard weighed deviation (Δϵ(σ)) and the uncertainties associated with the microturbulence (Δϵ(ξ)).

Appendix C: Fitting line models

Figures C.1–C.8 show, for each selected sample star, several H, He, Si and O spectral lines and the best fitting model derived from our abundance analysis.

Fig. C.1. Best-fit HHeSiO model (Teff = 33 000 K, log g = 4.0 dex, ξ(Si) = ξ(O) = 3.0 km s−1, log(Si/H) = −4.20 dex, log(O/H) = −3.30 dex) to the observed spectrum of J20314965+4128265. Si and O lines used in the analysis are indicated with dark and light blue vertical lines, respectively. H and He lines are indicated with dotted and solid black vertical lines for reference.

Fig. C.2. Best-fit HHeSiO model (Teff = 32 000 K, log g = 3.6 dex, ξ(Si) = ξ(O) = 7.5 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.30 dex) to the observed spectrum of J20295701+4109538. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.

Fig. C.3. Best-fit HHeSiO model (Teff = 30 000 K, log g = 3.9 dex, ξ(Si) = ξ(O) = 5.2 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.30 dex) to the observed spectrum of J20272428+4115458. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.

Fig. C.4. Best-fit HHeSiO model (Teff = 30 000 K, log g = 3.2 dex, ξ(Si) = ξ(O) = 17.0 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.40 dex) to the observed spectrum of J20334610+4133010. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.

Fig. C.5. Best-fit HHeSiO model (Teff = 30 000 K, log g = 4.0 dex, ξ(Si) = ξ(O) = 5.6 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.30 dex) to the observed spectrum of J20314605+4043246. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.

Fig. C.6. Best-fit HHeSiO model (Teff = 29 000 K, log g = 3.7 dex, ξ(Si) = ξ(O) = 6.8 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.30 dex) to the observed spectrum of J20292449+4052599. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.

Fig. C.7. Best-fit HHeSiO model (Teff = 19 000 K, log g = 3.8 dex, ξ(Si) = ξ(O) = 5.9 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.30 dex) to the observed spectrum of BD+40 4193. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.

Fig. C.8. Best-fit HHeSiO model (Teff = 18 300 K, log g = 2.20 dex, ξ(Si) = ξ(O) = 13.3 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.40 dex) to the observed spectrum of BD+40 4210. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.

All Tables

All Figures

	Fig. 1. Rotational velocity distribution of a sample of 57 OB-type stars in Cygnus OB2. Dark gray color indicates possible binary candidates. The red dashed line indicates υ sin i = 80 km s−1.
In the text

Fig. 2. Inverse Spitzer 8 μm image showing the location of the 57 stars (dots) observed in Cygnus OB2, and the selected OB slow rotators (yellow dots) suitable for the abundance analysis. The solid line circle indicates the 1 deg. Cygnus OB2 area adopted in this work, and the four dashed line circles the observed AF2/WYFFOS fields. For reference, the red star indicates the location of the Cyg OB2 #8 Trapezium-like system.

In the text

Fig. 3. Top: Example of the EW ratio of Si III-IV used as Teff indicator with FASTWIND models for the star J20272428+4115458. Black, gray, blue, green and purple lines correspond to log g values of 3.9, 4.0, 4.1, 4.2 and 4.3 dex respectively. The horizontal dashed line indicates the EW ratio of the star and its associated uncertainty is indicated by the dash-dotted lines. Bottom: Example of the H Balmer lines fitting to determine log g for the same star. Only two FASTWIND models are compared with the observed spectrum (black line) for clarity. In red, the Teff = 30 000 K, log g = 3.9 dex model. In green, the Teff = 31 000 K, log g = 4.2 dex model. Final fits can be seen in Figs. C.1–C.8.

In the text

	Fig. 4. Example of the silicon (top) and oxygen (bottom) abundance versus EW diagrams for the star J20272428+4115458. All the observed Si III-IV and O II lines are indicated in the plots. The red dashed line indicates the weighed mean abundance (fit to the data), and the gray dash-dotted lines the standard weighed deviation.
In the text

	Fig. 5. Derived silicon (left) and oxygen (right) abundance versus Galactic longitude in degrees. Dashed black line indicates the mean abundance value and dash-dotted red line the best fit to the data. Vertical lines indicate abundance errors, that include uncertainties related to the stellar parameters, microturbulence and the standard weighed deviation in the line-by-line abundances.
In the text

Fig. 6. Time evolution of the logO/H variation obtained from solar metallicity chemical evolution models, using three different SFHs (an IB and two extended SFH, both with a difference of a factor of six between the actual SF value and the value 30 Myr ago, but with different shapes), and two SNe cases: one where all massive stars produce a SNe explosion and another where stars more massive than 25 M⊙ suffer a direct collapse to black holes (where only the IB case is shown). In all cases a star formation efficiency of 10% is assumed. The arrow in the top right shows the maximum enrichment reached by the three SFHs.

In the text

Fig. 7. Spatial distribution of the oxygen abundance as a function of Galactocentric distances for Galactic H II regions (red dots), the mean value obtained for Orion OB1 (dark blue dot), the mean value obtained for the solar vicinity (green dot) and the mean value obtained for Cygnus OB2 (light blue dot). The red dashed line represents the least-squares linear fits to the H II regions. The Sun is assumed at 8 kpc.

In the text

	Fig. C.1. Best-fit HHeSiO model (Teff = 33 000 K, log g = 4.0 dex, ξ(Si) = ξ(O) = 3.0 km s−1, log(Si/H) = −4.20 dex, log(O/H) = −3.30 dex) to the observed spectrum of J20314965+4128265. Si and O lines used in the analysis are indicated with dark and light blue vertical lines, respectively. H and He lines are indicated with dotted and solid black vertical lines for reference.
In the text

	Fig. C.2. Best-fit HHeSiO model (Teff = 32 000 K, log g = 3.6 dex, ξ(Si) = ξ(O) = 7.5 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.30 dex) to the observed spectrum of J20295701+4109538. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.
In the text

	Fig. C.3. Best-fit HHeSiO model (Teff = 30 000 K, log g = 3.9 dex, ξ(Si) = ξ(O) = 5.2 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.30 dex) to the observed spectrum of J20272428+4115458. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.
In the text

	Fig. C.4. Best-fit HHeSiO model (Teff = 30 000 K, log g = 3.2 dex, ξ(Si) = ξ(O) = 17.0 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.40 dex) to the observed spectrum of J20334610+4133010. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.
In the text

	Fig. C.5. Best-fit HHeSiO model (Teff = 30 000 K, log g = 4.0 dex, ξ(Si) = ξ(O) = 5.6 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.30 dex) to the observed spectrum of J20314605+4043246. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.
In the text

	Fig. C.6. Best-fit HHeSiO model (Teff = 29 000 K, log g = 3.7 dex, ξ(Si) = ξ(O) = 6.8 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.30 dex) to the observed spectrum of J20292449+4052599. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.
In the text

	Fig. C.7. Best-fit HHeSiO model (Teff = 19 000 K, log g = 3.8 dex, ξ(Si) = ξ(O) = 5.9 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.30 dex) to the observed spectrum of BD+40 4193. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.
In the text

	Fig. C.8. Best-fit HHeSiO model (Teff = 18 300 K, log g = 2.20 dex, ξ(Si) = ξ(O) = 13.3 km s−1, log(Si/H) = −4.50 dex, log(O/H) = −3.40 dex) to the observed spectrum of BD+40 4210. Sort lines mark the positions of Si, O, H and He lines, as in Fig C.1.
In the text