The massive star population of Cygnus OB2 Free

Abstract

1 INTRODUCTION

Cygnus OB2 is the nearest region to the Sun that can lay claim to being amongst the most massive clusters or associations in our galaxy (e.g. Knödlseder 2000; Hanson 2003; Wright & Drake 2009). The region is known to harbour many tens of O-type stars and hundreds of OB stars (e.g. Massey & Thompson 1991; Comerón et al. 2002; Kiminki et al. 2007), though estimates of the total mass of the association from both high- and low-mass stars range over an order of magnitude from (2–10) × 10⁴ M_⊙ (Knödlseder 2000; Hanson 2003; Drew et al. 2008; Wright et al. 2010b). The association is itself the most obvious example of recent star formation in the massive Cygnus X complex, a region of ongoing star formation (Reipurth & Schneider 2008) that spans ∼200 pc at a distance of 1.4 ± 0.08 kpc (Rygl et al. 2012). The size and proximity of Cyg OB2 therefore make it a template for studies of distant, massive star clusters and OB associations in the Milky Way and in nearby galaxies.

The high-mass stellar content of Cyg OB2 has been well studied over the past half century since the association was first identified by Münch & Morgan (1953) and a number of mid- to late-O stars within it were classified by Johnson & Morgan (1954). Further members were identified and classified by Morgan, Johnson & Roman (1954) and Schulte (1956, 1958). While many of these classifications were revised over the following decades (e.g. Walborn 1973), it was not until the work of Massey & Thompson (1991) that a more extensive spectral survey of the region was conducted, resulting in the identification of 42 O-type stars and 26 B-type stars. A decade later the number of known OB stars increased again with the works of Comerón et al. (2002) and Hanson (2003) following the photometric study of Knödlseder (2000). Since then the number of very massive O-type stars has not increased significantly, with the more recent and deeper surveys uncovering predominantly B-type stars (e.g. Kiminki et al. 2007; Comerón & Pasquali 2012; Kobulnicky et al. 2012). This apparent dearth of new O-type stars discovered in Cyg OB2, despite deeper observations, particularly in the near-infrared (near-IR), and through spectroscopic observing campaigns, may suggest that our census of massive O-type stars is nearing completion.

The age of Cyg OB2 has proved to be a complex subject, due in no small part to the large spatial size of the association and the lack of any clearly recognizable limit to its extent. Hanson (2003) compiled a sample of 85 OB stars in Cyg OB2, from which the presence of a number of mid O-type dwarf stars and high-luminosity blue supergiants led her to conclude that the age of Cyg OB2 could not be much older than a few million years, and settled on an age of 2 Myr, with a spread of 1 Myr. The evolved giant and supergiant stars present were argued by Hanson (2003) to be contaminating non-members, brought into the association as the search for cluster members was extended over a larger area (e.g. Comerón et al. 2002). However, studies of A-type stars in Cyg OB2 have revealed a population of 5–7 Myr old stars towards and to the south of the association (Drew et al. 2008), while X-ray studies of low-mass stars in the association suggest an age of 3–5 Myr (Wright et al. 2010b). If the OB stars in Cyg OB2 are younger than the low-mass stars it could suggest that the massive stars formed after the low-mass stars, a feature observed in some other young clusters (e.g. Povich & Whitney 2010) that has interesting implications for theories of massive star formation and feedback.

The initial mass function (IMF) of high-mass stars provides a powerful diagnostic of the star formation process and is also therefore an area of considerable interest, particularly as to whether the IMF varies between OB associations and star clusters (e.g. Bastian, Covey & Meyer 2010), particularly given recent suggestions that Cyg OB2 has always been a low-density association and did not evolve from a star cluster. Massey & Thompson (1991) found the high-mass IMF in Cyg OB2 to be considerably flatter than found for other regions in the galaxy, though Kiminki et al. (2007) found it to be much steeper than the canonical Salpeter IMF (Salpeter 1955) from an expanded sample, while Wright et al. (2010b) found the lower mass IMF to be in agreement with the Salpeter value.

Work over the past decade has assembled the data needed to form a much fuller view of the OB population in Cyg OB2, thanks to a combination of new spectroscopic surveys of both primary OB stars and their binary companions (e.g. Kiminki et al. 2007; Kobulnicky et al. 2012), improvements in the stellar effective temperature scale that consider line-blanketing and non-LTE (local thermodynamic equilibrium) effects (Martins, Schaerer & Hillier 2005), and the availability of revised stellar evolution models for massive stars that take into account the important effects of stellar mass-loss and rotation (Ekström et al. 2012). Weidner & Vink (2010) found that masses derived from rotating stellar evolutionary models were in much better agreement with dynamical masses for stars in binary systems than were masses derived from non-rotating stellar models, which may go some way towards addressing previous reports of a non-standard IMF in Cyg OB2. Finally, recent high-spatial resolution optical and near-IR Galactic plane surveys (Drew et al. 2005; Lucas et al. 2008; Groot et al. 2009) provide the necessary photometry to resolve close and blended stars in the commonly used catalogues of Massey & Thompson (1991) and the 2 Micron All Sky Survey (2MASS; Skrutskie et al. 2006) – both of which suffer from some degree of blending that is occasionally ignored by users of the data.

Our objective for this work is to bring together all available information on the massive star population in Cyg OB2, gathering spectroscopy from across the literature and selecting the best available photometry for these sources. Our hope is that this will allow a quantitative analysis of the star formation history, IMF, and total mass of Cyg OB2. In Section 2 we outline the compilation of spectroscopic and photometric data from the literature for our census of 169 OB stars, almost doubling the sample of 85 OB stars studied by Hanson (2003). In Section 3 we use this data to calibrate the extinction towards Cyg OB2 and verify the extinction law for this sightline that is then used to derive individual extinctions towards the massive stars. In Section 4 we compile a Hertzsprung–Russell (HR) diagram for the population of massive stars, from which individual stellar masses and ages are derived. In Section 5 we study the age distribution of these stars to infer the star formation history of the association and in Section 6 we study their mass function. In Section 7, we use these results to estimate the total stellar mass of the association, and finally in Section 8 we summarize our conclusions.

2 COMPILATION OF THE SAMPLE

In this section we compile the available photometry and spectroscopy from the literature. Fig. 1 shows the distribution of known OB stars in Cyg OB2 and the wider Cygnus X region centred on the densest part of the association. One of the difficulties with compiling a census of Cyg OB2 is that as an association and not a centrally condensed cluster it is not immediately obvious where to set the limits of such a study. We focus on a 1-deg² circular area centred on Cyg OB2 #8 – the trapezium of O stars at RA 20^h33^m16^s, Dec. +41°18′45′′ long regarded as being at the heart of the association (e.g. Hanson 2003; Vink et al. 2008). This is equivalent to a radius of ∼14 pc at the distance of Cyg OB2. This choice has a number of advantages, most notably that it includes the clear overdensity that constitutes Cyg OB2, but does not extend too far into the ‘field’ population of OB stars that exists across Cygnus X (e.g. Comerón et al. 2008).

Figure 1.

Mid-IR 12 μm image of the Cygnus X giant molecular cloud compiled from Wide-Field Infrared Survey (Wright et al. 2010a) data and spanning 4 × 3 deg² (100 × 75 pc² at the distance of Cyg OB2). Known O, B, and Wolf–Rayet stars are shown as red, green, and yellow dots, respectively (compiled from the literature as outlined in Section 2.1). The white circle shows the 1 deg² region studied in this work, which is centred on the ‘trapezium’ of O stars known as Cyg OB2 #8, itself marked with a white diamond, and as broadly cospatial with the Cygnus OB2 Chandra Legacy Survey (Wright et al. 2014a).

Open in new tabDownload slide

2.1 Spectroscopic data

The census was compiled by searching the literature for all spectroscopically classified stars of spectral type B5 or earlier, a range of spectral types that can provide a solid diagnosis of the age, mass function, and spatial extent of the association. This resulted in 167 stars (22 with known binary companions) from 18 different publications. To limit any potential discrepancies arising from using stars classified by different authors wherever possible, we use classifications made by comparison with the stellar atlas of Walborn & Fitzpatrick (1990). Fortunately, this includes the majority of known OB stars in Cyg OB2. When a star is classified in this way in multiple papers, the practice here is to adopt the classification made at the highest spectral resolution.

The majority of stars come from the studies of Massey & Thompson (1991, 50 stars, or 26 per cent of both primary and secondary stars) and Kiminki et al. (2007, 78 stars or 41 per cent), which show excellent agreement with each other when their classifications overlap. A further 27 stars (14 per cent) were taken from a series of papers by the latter team of authors (Kiminki, McSwain & Kobulnicky 2008; Kiminki et al. 2009, 2012; Kobulnicky et al. 2012). The remaining 37 spectral classifications (19 per cent) used in this work come from 12 different studies with between one and seven classifications each, including many studies devoted to a single particularly interesting star. Classifications from these focused studies were favoured where possible as they often uncovered previously unknown binary companions that changed the classification of the primary star.

Spectral classifications are presented in Table B1 with references, known binary companions, and source numbers under the systems initiated by Johnson & Morgan (1954, the Schulte system, e.g. Cyg OB2 #5), Massey & Thompson (1991, e.g. MT426), and Comerón et al. (2002, e.g. A26). We note that the star listed as ‘A11’ by Comerón et al. (2002) is actually MT267 from Massey & Thompson (1991), and is listed as the latter in this work. By comparing spectral classifications of the same star performed by different authors, we estimate the uncertainty in the spectral types presented here to be approximately half a subtype for O and early B-type (B0–B1) stars and one subtype for later B-type stars. The standard deviation in luminosity class is approximately half a luminosity class for O-type stars, one luminosity class for early B-type stars (B0–B2), and may be higher for later B-type stars, but there is insufficient data for a meaningful comparison. These uncertainties are used in all following calculations to derive uncertainties on other parameters.

2.2 Model atmosphere fits

In this work physical parameters for stars will mainly be based on spectral types from the literature, using the effective temperature scales of Martins et al. (2005), Trundle et al. (2007), and Humphreys & McElroy (1984). However, effective temperatures derived from model atmosphere fits should be more reliable, both because they are not quantized in spectral subtypes (or half subtypes), but also because they take into account other parameters such as stellar mass-loss and chemical abundances.

We searched the literature and identified 11 stars with suitable model atmosphere fits from the works of Herrero, Puls & Najarro (2002) and Negueruela et al. (2008). We did not include the model atmosphere fits from Herrero et al. (1999) since these used mass-loss rates and line blanketing analyses that have since been superseded (Herrero, private communication). We have also excluded stars that have since been found to be binaries with similar luminosities, as these will have blended spectra that will affect the model atmosphere fits.

The 11 stars with model atmosphere fits are listed in Table 1 with the fitted effective temperatures and gravities. Compared to the effective temperatures and gravities determined from the spectral types (also listed in Table 1) there is good agreement within the uncertainties. Given that both methods employ the same effective temperature scale, this is perhaps not surprising.

For the remainder of this work, for the 11 stars listed in Table 1, we used the effective temperatures derived from the model atmosphere fits and then calculated intrinsic colours and bolometric corrections (BC) for these stars as a function of effective temperature (instead of spectral type), as described below.

2.3 Photometric data

Photometric data in the optical and near-IR was gathered from the literature to complement the spectroscopy. In the optical the primary source of photometry was the UBV data from Massey & Thompson (1991), which is commonly used for studies of Cyg OB2. However, some of these photometric measurements are tagged by Massey & Thompson (1991) as ‘blended’, including 11 stars in our sample. We therefore supplemented the UBV data with Ug′r′ photometry from the Ultra-violet Excess Survey (UVEX, Groot et al. 2009) and converted this on to the UBV system (see Appendix A). For the 11 ‘blended’ sources, 7 had sufficient photometry to calculate full UBV photometry while 4 had only U and g′ UVEX photometry (r′ was saturated) and so only UB photometry was calculated. In addition, 12 sources in our sample were not from the original Massey & Thompson (1991) sample (predominantly from the near-IR sample gathered by Comerón et al. 2002), and we therefore sought UVEX photometry to homogenize our sample. Of these 12 sources, 9 had sufficient photometry to derive U and B photometry, while the remaining 3 were either saturated in the UVEX images or the sources were not observed. For all these sources the photometric uncertainties used were propagated from the UVEX photometric uncertainties and the dispersions on the transformations used.

Near-IR photometry was taken from 2MASS for the majority of stars. In view of the low spatial resolution of the 2MASS observations that could lead to blending in some cases, we sought replacement photometry from the UKIDSS-GPS (United Kingdom Infrared Deep Sky Survey Galactic Plane Survey; Lucas et al. 2008), but were only able to obtain unsaturated photometry for three sources. For these sources, we replaced their 2MASS photometry by UKIDSS photometry using the transformations in Lucas et al. (2008).

The final photometric data used in this study is presented in Table B1 on the UBV and 2MASS JHK_s photometric systems. The vast majority (93 per cent or 155 out of 167) of the primary stars have optical photometry in at least two of the UBV bands, and all of the targets have near-IR photometry in J, H, and K_s.

2.4 Binary stars and their influence on stellar parameters

Unresolved binary systems can lead to the miscalculation of the properties of the primary star. The contribution of light from a binary companion, which in involved systems is cooler than the primary, will make the primary appear cooler and more luminous, which could potentially lead to the star being diagnosed as more evolved and higher mass. The binary population of massive stars in Cyg OB2 has recently been well studied (e.g. Kiminki et al. 2007, 2009, 2012; Caballero-Nieves et al. 2014) and 22 of our 167 primary stars have secondaries with known spectral types, although it is likely that there still remain a number of undiscovered binary companions in this sample.

For objects with known companions, it is possible to correct for the light of the secondary when using the observed photometry as long as the spectral type of both stars is known. This is relatively simple when the luminosity ratio of the binary system is known in all the relevant bands, or if the bolometric luminosity ratio is known and the stars are of similar spectral type (e.g. Cyg OB2 #5 is an O7 I + O6 I binary with a bolometric luminosity ratio of 3.1; Linder et al. 2009). For other binary systems this is trickier and we must make an assumption about the relative in-band fluxes of the two stars, for which the simplest approach and that adopted here is to assume that the ratio of fluxes is equal to the ratio of intrinsic fluxes according to spectral type from Martins et al. (2005) for O-type stars and Conti, Crowther & Leitherer (2008) for B-type stars.

Doing this for our sample, we find that the typical (rms) adjustment to the colour of the primary star is very small, 0.006 mag, having a negligible effect on the extinctions. This is because the strongest effect occurs when two stars have significantly different intrinsic colours, but in such systems the stars usually have very different luminosities and therefore the influence of the secondary is small. The rms increase in the observed magnitude is more significant at 0.5 mag, equivalent to a 37 per cent reduction in the bolometric luminosity. This luminosity reduction (magnitude increase) varies from 8 per cent (0.09) for Cyg OB2 #11 (O5 If + B0 V) to 50 per cent (0.75) for very similar pairs of stars such as Cyg OB2 #73 (O8 III + O8 III) and B17, an O7 I + O9 I binary.

Cyg OB2 does contain one known triple system, MT429, which Kiminki et al. (2012) find to be a B0 V + B3 V + B6 V system. For simplicity we consider only the two brightest components of the system because the B6V star is too faint to affect the colour or luminosity of the primary star. We also note that while Cyg OB2 #22, an O3 If + O6 V binary (Walborn 2002), has been resolved in high angular resolution images with a separation of ∼1.5 arcsec between the two components (Mason et al. 2009; Maíz Apellániz 2010), these observations do not provide sufficient photometry to place both components on the HR diagram separately, therefore we only consider the primary star in this work (after applying the corrections for the secondary noted above).

2.5 Completeness

It is difficult to estimate the completeness of a sample that is compiled from so many different sources of spectroscopy, each of which is based on a different method of photometric selection, area of study, and spectroscopic completeness. However, one possible indicator of the level of completeness of our sample is the fact that the most recent spectroscopic studies (e.g. Kiminki et al. 2007; Kobulnicky et al. 2012) have discovered only B-type stars and not new O-type stars, possibly suggesting that the majority of O-type stars in Cyg OB2 have already been detected. The fact that the most recent studies have found stars with similar extinctions to previous studies (Comerón & Pasquali 2012) also suggests that a previously undiscovered population of highly obscured (A_V > 10 mag) O-type stars does not exist. Based upon this, we will make the assumption that our census of massive stars is complete for O-type stars in Cyg OB2, which is broadly equivalent to being complete down to 15 M_⊙ for zero-age main-sequence (MS) stars, or 20 M_⊙ if one takes into account some further evolution. In a future paper, we will attempt to test this assumption about our completeness using new spectroscopic observations of Cyg OB2.

3 THE EXTINCTION TOWARDS CYGNUS OB2

In this section, the spectral types and photometry are used to rederive the form of the extinction law towards Cyg OB2 and then calculate individual extinctions for all the stars in our sample. To achieve this, we use the unreddened colours from Martins & Plez (2006) for O-type stars (the ‘observational’ T_eff scale) and from Fitzgerald (1970) and Koornneef (1983) for B-type stars, adjusted for B0–B1 stars to provide a smooth transition with the O star colours.

3.1 The extinction law towards Cygnus OB2

To determine the form of the extinction law towards Cyg OB2, we selected a high-confidence subset of 18 stars from within our sample to which we perform spectral energy distribution (SED) fitting. We selected only single O-type MS stars classified by Massey & Thompson (1991) and with no evidence of having a similar-luminosity binary companion. By using only O-type stars, we reduce the uncertainty due to the effective temperature and the binarity of the stars, both of which increase towards later spectral types. In addition we required all stars to have been detected in, and have low errors (<0.05 in UBV and <0.02 in JHK_s) in, all six photometric bands. These criteria resulted in the sample of 18 stars in Table 2.

Reddened stellar spectra were calculated using CMFGEN non-LTE model spectra (Hillier & Miller 1998), as used by Martins et al. (2005), reddened using the extinction curves presented by Fitzpatrick & Massa (2007) as a function of R_V and A_V. Reddening coefficients were determined using Kitt Peak UBV filter profiles (appropriate for the observations of Massey & Thompson 1991) and 2MASS JHK_s filter profiles (Skrutskie et al. 2006). Combining these with the intrinsic colours for O-type dwarf stars presented by Martins & Plez (2006) results in reddened photometry as a function of spectral type, R_V, and A_V.

χ² fits were performed between the observed photometry and the reddened model photometry allowing R_V to vary in steps of 0.01 in the range 2.5–3.5 and with A_V varying between 4.0–7.0 mag in steps of 0.01 mag. The distance modulus was left as a free parameter to account for variations in the intrinsic distance and luminosity of stars. The best fits are shown in Fig. 2 and listed in Table 2.

Figure 2.

R_V versus A_V for our high-confidence sample of 18 single O-type dwarfs derived using χ² fits between the observed photometry and reddened model spectra using the Fitzpatrick & Massa (2007) reddening curves. 1σ error bars are shown derived from the χ² fits. The weighted mean R_V = 2.91 ± 0.06 is shown as a dashed line.

Open in new tabDownload slide

The weighted mean extinction law from our χ² fits is R_V = 2.91 ± 0.06 (taken as such for the remainder of this work), in good agreement with previous studies based on optical colours (e.g. Hanson 2003, measured R_V = 3.0). However, this value is in disagreement with the findings of Comerón & Pasquali (2012) who use a combination of optical and near-IR photometry to derive R_V = 4.0 using a Cardelli, Clayton & Mathis (1989) extinction law for OB stars over a wider area across the Cygnus X molecular cloud. This disagreement is likely explained by there being two different extinction laws at work along different sight lines: a ‘standard’ R_V ∼ 3 extinction law towards Cyg OB2 (due to obscuration from the Cygnus Rift, which lies predominantly in the foreground; Guarcello et al. 2012), and an ‘anomalous’ R_V ∼ 4 extinction law away from Cyg OB2, possible due to a change in the grain size distribution favouring larger dust grains within the Cygnus X molecular cloud (which has already been dispersed from the vicinity of Cyg OB2; see Fig. 1). A similar picture of two extinction laws has been observed by Povich et al. (2011) towards the Carina molecular cloud.

For reference, this new extinction law is presented in Table 3 for extinctions between A_V = 4.0–7.0 (the typical range observed in Cyg OB2; see Section 3.2) and compared to that found previously by Hanson (2003) from comparison of UBVJHK_s colour excesses compared to the Lejeune & Schaerer (2001) intrinsic colours. The extinction law is relatively unchanged in the optical, but has changed significantly in the near-IR, due mostly to shifts in the near-IR intrinsic colours of O-type stars between the works of Lejeune & Schaerer (2001) and Martins & Plez (2006).

3.2 The extinction of Cygnus OB2 members

Using the Fitzpatrick & Massa (2007) R_V = 2.91 extinction law derived above, we calculated individual extinctions for the 164 primary O- and B-type stars in our sample (not including the WR stars) from the comparison between observed and intrinsic colours. We favoured the B − K_s colour for this due to its long baseline and the availability of the relevant photometry for the vast majority (160 or 98 per cent) of our sources. For one star B-band photometry was unavailable, so we used the V − K_s colour, and for three stars only near-IR photometry was available, so we used the J − K_s colour. Uncertainties in the derived extinction were calculated using a Monte Carlo (MC) simulation taking into account the uncertainties in the photometry, spectral classification, the intrinsic colours (assumed to be half the difference between the colours of neighbouring subtypes), and the reddening law of R_V = 2.91 ± 0.06. The typical uncertainty in the extinction is ∼0.3 mag in the V band. Extinctions and uncertainties are presented in Table B1.

Fig. 3 shows the extinction distribution, which varies from A_V = 2.2 mag for MT140 (an O9.5 I star on the western side of the association) to A_V = 10.2 mag for the central blue hypergiant Cyg OB2 #12 (a B3.5 Ie star, and one of the most luminous known stars in the Milky Way; Morgan et al. 1954; Clark et al. 2012), though the majority lie in the range A_V = 4–7 mag with a median of A_V = 5.4 mag and 25 per cent and 75 per cent quartiles of 4.7 and 6.5 mag. Error bars on the extinction distribution histogram were derived from the individual uncertainties using an MC simulation. The large spread in extinction compared to the typical uncertainty in A_V suggests that the spread is real and is either due to foreground extinction from the Cygnus Rift or intra-association extinction within Cyg OB2.

Figure 3.

Histogram of the extinctions calculated for all 164 massive OB stars in this work with 1σ error bars derived from MC simulations.

Open in new tabDownload slide

Fig. 4 shows a spatial map of the extinction across Cyg OB2 in the form of a Voronoi tessellation. The extinction distribution varies smoothly across the association, from A_V ∼ 4–5 mag in the north-west (where a known ‘reddening hole’ exists; Reddish, Lawrence & Pratt 1967) to A_V ∼ 6–7 to the south of the association and on the outskirts.

Figure 4.

Voronoi tessellated extinction map (one star per Voronoi cell) for all 164 massive stars in Cyg OB2 with the extinction indicated by the grey-scale colour. The circular border of our study is approximately at the border of this image.

Open in new tabDownload slide

These results are in good agreement with previous studies of both the high- (Hanson 2003) and low-mass populations (Wright et al. 2010b), albeit with a smaller spread than the latter, likely due to the higher precision afforded by having known spectral types and a calibrated extinction law. The extinction is larger than that derived by Guarcello et al. (2012) from their riz study by ∼1 mag (they derive a median extinction of A_V = 4.3 mag), though we observe the same spatial variation in extinction that they do. The anomalously low extinction derived by Guarcello et al. (2012), particularly compared to many other studies, may be due to the difficulties deriving extinction from isochrone fits in the riz filter system that requires a complex transformation from BVI colours.

The high extinction for Cyg OB2 #12, >2 mag larger than any of the stars surrounding it, is particularly apparent. Some of this extinction may be due to circumstellar material, given the advanced evolutionary stage of this object and that material is known to surround other very massive stars (e.g. Wright et al. 2014b); however, Clark et al. (2012) could find no evidence for a near- to mid-IR excess in the SED of this object, arguing against a significant amount of warm circumstellar dust. It is also possible that the SED we have assumed is wrong, given the unusual character of the star.

4 HR DIAGRAM FOR Cyg OB2

In this section, we place the stars on the HR diagram and compare their positions with both rotating and non-rotating stellar evolutionary tracks to derive stellar masses and ages. Cyg OB2 also includes four evolved massive stars (such as Wolf–Rayet stars) whose positions on the HR diagram and resulting stellar properties cannot be derived in the same way as the less-evolved OB stars. We have gathered the properties of these stars from the literature, adjusting for a distance of 1.45 kpc if necessary, and derive additional physical properties (such as upper or lower limits on stellar ages of masses) by comparison with the evolutionary models used in this work or from our current, but far from complete, understanding of the evolution of massive stars. These stars and their properties are included in Table B1 and included in all later analysis where possible.

4.1 Calculation of L_bol and T_eff

Bolometric luminosities were calculated using the observed K_s band magnitudes because of the high extinction towards Cyg OB2 and the typical uncertainty of 0.3 mag in A_V resulting from this. The considerably lower uncertainty in A_K compared to A_V is preferable despite the higher uncertainty in BC_K over BC_V (due to the strong sensitivity of the SED to wind parameters at longer wavelengths; Martins & Plez 2006). BC were taken from Martins & Plez (2006) for O-type stars, Crowther, Lennon & Walborn (2006) for B-type supergiants, and Humphreys & McElroy (1984) for less luminous B-type stars, although there is very little variation with luminosity class. Effective temperatures were determined as a function of spectral type using the tabulations of Martins et al. (2005, ‘observed’ scale) for O-type stars, Trundle et al. (2007) for early (B0–B3) B-type stars and Humphreys & McElroy (1984) for later B-type stars, all of which show a good agreement between each other.

The accepted distance to Cyg OB2 has been reduced over the last few decades due to improved observations and calibration data. Reddish, Lawrence & Pratt (1966) originally placed the association at a distance of 2.1 kpc, while Massey & Thompson (1991) estimated a nearer distance of ∼1.75 kpc [distance modulus (DM) = 11.20] following their spectroscopic study of massive stars in the association. The recent recalibration of OB star absolute magnitudes led Hanson (2003) to determine an average distance of 1.2 kpc (DM = 10.44) but concluded based on inconsistencies with current stellar wind theory that this was too low and suggested a value of 1.45 kpc (DM = 10.80). This is the value we use in this work, which is in agreement with the parallax distance of 1.40 ± 0.08 kpc for the Cygnus X complex (Rygl et al. 2012).

The uncertainties in our values of T_eff and L_bol were derived from MC simulations. For the former we use the quoted uncertainties in T_eff from Martins et al. (2005), Trundle et al. (2007), and Humphreys & McElroy (1984), while for the latter we take into account the uncertainties in the observed photometry, derived extinctions, BC (σ = 0.1 mag for O-type stars in the K band, Martins & Plez 2006, and assumed to be σ = 0.2 mag for B-type stars), and spectral classification. For the distance to Cyg OB2, we assume that the distance of 1.45 kpc is correct (small variations in this value up to a few hundred parsecs do not significantly affect the results of this paper, though larger values do).

4.2 The HR diagram

The new HR diagram for Cyg OB2 showing 167 stars is shown in Fig. 5 compared to the isochrones and evolutionary tracks from the non-rotating models of Lejeune & Schaerer (2001) and in Fig. 6 compared to the rotating models of Ekström et al. (2012). The non-rotating models at a metallicity of Z = 0.02 are shown partly for comparison with previous studies (e.g. Hanson 2003; Comerón & Pasquali 2012), and also to contrast with the more recent, rotating models at the revised solar metallicity of Z = 0.014. An indication of the uncertainties on L_bol and T_eff are shown in Figs 5 and 6 for each spectral subclass (where possible) and in Fig. 6 for individual stars to illustrate the spread they induce in the HR diagram. For O-type stars the uncertainties are typically small, ∼0.02 dex in log T_eff and ∼0.1 dex in log L/L_⊙, though for the B-type stars the errors rise significantly to ∼0.07 dex in log T_eff and ∼0.2 dex in log L/L_⊙ due to classification uncertainties and the larger difference in log T_eff between spectral subtypes.

Figure 5.

HR diagram for Cyg OB2 showing all the stars compiled in this work colour-coded by their luminosity class (dwarfs are black, subgiants are blue, giants are green, bright giants are orange, and supergiants are red, while WR stars are shown with a star symbol) and with source numbers noted according to the scheme adopted in this paper. The figure shows isochrones and evolutionary tracks from Lejeune & Schaerer (2001) with ages and masses labelled. The effective temperature scales used for MS O- and B-type stars are shown along the top of the diagram and mean 1σ errors bars for each subtype are shown to the bottom left.

Open in new tabDownload slide

Figure 6.

HR diagram for Cyg OB2 showing all the stars compiled in this work as per Fig. 5, with the difference that the main panel shows isochrones and evolutionary tracks for rotating (v_ini/v_crit = 0.4) models from Ekström et al. (2012, smoothed for clarity) with the ages and masses labelled. The bottom-left figure shows the error bars for each source and the isochrones and evolutionary tracks. The bottom-right figure shows the distribution of all stars without labels and evolutionary tracks.

Open in new tabDownload slide

As noted by previous authors (e.g. Massey & Thompson 1991; Hanson 2003), there is a well-populated MS extending up to O5.5V,¹ for which there is now the data to extend it down to mid-B spectral types.

4.3 A large luminosity spread amongst the B stars

Many evolved massive stars are evident, both at masses M > 25 M_⊙ where such evolution can occur within 2–10 Myr, and at lower masses where evolved B-type stars occupy positions in the HR diagram suggesting ages >10 Myr. Hanson (2003) noted in her sample a large luminosity spread in the B1–B2 dwarfs and this is also apparent here, although it is now more concentrated towards the MS (most likely because the massive stars discovered recently have been fainter as the depth of spectroscopic observations has increased). Hanson (2003) and Negueruela et al. (2008) have suggested that their position in the HR diagram may be due to them being foreground contaminants from a separate population. However, we note that for their luminosities to be consistent with foreground, MS B-type stars would require some of the most overluminous stars to lie at distances as low as ∼500 pc, which would be inconsistent with their high extinctions (the five most overluminous B-type stars all have A_V > 5.5 mag) and the estimated distance to the extinguishing Cygnus rift of 800–1400 pc (Guarcello et al. 2012; Sale et al. 2014).

Alternatively, the apparently large luminosity spread amongst the B-type stars could be due to the large uncertainties in L_bol and T_eff at these spectral types that are smearing out the MS in this range. Observational uncertainties in spectral type combined with the steep relationship between T_eff and spectral subtype in this range creates large uncertainties in T_eff amongst the early B-type stars. Comparing the observed spread in the HR diagram with our estimated uncertainties suggests there is no significant evidence favouring a foreground population of B-type stars or a cospatial population of mature (>10 Myr) stars to explain the observed spread. This argument is supported by the fact that the seemingly high-luminosity B-type stars have a similar spatial distribution to the low-luminosity B-type stars, suggesting they are unlikely to be a separate population.

4.4 Possible foreground and background contaminants

A small number of stars fall significantly below the MS in the HR diagram such as MT170 (B5 V, log L/L_⊙ = 1.74), MT426 (B0 V, log L/L_⊙ = 2.51), and MT427 (B4 III, log L/L_⊙ = 2.19). MT170 and MT427 appear to be background sources, most likely at distances of ∼ 2.19–3.5 kpc if they are on the MS (a number of Hii regions and massive young stellar objects have been identified at similar distances by Xu et al. 2013). Their moderate extinctions (A_V = 3.6 and 3.8 mag, respectively) are acceptable given that they are projected against the ‘reddening hole’ in the north-west of the association (see also Drew et al. 2008, who find evidence of a population for A-type stars at such distances in this direction).

The low luminosity of MT426 is harder to explain as it would imply that the star is ∼10 kpc distant yet has an extinction of only A_V ∼ 4 whilst being projected against the centre of the association where the extinction is typically A_V = 6–7 mag (Fig. 4). MT426 is listed by Massey & Thompson (1991) as having ‘blended’ photometry (see Appendix A) with V = 14.05 mag and B − V = 1.95 mag. These values have led previous authors to derive reasonable extinctions and luminosities for this star, yet such values would imply a K_s-band magnitude of ∼8.4, yet the source is not detected by 2MASS. Inspection of the UVEX and UKIDSS images reveals a second point source that is highly blended in the lower spatial resolution Massey & Thompson (1991) and 2MASS images with its neighbour MT425. The UVEX and UKIDSS photometry, which we have used, does not appear discrepant and we believe it to be sound. It is possible that the spectral type of B0 V is erroneous, given the blending reported by Massey & Thompson (1991) and the fact that those authors report an identical spectral type for its neighbour (though Kiminki et al. 2007, report the same spectral type, although it may also have suffered blending).

Due to the unresolved nature of MT426 and the fact that MT170 and MT427 appear to be background sources, we argue that all three are non-members of Cyg OB2 and exclude them from further analysis of the association. They are noted in Table B1 as ‘Unresolved’ or ‘Background’ sources, and physical properties are not derived for them.

Finally, we note that the membership status of Cyg OB2 #12 has been questioned in the past because of its advanced evolutionary state compared to the early O-type stars in the association and its high luminosity (Walborn 1973), with the suggestion that the star might be a foreground, though it could also be an luminous blue variable (Clark et al. 2012) or a similar-luminosity binary (Caballero-Nieves et al. 2014 recently identified a companion to the star, but note that it is too faint to affect the calculated properties of the star).

4.5 Physical quantities derived from the HR Diagram

Stellar ages and initial masses were calculated based on the position of each star in the HR diagram relative to the non-rotating Z = 0.02 models of Lejeune & Schaerer (2001, commonly used for previous studies of Cyg OB2) and the rotating Z = 0.014 (the revised solar metallicity value; Asplund, Grevesse & Sauval 2005) models of Ekström et al. (2012), both from the Geneva group. The rotating models have an initial rotation rate of v_ini/v_crit = 0.4, which corresponds to typical MS velocities of 110–220 km s⁻¹ for massive stars (Ekström et al. 2012), in agreement with observations (e.g. Dufton et al. 2006). However, it is not currently clear whether stars are born with a narrow or broad range of initial rotation velocities.

Stellar mass-loss is an important ingredient in stellar models, highly influenced by the metallicity (Meynet et al. 1994), rotation rate (Maeder & Meynet 2000), and the exact mass-loss prescription used (e.g. Vink, de Koter & Lamers 2001). The non-rotating evolutionary models of Lejeune & Schaerer (2001) used the radiative mass-loss rates from de Jager, Nieuwenhuijzen & van der Hucht (1988) for stars in the blue part of the HR diagram, while Ekström et al. (2012) update this with the prescription from Vink et al. (2001) where possible, and also take into account supra-Eddington and mechanical mass-loss processes in their models. Different mass-loss rates can play a significant role in the evolution of massive stars, their positions in the HR diagram, and therefore the physical parameters derived in such a way. These differences are most pronounced for red supergiants and the most massive stars (40–120 M_⊙; Ekström et al. 2012), with the latter having considerably increased MS lifetimes when mass-loss rates are increased (Meynet et al. 1994). For moderately massive stars (20–40 M_⊙), the effects of mass-loss can be less notable (e.g. Meynet et al. 1994 note that the MS lifetime of 25 M_⊙ stars increased by only 5 per cent when the mass-loss rate was doubled).

The ages and masses derived from these models are listed in Table B1 with uncertainties. For the vast majority of objects this is a simple process and we do this on the assumption that all the stars are on their first journey across the HR diagram from hotter to cooler temperatures (generally a valid assumption since stars move back to hotter temperatures much faster and therefore it is rarer to observe stars during this transition). For stars that fall to the left of the zero-age MS in the HR diagram, we assume that their positions are due to uncertainties in their luminosity or effective temperature and that these are zero-age MS objects.

To derive uncertainties on the age and initial mass of stars, we run an MC simulation that takes into account both the observational uncertainties on L_bol and T_eff as well as uncertainties due to unresolved binary companions and an intrinsic spread in both luminosity and temperature for a star of a given initial mass and age. The latter simulates the spread in the HR diagram at a given initial mass and age caused by variations in fundamental parameters such as rotation and metallicity, and was estimated to be 0.12 dex in log L/L_⊙ and 0.02 dex in log T_eff from inspection of the spread in the upper MS in rotating stellar models, and an extrapolation of the spread measured by Houk et al. (1997) in Hipparcos data. The effects of unresolved secondaries were modelled in our MC simulation using the binary characteristics found by Kiminki & Kobulnicky (2012) for Cyg OB2: a binary fraction of 90 per cent for O-type stars (we assume a slightly lower binary fraction of 80 per cent for the B-type stars) and a flat mass ratio distribution over the range q = 0.005–1.0. The added luminosity uncertainty due to unresolved binary companions is generally small, typically adding an uncertainty of 0.1 dex, and therefore we do not consider more complex multiple systems. The uncertainties shown in Figs 5 and 6 (for both individual stars and representative values for each spectral subtype) are the uncertainties on the calculated values of L_bol and T_eff only and do not take into account these additional uncertainties used when calculating the stellar mass and age. The median fractional uncertainty on the stellar mass is found to be σ_M/M ∼ 0.14 and for stellar age is σ_τ/τ ∼ 0.3 (for stars in the range 0–10 Myr).

5 THE AGE OF Cyg OB2

Using the stellar ages derived from the HR diagram, we can assess the age distribution of stars in Cyg OB2 and hence the star formation history of the region. This is shown in Fig. 7 using both the non-rotating and rotating evolutionary models to highlight the different age distribution that they each imply. We only show the 47 stars more massive than 20 M_⊙ (down to a MS spectral type of ∼B0.5 V) based on our estimated completeness level and because the positions of stars less massive than this in the HR diagram do not provide an accurate indication of their ages (due to the closeness of the isochrones and the large uncertainty on T_eff for these stars). This excludes all the apparently overluminous B-type stars discussed earlier. The uncertainties shown in Fig. 7 were calculated from an MC simulation of the data analysis and binning process.

Figure 7.

Histograms showing the distribution of derived ages for stars with reliably determined ages from the non-rotating (left) and rotating (right) models. The ages of all stars with M ≥ 20 M_⊙ are shown with a black line. In the left-hand panel, we also show in red the age distribution of stars with 20 ≤ M/M_⊙ ≤ 40, which should be complete and unbiased in the age range of 0–4.5 Myr. In the right-hand panel, we show in red the age distribution of stars with 20 ≤ M/M_⊙ ≤ 35 that are unbiased in the age range of 0–6 Myr. 1σ uncertainties are shown for all histograms derived from an MC simulation of the data analysis and binning process.

Open in new tabDownload slide

5.1 Comparison of age distributions from rotating and non-rotating models

The age distribution of massive stars derived from the non-rotating models shows a strong peak at 2–3 Myr with the majority of stars in the age range of 0–5 Myr, in good agreement with previous studies of the massive star population of Cyg OB2 (e.g. Hanson 2003; Comerón & Pasquali 2012). However, since a 20 M_⊙ star has a lifetime of only ∼8 Myr (in the Lejeune & Schaerer 2001 models) then our sample of massive stars with masses >20 M_⊙ is only sensitive to stars <8 Myr old. We therefore cannot reliably infer the star formation history greater than 8 Myr ago when using these models. Furthermore this distribution is biased by the fact that the most massive stars have total lifetimes <4 Myr. For example a 60 M_⊙ star is estimated to have a total lifespan of ∼3.5 Myr, and therefore any such stars born more than 4 Myr ago will not be seen now. To rectify this well-known bias, we also show in Fig. 7 the distribution of ages of stars in the mass range 20 ≤ M/M_⊙ ≤ 40, which should be unbiased in the age range of 0–4.5 Myr. This presents a slightly different picture with a smaller peak suggesting that this peak age is an observational bias caused by the short lifetimes of the most massive members of the association.

The rotating stellar models of Ekström et al. (2012, Fig. 7) show a different picture with a prominent but wide peak at an age of 4–5 Myr. The increase in the typical age of massive stars in Cyg OB2 when using the rotating stellar models can be attributed to the changes in evolutionary time-scales induced by rotation, which increases the MS lifetime by ∼25 per cent (Ekström et al. 2012). The marked influence of this on stellar ages is evident in the two HR diagrams (Figs 5 and 6). For example, the group of O4–7 supergiants with L/L_⊙ ∼ 5 × 10⁵ L_⊙ have ages of 1–3 Myr according to the non-rotating models, but have ages of 3–5 Myr from the rotating models. The older ages are in better agreement with those derived from lower mass stars in Cyg OB2, e.g. Drew et al. (2008) estimate an age of 5–7 Myr from A-type stars, while Wright et al. (2010b) derive an age of 3.5–5.25 Myr from an X-ray selected sample of approximately solar-mass stars.

The increased stellar lifetimes of rotating stellar models means that the mass range over which our sample is unbiased has also changed relative to the non-rotating models. Fig. 7 shows the age distribution of stars in the mass range 20 ≤ M/M_⊙ ≤ 35, which should be complete in the age range 0–6 Myr. This distribution is remarkably flat over the unbiased age range, suggesting a relatively constant level of star formation over the last 6 Myr. The peak stellar age previously observed at 4–5 Myr is mostly due to stars in the HR diagram (Fig. 6) at masses of 40–60 M_⊙ (the O4–7 supergiants noted earlier) that have maximum stellar ages of <5 Myr. The large number of slightly lower mass (25–50 M_⊙) stars with ages of 5–6 Myr (the O7–9 supergiants observed in the HR diagram) could suggest that many of the more massive members have already evolved to their end states and their loss biases the observed age distribution.

The comparisons with both sets of evolutionary models indicate a number of stars in Cyg OB2 with very young ages, ≤1 Myr. The existence of such young stars is not impossible but appears unlikely given that long wavelength images (e.g. Fig. 1, and Schneider et al. 2006) show that the majority of dust and molecular gas been evacuated from the cavity in which Cyg OB2 resides. Furthermore both Wright et al. (2010b) and Guarcello et al. (2013) found a paucity of Class I or heavily embedded low-mass stars within the association, suggesting the low-mass stellar population is not so young (though this could also be due to the high levels of feedback in the vicinity of the association; Wright et al. 2012; Guarcello et al. 2014). The typical uncertainties in the HR diagram shown in Figs 5 and 6 suggests that the positions of these stars are not incompatible with them being slightly older (>1–2 Myr) and the lack of very massive stars (>30 M_⊙) with ages <1 Myr (where the large separation of isochrones allows young ages to be estimated more accurately) supports this idea.

To conclude, our unbiased age distribution derived from rotating stellar models suggests a relatively constant level of star formation over the past ∼6 Myr (Fig. 7) with some evidence for a peak in the star formation intensity between 4–5 Myr ago. We cannot be certain from this sample when star formation in the region began, but there is some evidence that it started at least ∼6–8 Myr ago. Addressing the question of when star formation ceased is harder, but our best estimate is that it ended ∼1 Myr ago. As will be discussed in Section 6.1, these results are in approximate agreement with our star formation history fits to the mass function of stars in Cyg OB2. The use of different mass-loss prescriptions is not expected to significantly affect the evolution of stars in our ‘unbiased’ mass range of 20–35 M_⊙ (Meynet et al. 1994), and therefore this is unlikely to explain the age spread observed. It is possible that the large spread in the HR diagram could instead be due to a single-aged population with an extremely large range of rotation velocities, which has been shown to lead to large spreads in the HR diagram (Brott et al. 2011) and could reproduce the entire spread seen in our HR diagram. However, for this to be the case the stars nearest the zero-aged MS in the HR diagram would have to be rotating very rapidly with speeds of 400–450 km s⁻¹, including Cyg OB2 #20 and #29, for which Herrero et al. (1999) measured V sin i = 25 and 180 km s⁻¹, respectively. This would therefore require Cyg OB2 #20 to have been observed virtually pole-on, which makes this scenario unlikely, though further high-resolution spectroscopy of other stars in Cyg OB2 would be useful to test this hypothesis in more detail. It is also worth noting though that this scenario is not supported by measurements of rotation rates in other clusters (e.g. Dufton et al. 2006).

The uncertainties in the stellar ages we have calculated could explain some of this spread, but do not appear to be sufficiently large enough to explain the entire ∼5 Myr spread. Furthermore, the age distribution of our ‘unbiased’ sample shown in Fig. 7 does not resemble a single-aged population broadened by observational uncertainties (which should have a more centrally-concentrated peak) and does appear to represent a spread of ages amongst the OB star population.

5.2 Implications for large age spreads in OB associations

Large age spreads have been reported in many star clusters and associations, typically diagnosed from luminosity spreads in colour–magnitude diagrams for low-mass pre-MS stars. Some authors have suggested that such spreads in luminosity do not equate to age spreads and may instead be due to intrinsic luminosity spreads amongst a coeval population (e.g. Hillenbrand, Bauermeister & White 2008; Jeffries et al. 2011) potentially caused by different accretion histories (e.g. Baraffe, Chabrier & Gallardo 2009). Age spreads diagnosed from post-MS massive stars are not believed to suffer from such issues since the intrinsic spread and measurement uncertainties in L_bol and T_eff are typically smaller than the separation of isochrones (see also Massey 2003). Therefore if the HR diagram is sufficiently well sampled to overcome statistical and observational uncertainties (as it is here), age spreads derived from high-mass post-MS stellar evolution models can provide more accurate age diagnostics than low- and intermediate-mass pre-MS stellar evolution models.

The reported age spreads for low-mass stars are often argued as representing a continuous period of star formation within a single cluster or association and in a small area of space (e.g. Palla & Stahler 2000). Cyg OB2 however occupies a much larger region than typical young compact star clusters such as the Orion Nebula Cluster, with the stars spread over >10 pc. It has recently been argued that Cyg OB2 has always been a large and low-density OB association and was never a dense and compact star cluster (Wright et al. 2014c, see also Parker et al. 2014)). It is therefore possible that the observed age spread is due to either an age gradient or a series of discrete star-forming events that have since mixed. This would imply that Cyg OB2 was not born as a single star cluster but as a ‘distribution’ of smaller groups or clusters of stars with a range of stellar ages. We searched for evidence of a spatial variation in the stellar age distribution to support this theory, but could not find any such evidence. If the association is composed of multiple populations of different ages then they are well mixed.

6 THE UPPER END OF THE IMF IN Cyg OB2

In this section we analyse the IMF of the massive stars in Cyg OB2, as derived using rotating stellar models (because Weidner & Vink 2010 find that such masses are in much better agreement with dynamical masses than those derived from non-rotating stellar models). Fig. 8 shows a histogram of the masses of all 167 stars in our sample, showing a rapid increase as one goes from high to low masses, with the data flattening off for stars less massive than our estimated completeness limit of ∼20 M_⊙. To calculate the power-law slope of the mass function, we fitted these data using an ordinary least-squares regression of Y on X (OLS (Y|X); Isobe et al. 1990) of the form dN/d log m ∝ m^−Γ, and for the mass bins >20 M_⊙. Because the uncertainty on this fit is just a fitting uncertainty, we also performed an MC simulation to determine the true uncertainty, varying the masses of the individual stars according to their uncertainties and adjusting the number of stars in each bin as appropriate. The 1σ spread in the best-fitting slope that results is then combined with the fitting error to produce the true uncertainty on the mass function slope. The best-fitting mass function slope for all stars in our sample is Γ = 2.06 ± 0.06, significantly steeper than the canonical Salpeter slope of Γ ∼ 1.35 (Salpeter 1955), although it agrees with that determined by other authors for the massive stars in Cyg OB2 (e.g. Kiminki et al. 2007; Comerón & Pasquali 2012). Note however that this uncertainty is still falsely low because it does not take into account the random sampling of the IMF.

Figure 8.

Mass functions for massive stars in Cyg OB2 derived using the rotating stellar models of Ekström et al. (2012). Left: histograms showing the distribution of derived masses for all 167 stars in our sample (black histogram, with 1σ error bars determined from MC simulations) and for an age-unbiased sample of 82 stars aged <3.5 Myr (red histogram, with 1σ error bars). Both histograms have been fitted using an ordinary least-squares regression of Y on X (OLS (Y|X), Isobe et al. 1990) for the mass bins >20 M_⊙, giving slopes of Γ = 2.06 ± 0.06 and 1.39 ± 0.19 for the complete and age-unbiased samples, respectively (note that these uncertainties do not take into account random sampling of the IMF). Right: a cumulative mass function showing N (M_⋆ < M) as a function of M for all 98 stars in our sample more massive than 10 M_⊙ (black line) with the 1σ uncertainty (calculated from MC simulations) shown in grey. Also shown are the results of simulated IMFs for Γ = 1.3 (red line, with 1σ random sampling uncertainty shown in pink) and Γ = 2.2 (green line) calculated using the equations of Maschberger (2013, see full explanation in the text). The blue line shows a simulated and evolved IMF for Γ = 1.3 for stars aged 3–7 Myr (i.e. a constant star formation rate from 3 to 7 Myr ago, and no star formation since then). The uncertainties for three simulated IMFs are similar and so only one is shown for clarity.

Open in new tabDownload slide

The mass distribution of stars in our sample is likely to be altered by the loss of the most massive stars that have already evolved to their end states (Section 5), steepening the observed mass function. To counter this, we have attempted to identify an unbiased sample of stars from which to study the IMF in Cyg OB2, selecting all stars with ages <3.5 Myr, which we estimate from the rotating stellar models (Fig. 6) as the age at which the most massive stars in Cyg OB2 will have evolved to their end states. The mass function of this unbiased sample of 82 stars is shown in Fig. 8, which shows a slightly shallower slope to that of the full sample. An OLS (Y|X) fit gives a slope of Γ ∼ 1.39 ± 0.19, in good agreement with the concept of a universal IMF at high masses with a slope of Γ = 1.3 (Kroupa, Aarseth & Hurley 2001; Bastian et al. 2010). This value also agrees with the mass function measured at lower masses (Γ = 1.09 ± 0.13; Wright et al. 2010b), and the mass function of massive stars found by Massey & Thompson (1991) of Γ = 1.0 ± 0.3, although their sample was considerably smaller than the current sample and may therefore have been incomplete.

6.1 The cumulative mass function

To avoid the well-known and negative effects of binning our data, and to overcome any uncertainties arising from the narrow age range used for the unbiased mass sample, we have also considered the mass function of stars in Cyg OB2 represented as a cumulative function (CF). This is shown in Fig. 8 with 1σ uncertainties calculated from an MC simulation using the full uncertainties on all stellar masses (the 1σ uncertainties shown are very small due to the cumulative nature of the function and the lack of binning).

The CF shows a steady increase up to ∼30 M_⊙, and then flattens off, reaching a total of 98 stars with masses >10 M_⊙. We compare our CF with simulated IMFs generated using the equations described in Maschberger (2013), varying α (=Γ + 1), but fixing β = 1.4 (the low-mass exponent has no effect on the high-mass IMF, but is required for the calculations) and drawing individual stellar masses from the mass range 0.01–150 M_⊙. We continue to do this until the number of stars with M > 20 M_⊙ reaches the observed number of 98 and use this to produce a CF scaled to the observed number of stars with M < 20 M_⊙. The simulated IMF is generated 10 000 times and the median CF shown in Fig. 8 with ±1σ range illustrated.

We initially simulated an unevolved (i.e. counting all stars) ‘universal’ IMF (Kroupa et al. 2001) with Γ = 1.3, but found that this overpredicted the number of very massive stars (>50 M_⊙) relative to the moderately massive stars (20–50 M_⊙) and we can reject a standard, unevolved mass function of this type with a confidence of ∼3σ. Varying Γ we found a best fit with an unevolved IMF for |$\Gamma = 2.2^{+0.6}_{-0.3}$|⁠, in approximate agreement with that found fitting a binned IMF for our entire sample (but with larger uncertainties that also take into account the random sampling of the IMF).

We then simulated a number of evolved mass functions using the ‘universal’ Γ = 1.3 IMF and different star formation histories. For a star of a given mass and age, we compute its position in the HR diagram according to the models of Ekström et al. (2012) and remove those stars that have evolved beyond their reddest position in the HR diagram. This excludes WR stars from the simulated mass functions, matching our observed CF that also leaves them out (because of the difficulty of determining their initial masses).

We found that a number of different star formation histories provided satisfactory fits (within 1σ) to the observed mass function, preventing us from accurately determining it in this way, though we were able to place some partial constraints on it. We find fits using a ‘universal’ Γ = 1.3 IMF requires the star formation rate over the last ∼3 Myr to be very low, <20 per cent of the total star formation in Cyg OB2, with the fit improving slightly as the fraction of star formation in this time period decreases. Exploring models with this feature, we find acceptable fits when the duration of star formation is at least 2 Myr (i.e. a constant level of star formation from 3 to 5 Myr ago), with the best fit obtained for a model with a constant star formation rate over the time period of 3–7 Myr ago, in reasonable agreement with that found from the distribution of stellar ages studied earlier. Significant star formation prior to ∼7 Myr ago degrades the fit to the cumulative mass function. These results are in good agreement with the stellar age distribution discussed in in Section 5.

To conclude, the current mass function of massive stars in Cyg OB2 is relatively steep, with a slope of Γ ∼ 2.1, but the stellar age distribution suggests that it has been steepened by the most massive stars evolving to their end states. From both an unbiased sample of stars and modelling an evolved cumulative mass function, we find that the mass function is consistent with the canonical Γ = 1.3 ‘universal’ IMF (Salpeter 1955; Kroupa et al. 2001).

6.2 Have there been any supernovae in Cyg OB2?

There is considerable evidence, both in the age distribution of stars and in the steepened IMF, to suggest that some of the most massive stars in Cyg OB2 have already evolved to their end states and therefore that there has already been a supernova within the region.

Despite this there are not any known supernova remnants (SNRs) within Cyg OB2, although Green (2009) does list a number of SNRs in the wider Cygnus X region (including γ-Cygni, G78.2+2.1) but these are predominantly to the west in the Cyg OB1 and OB9 associations. Butt et al. (2003) identified an expanding shell of molecular and dusty material from CO J = 1–0 and IRAS observations that they suggested could be due to an SNR, but this cavity could equally arise from the combined effects of ionizing radiation and stellar winds from the massive stars in Cyg OB2. The lack of observed SNRs in Cyg OB2 does not however rule out any previous supernovae in the association, since SNRs typically remain visible for only ∼10 000 yr (e.g. Levenson et al. 1998) and may not leave any detectable signature when expanding into hot, low-density media such as the cleared-out cavities in OB associations (Chu 1997).

There have been a number of very-high-energy sources identified towards Cyg OB2 that have been considered as possible compact stellar remnants (e.g. Bednarek 2003; Butt et al. 2006), but confirming their nature or identifying counterparts at other wavelengths has proved difficult. Abdo et al. (2009) discovered an ∼115 800 yr old γ-ray pulsar (J2032+4127) with a frequency of 6.98 Hz towards Cyg OB2, which has since been detected in the radio with the same frequency (Camilo et al. 2009) and is also believed to be associated with the previously unidentified TeV source J2032+4130 (Aharonian et al. 2002) as its pulsar wind nebula. Camilo et al. (2009) argue that the pulsar is likely to be located within the Cyg OB2 association, and combined with our newly derived age of ∼1–7 Myrs this is consistent with the pulsar originating within Cyg OB2 supporting the theory that the association has already seen it's first supernova.

Another indication of past supernovae is the presence of high-velocity ‘runaway’ OB stars that are ejected when their binary companion explodes as a supernova (Blaauw 1961). Comerón & Pasquali (2007) identified an O4 If runaway star with a proper motion of ∼40 km s⁻¹ moving away from Cyg OB2 and argued that the star was likely ejected from the association ∼1.7 Myr ago. The authors estimate a mass of 70 ± 15 M_⊙ and an age of ∼1.6 Myr (by comparisons with non-rotating stellar isochrones) suggesting the star was ejected from Cyg OB2 shortly after it was born, and therefore not from a supernova explosion. The other mechanism for stellar ejection is via dynamical encounters in dense stellar systems (Poveda, Ruiz & Allen 1967), yet this process requires much higher stellar densities than are believed to have existed in Cyg OB2 (Wright et al. 2014c) making it unlikely. Using the rotating stellar models of Ekström et al. (2012), and with updated photometry we recalculated the properties of this star. We find B − K_s = 5.27 mag (using Tycho-2 and 2MASS photometry; Høg et al. 2000; Cutri et al. 2003), which for a star of type O4 If with (B − K_s)₀ = −1.15 mag (Martins & Plez 2006) gives A_V = 5.22 mag (using the reddening law determined above) and therefore log L_bol = 5.85 L_⊙ and T_eff = 40 422 K. This implies an initial mass of ∼55 M_⊙ and an age of 4–5 Myr. This would suggest the star was ejected from Cyg OB2 when it was about ∼3 Myr old, the age at which the most massive star in a population might explode as a supernova. This is therefore consistent with a scenario whereby the star's very massive (>100 M_⊙) companion exploded as a supernova and ejected its companion from the association. The remnant of the massive star was presumably ejected from Cyg OB2 in the opposite direction and could be searched for to verify this scenario.

7 THE TOTAL MASS OF Cyg OB2

Using a sample of massive stars that is believed to be complete within a given mass range, and the assumption that the IMF has a ‘universal’ form in all star clusters and OB associations, it is possible to estimate the total mass of a group of stars beyond that which have been detected. To do this for Cyg OB2, we have selected a mass range of 20–40 M_⊙ for which we believe our sample is complete (at least for stars aged <6 Myr). Above 40 M_⊙ our sample will be incomplete for stars aged <6 Myr and below 20 M_⊙ we have argued that our sample is likely to be observationally incomplete. Based on the stellar masses determined above the number of stars in this mass range is |$36^{+1}_{-4}$| (uncertainties calculated from MC simulations).

To estimate the total mass of Cyg OB2, we ran an MC simulation by sampling stellar masses from a ‘universal’ IMF with α = 2.30 (Γ = 1.3) and β = 1.40 (Kroupa et al. 2001; Maschberger 2013) across the entire mass range. We proceeded to count the total mass of stars produced and the number of stars in the mass range of 20–40 M_⊙, until the latter number reached the observed number of stars in that mass range. We performed this simulation 10 000 times, recording the total stellar mass that was necessary in each iteration to produce the observed number of stars. Furthermore, we varied the observed number of stars, |$36^{+1}_{-4}$|⁠, in each iteration according to the uncertainties on that number. This MC simulation therefore takes into account both the uncertainty in the total number of stars observed, and the variation in the number of stars in a cluster of a given mass due to the random sampling of the mass function.

We find that the observed number of massive stars in Cyg OB2 in the mass range of 20–40 M_⊙ can be reproduced if the association has a total stellar mass of |$16500^{+3800}_{-2800}$| M_⊙. This is based on three assumptions. First that we are complete in the chosen mass range, as argued above. Secondly, that the vast majority of star formation in Cyg OB2 occurred in the last 6 Myr (broadly supported by the age distribution and the mass function simulations, though it is difficult to test this assumption without extending the completeness of the sample to lower masses). The final assumption is that the IMF in Cyg OB2 is well represented by the ‘universal’ IMF of Kroupa et al. (2001). This is supported by our IMF simulations at high stellar masses and that found at lower masses (Wright et al. 2010b); however, a small variation can lead to large differences in the total stellar mass. For example increasing (decreasing) the high-mass exponent α by 0.1 increases (decreases) the total stellar mass by 5600 (3900) M_⊙. However, for comparison with other young clusters and associations that have also adopted such an IMF, we believe this is the most appropriate choice.

Our mass estimate is smaller than some previous mass estimates of Cyg OB2, notably that of Knödlseder (2000) who estimated a mass of (4–10) × 10⁴ M_⊙ for the association. Knödlseder (2000) based their estimate on the number of stars brighter than a given magnitude (equivalent to an F3 V star) from a background-subtracted 2MASS star counts study. The choice of a background region in this crowded and complex area of the Galactic plane is fraught with difficulties and could lead to large uncertainties in any resulting star counts or mass estimates. Knödlseder (2000) also estimate the total number of O stars in Cyg OB2 to be 120 ± 20, an estimate almost double the size of our census and one which other studies have not supported (e.g. Hanson 2003). Our mass estimate is more consistent with recent estimates based on studies of lower mass stars, e.g. Wright et al. (2010b) estimated (2–4) × 10⁴ M_⊙ from an X-ray study and Drew et al. (2008) estimated (1–4) × 10⁴ M_⊙ from a photometric study of A-type stars in the association.

8 CONCLUSIONS

We have compiled a census of all known massive stars in the nearby OB association Cygnus OB2. We have gathered data from across the literature to achieve a census containing 169 primary OB stars, including 52 O-type stars and 3 WR stars.

We have used the sample to measure the extinction law towards Cyg OB2, fitting the available photometry to a Fitzpatrick & Massa (2007) R_V = 2.91 ± 0.06 extinction law. This is in good agreement with that determined by Hanson (2003) in the optical and differing only in the near-IR due to shifts in the intrinsic near-IR colours of O-type stars in recent years. We use this new law to derive individual extinctions for all the stars in our sample, the majority of which have extinctions A_V = 4–7 mag, with a typical uncertainty of ∼0.3 mag.

An HR diagram was compiled and compared to both non-rotating and rotating stellar evolution tracks, from which stellar masses and ages are calculated. Uncertainties are derived for all quantities using an MC simulation that takes into account all available uncertainties including quantified classification uncertainties, photometric uncertainties, extinction law uncertainties, estimated spreads in fundamental parameters, and unresolved binaries.

The distribution of stellar masses is assessed using both a binned mass function and a cumulative mass function, both of which are affected by the loss of the most massive stars to stellar evolution. The current mass function is relatively steep with a slope of Γ = 2.2 (in agreement with that found by Kiminki et al. 2007), though when considering either an unbiased sample of stellar masses or by modelling the effects of stellar evolution on the mass function we find it to be consistent with a Salpeter-like ‘universal’ slope of Γ = 1.3.

The non-rotating stellar evolution models suggest a peak stellar age of 2–3 Myr, as derived by previous authors, while the rotating stellar models suggest a slightly older peak age of 4–5 Myr. The results from the rotating stellar models are favoured, both for physical reasons, and because the age distribution that they imply is in better agreement with that found from studies of low- and intermediate-mass stars. When considering only stars in a mass range unaltered by evolutionary effects the age distribution is not inconsistent with a relatively constant level star formation over at least the past ∼6 Myr. Our star formation history fits to the mass function of stars in Cyg OB2 suggest star formation most likely began ∼7 Myr ago. The lack of molecular material or dust in the immediate vicinity of the associations suggests star formation likely ceased at least ∼1 Myr ago, and this is also supported by our mass function models. Our best estimate is that the majority of star formation occurred during the period of 1–7 Myr ago, with a possible peak in the star formation intensity between 4–5 Myr ago.

Given the evidence from both the age distribution and the mass function that some of the most massive stars have already evolved to their end states, we searched the literature for evidence of this and found both a young pulsar and a runaway O star that may have originated in Cyg OB2, consistent with the view that Cyg OB2 has already seen its first supernova.

Finally, we use an MC simulation to estimate the total mass of the association, taking into account the loss of the most massive stars, and estimate a total stellar mass of |$16\,500^{+3800}_{-2800}$| M_⊙ (assuming a Kroupa et al. 2001 IMF) that has formed over the past 6 Myr. This is at the lower end of many previous estimates, based on the number of low- and intermediate-mass stars, but consistent given the uncertainties involved in extrapolating the IMF.

In conclusion, our new census of massive stars in Cyg OB2 provides an updated view of the high-mass stellar population in one of the largest groups of young stars in our Galaxy and for the first time these results show a good agreement with those determined for low- and intermediate-mass stars, thereby providing a coherent picture for this massive OB association. In a future paper, we will use new photometry and spectroscopy to search for previously undiscovered massive stars, extending this census to both lower luminosity B-type stars in Cyg OB2, and O- and B-type stars over a wider area and further along the sightline.

We would like to thank Facundo Albacete Colombo, Nate Bastian, Mario Guarcello, Artemio Herrero, Ian Howarth, Raman Prinja, and Nolan Walborn for useful discussions and comments on this paper. We would like to thank the referee, Matthew Povich, for a helpful and insightful report that improved the content of this paper. NJW was supported by a Royal Astronomical Society research fellowship. This research has made use of NASA's Astrophysics Data System and Simbad and VizieR data bases, operated at CDS, Strasbourg, France.

REFERENCES

APPENDIX A: UVEX PHOTOMETRY AND TRANSFORMATIONS TO THE UBV SYSTEM

Figure A1.

Calibration of transformations between the Massey & Thompson (1991) UBV photometry and the UVEX Ug′r′ photometry. All stars are shown as black dots, with the stars listed by Massey & Thompson (1991) as ‘blended’ shown with red dots. The best-fitting transformation obtained using a least-squares fit that excludes the ‘blended’ stars are shown with red dashed lines and are listed in each panel with the number of stars used for the fit and the dispersion on the best fit.

Open in new tabDownload slide

APPENDIX B: THE FINAL CATALOGUE

The full census of known OB stars in Cyg OB2 is presented in Table B1 in order of increasing primary spectral subtype and with derived observational (A_V), physical (T_eff, L_bol), and stellar (M, age) quantities (the latter calculated using the Ekström et al. 2012 rotating evolutionary models). Stars that fall to the left of the zero-aged MS in the HR diagram are listed as having ages of 0 Myr, while those whose positions in the HR diagram suggest ages >10 Myr are listed as such (see Section 4.5 for discussion of these issues). References for spectral classifications are noted at the end of the table. The full table, which also contains UBV and JHK_s photometry, stellar quantities derived using the Lejeune & Schaerer (2001) non-rotating models, along with 1σ uncertainties on all derived properties, is available from Vizier.