© ESO, 2015

1. Introduction

The change of a period (Ṗ) in a binary system has been commonly evaluated by measuring the times of minimum of the primary eclipses and then comparing these times to the calculated ephemerides for the eclipses. The period change value is then acquired by fitting a polynomial model through the observed minus calculated (O−C) data points. The value of this period change can indicate mass loss through stellar winds or mass transfer between the binary components (Singh & Chaubey 1986).

Mass loss in binary systems mostly comes from the mass transfer via Roche-lobe overflow when the primary star’s photosphere crosses the inner Lagrange point and starts to exchange mass with its companion. While the mass is transferred, the centre of mass of the interacting binary system shifts towards the mass-receiving star. As a result of the conservation of angular momentum, this mass transfer causes a change in the system’s orbital period: the period increases while the mass is transferred to the more massive star, and vice versa. The period also increases when mass is lost from the system through stellar winds.

Sana et al. (2012) studied a sample of massive O stars and found that over 70% of them exchange mass with a companion. The speckle interferometric survey of Galactic massive stars (Mason et al. 2009) expanded upon a decade-old survey and noted that 75% of the massive stars in clusters and 59% in the field have either a spectroscopic or a visual companion. As low inclination or long periods make the binary systems hard to detect, these values are only the lower limit and, when corrected to the observational bias, the upper limit for binarity fraction in bound systems can reach as high as ~ 90 ± 10% (Kiminki & Kobulnicky 2012).

Massive stars have a relatively short lifetime that is influenced by their mass-loss rates. This mass loss determines the evolution of their temperature and luminosity, and also the supernova explosion in the final evolution phases. Massive stars provide feedback in the form of ultraviolet (UV) radiation, stellar winds, and supernovae, which affect the evolution of disc galaxies by mixing and driving turbulence of interstellar gas and by triggering star formation. They also enrich the interstellar medium with heavy elements and drive the chemical evolution of galaxies and the metallicity evolution of the universe. Therefore, knowing the evolution of massive stars gives us insight into the evolution of galaxies. The uncertainties in determining the mass-loss rates for massive stars are still large even in the local Universe (Smith 2014).

Constraining orbital parameters for massive binary systems is important to test the predictions of theoretical models. To determine all the parameters of a single binary system, one needs to obtain photometric and spectroscopic data at various orbital phases, which requires great observational effort. As it is time-consuming to obtain those observations, and as the eclipsing nature is dependent on the inclination of the system, there are not many massive binaries with fully constrained parameters known. In this work we present the stellar parameters as well as the period and the period change values of seven selected systems in Cygnus OB2.

The stellar association Cygnus OB2 is one of the most massive and richest associations known in our Galaxy, containing at least 60 O-type stars (Negueruela et al. 2008). Determined from the stellar evolution models, the age of the cluster is found to be 1−7 Myr (Wright et al. 2015). It is a dynamically young sub-structured and unbound association that is located in the molecular cloud complex Cygnus X, and its total estimated mass is (1.6 ± 0.3) × 10⁴M_⊙ (Wright et al. 2014).

The Cygnus OB2 Radial Velocity Survey (RVS) of 128 selected O and early-B stars is underway (Kiminki et al. 2007, 2008, 2009, 2012; Kobulnicky et al. 2014) to assess the binary properties of massive stars within a common formation environment at a similar age. There are now 48 massive binary systems known in Cygnus OB2.

We have chosen Cygnus OB2 to study the mass loss in massive binaries because of the previously published data for massive stars in that association. Based on the radial-velocity (RV) data and our own photometric data that we obtained during 2011−2014, we have chosen seven massive binary systems with eclipsing light curves. To assess the change of period in those systems, we use a Bayesian parameter estimation method. We also use a binary modelling software to acquire orbital and physical parameters. For several stars that were previously modelled using RV data only, we present models with additional photometric data.

2. Data

We here combine new photometric observations with photometry and RV data from the literature.

References to the data that we used for each star are given in Table 1. We use the same naming conventions as in the Cygnus OB2 RVS papers: the MT names from Massey & Thompson (1991), Schulte from Schulte (1958), and A36 and B17 from Comerón et al. (2002). The data sets are described in more detail below.

We monitored Cygnus OB2 with two telescopes during 2011−2014. The first data set was obtained with the commercial iTelescope.net 0.25-m Takahashi Epsilon telescope (T4), located in Mayhill, New Mexico, USA. An SBIG ST-10 CCD camera with Johnson-Cousins BVI filters was used, providing a field of view (FOV) of 60.5 × 40.8 arcmin. Observations were carried out in 344 nights during 2011−2013. An observation consisted of three consecutive exposures in each filter. Data reduction was carried out using the PHOTWORK software (Tuvikene 2012). Magnitudes were extracted by means of aperture photometry, except for MT720, which required PSF photometry because of a neighbouring star. To minimise observational errors, we used the nightly mean instrumental differential photometry values (for details see Laur et al. 2012).

The second data set was obtained in 103 nights during two seasons, 2013−2014, using the Tartu Observatory 0.31-m RAITS telescope, located in Tõravere, Estonia. The telescope was equipped with an Apogee Alta U42 CCD camera and Johnson-Cousins BVRI filters. The FOV was 37 × 37 arcmin. Five consecutive exposures were taken with each filter during an observation. The CCD frames were reduced and magnitudes were extracted by aperture photometry, as described above.

Figure 1 shows the FOV of the two telescopes we used. Two of the seven programme stars, A36 and B17, fell outside of the RAITS FOV.

Fig. 1 Full frame of the Cygnus OB2 region from our first data set (iTelescope) with the selected binaries shown in red. The dashed line indicates the FOV of our second data set (RAITS).

The Cygnus OB2 RVS is a spectroscopic survey of unbiasedly chosen 128 stars of which 83 are B-type and 45 are O-type stars. They have used Keck+HIRES spectrograph (1999−2001), the Lick+Hamilton echelle spectrograph (1999−2000), the WIYN+Hydra spectrograph (2001−2008), and a longslit optical spectrograph with a 2.3 m telescope at the Wyoming Infrared Observatory (2011−2013, 2014 May). The Cygnus OB2 RVS data that we used were attained between 1999 to 2009. The longest data set of 43 measurements (1999−2007) is that for MT059. Schulte 3 and A36 have the fewest observations (~15) and a time base of one year, while MT720 has 32 observations in two years and MT696 has 31 measurements with the time base of 30 days.

The Northern Sky Variability Survey (NSVS) covers the entire northern sky in an unfiltered optical region, ranging from 8 to 15.5 magn from the Robotic Optical Transient Search Experiment (ROTSE-I; Woźniak et al. 2004). Except for MT696, 170 data points were obtained within one year (April−December, 1999) for all of our selected stars. The data set for MT720 was available, but due to large photometric errors, the signal of its eclipses was unattainable, and we had to exclude it from our analysis. The errors for NSVS data sets are typically larger than our own data sets because their data quality is limited by a severe blending that is due to the 14 arcsec pixel size.

The Hipparcos catalogue data (van Leeuwen 2007) are available only for the brightest star in our selection, Schulte 5, because of the limiting magnitude of Hipparcos observations. The available data span over three years, from the end of 1989 to the beginning of 1993, with one hundred broad-band measurements.

The Optical Monitoring Camera (OMC) on board the ESA International Gamma-Ray Astrophysics Laboratory (INTEGRAL) mission consists of 5° × 5° FOV photometric camera observing in the Johnson V band (Mas-Hesse et al. 2003). The only stars in our selection that have been observed by IOMC are Schulte 3 and Schulte 5. There have been more than 5000 observations made between 2002 and 2009 (which is our longest time series from a single source), and their optical accuracy is roughly on the same level with NSVS photometry.

Souza et al. (2014) have been monitoring MT696 since 2010 in two 5 nm wide filters centred on the continuum and Hα (645 and 656 nm, respectively). Their photometric data span over three years from 2010 to the end of 2013 and consist of 105 observations. These data and the RV data of Kiminki et al. (2012) are the only published time series of MT696 to date.

For Schulte 5 we have the longest time base, with measurements from 12 different sources. The first published measurements for Schulte 5 are the times of minimum from 1898, although the accuracy of those measurements was doubtful for the authors themselves (Sazonov 1961). During a five-month period in 1952, Miczaika (1953) obtained 36 measurements with a photoelectric equipment on the 32 cm refractor at the Heidelberg Observatory. Hall (1974) observed Schulte 5 photoelectrically for 142 times with a 60 cm Seyfert reflector during the period of 1967 to 1970 using a set of closely matching UBV passband filters. Linder et al. (2009) acquired over 300 exposures during a campaign of 14 nights in 1998 with five filters, three of which were centred on different emission-line wavelengths found in various subsets of Wolf-Rayet objects. The other two were centred on line-free areas (around 5057 and 6051 Å). For our period change analysis we chose only the measurements taken with the continuum narrow-band filters because the emission line data also indicate variability that is not directly related to the eclipses of the binary.

RV data for Schulte 5 are available from three sources (Wilson & Abt 1951; Bohannan & Conti 1976; Rauw et al. 1999). All the observations span three years (1946−1949, 1971−1974, and 1994−1997) and consist of 42, 7, and 25 measurements, respectively. Schulte 5 has the longest time base for RV data because it is the only binary in our selection on which more than one RV observation campaign has been carried out.

Fig. 2 Available data shown in time domain with reduced Julian dates at the bottom and calendar years at the top. Point sizes are indicating the number of measurements in a one-day bin. We have changed the scale of the plot in the middle (dashed line) to better show the data available from 1990 onwards. Green points indicate our photometric data sets, red points are the available RV data, purple points show the times of minimum from Sazonov (1961), and black points are all other available photometric data listed in Table 1.

All the available data points are shown in Fig. 2 to illustrate the time base for each star. It can be seen from the figure that Schulte 5 is by far the most frequently observed star in our selection, while the other massive binaries have scarcely been studied. Our analysis mostly depends on three sources: our own two data sets from 2011 to 2014, NSVS, and Cygnus OB2 RVS. While Schulte 5 has a time base of 116 years, all the other binaries have a time base shorter than 15 years.

The time base is important in period change studies because a constant change in period introduces a gradual shift in eclipse timings for each consecutive orbit. This causes a measurable cumulative phase shift over a long time base. The longer the time base, the more apparent the shift becomes; this increases the accuracy of the period change assessments. Figure 3 illustrates the cumulative effect of the period change for Schulte 5 with two different data sets that were observed 45 years apart.

Fig. 3 V-filter phase diagram of Schulte 5, folded with a constant period. The iTelescope data (black points) and Hall (1974) data (blue points) have been obtained 45 years apart. Solid lines are smoothed empirical light curves. The visible phase shift is due to the cumulative effect of period change (Ṗ).

3. Method

To estimate the variations in orbital period of our selected massive binary systems, we used a Bayesian approach of simultaneously fitting all available observational data: light curves, times of minimum, and RV data.

This method consists of two distinctive steps: acquiring the model curves and then assessing the period (P) and period change (Ṗ) values that simultaneously fit the data to these curves (Table 2). We assumed that the shape of the model curves remain constant throughout the whole data series¹. For the initial assessment of P, we used a phase dispersion minimisation method (Stellingwerf 1978). This initial period was necessary to fold the photometric and RV data for input to the modelling software. For the modelling we used PHysics Of Eclipsing BinariEs (PHOEBE), (legacy version 0.31a), a binary modelling software based on the Wilson-Devinney code (Wilson & Devinney 1971). The modelling was carried out using RV data from Kiminki et al. (2008, 2009, 2012), Bohannan & Conti (1976), Rauw et al. (1999) and photometric data from our 2011−2013 data set. We used this particular photometry data set as it had relatively smaller errors and included all of our selected stars. The spectral classes for the objects were taken from Kobulnicky et al. (2014), and they were used as an input for the initial temperature (Gray 2001) and mass values (Martins et al. 2005) for the models. The resulting model parameters for masses (M₁,M₂), mass ratio (q), radii (R₁,R₂), temperatures (T₁,T₂), and orbital parameters are presented in Table 3.

The phased curves from PHOEBE software were used as input to the Bayesian inference tool called MultiNest (Feroz & Hobson 2008; Feroz et al. 2009, 2013). The Bayesian parameter estimation method is a Monte Carlo algorithm that calculates the ephemerides of the star as $\mathrm{HJD}\mathrm{=}\mathrm{HJ}{\mathrm{D}}_{\mathrm{0}}\mathrm{+}\mathit{PE}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{Ṗ}\mathit{P}{\mathit{E}}^{\mathrm{2}}\mathit{,}$(1)where E is the cycle number.

Using a random walk approach, we assessed the values of P and Ṗ that best set the data to the model curve. By comparing the folded light curve and folded RV curve that were calculated using the ephemerides to the model curves, we obtained the likelihood of the parameter fit. The global maximum of the likelihood posterior is the best fit of the parameters in conjunction with the model curves. We assumed the posterior to be Gaussian, and thus the parameter fit errors are given as a standard deviation. To minimise the errors, we fixed the value of HJD₀ using the minimum of the modelled PHOEBE light curve.

For the period change analysis, we used one model light curve for all the photometric filters because the fitting of the ephemerides mainly depends on the time domain as opposed to the small differences in eclipse depths from different filters.

To test whether the PHOEBE model phased light curve causes any systematic errors in the period change determination, we constructed empirical model light curves for our selected stars. All the results for P and Ṗ were obtained using both the PHOEBE model and the empirical model. The empirical phased light curve, which was sigma clipped and smoothed, was derived from all the photometric measurements where the first estimation of the period change was already accounted for. As there were not enough RV data distributed over the full orbital phase, no accurate empirical RV model curve could be derived.

The errors for each photometric data set were obtained separately as the standard deviation of the difference between the data points and the empirical light curve. The errors for RV data were taken from the respective papers, and the errors for the times of minimum in Sazonov (1961) were given as the timing error of the observations. These errors determined the weight of each data set in the final period change analysis.

Conventionally, O−C diagrams are constructed using only precisely measured photometric times of the eclipses where the times are obtained by continuous monitoring of the selected star during an eclipse. Any RV data are usually completely ignored or only the times of conjunctions are used in the O−C analysis. In our case, however, RV data can be fitted together with any photometric data. This method is especially useful for sparsely covered time series of the eclipses and only approximate times of the minima. This can either be due to observational difficulties or when working with older data that have not been optimised in time domain for this type of studies. In general, our method can be used for any data coming from all-sky-covering campaigns. They usually provide unevenly spaced observation epochs that make it harder to detect the times of minimum but on the other hand have a long time-series for many objects. Our simultaneous fitting of data with the appropriate weights provides a solution that is not dependent on the measurement errors of any single minimum and also improves the statistical significance of the period assessment because more data points are used. Previous works for finding period change using all of the available data have been carried out for RR_d type stars (Paparo et al. 1998) and for selected eclipsing binaries (Wilson & Van Hamme 2014).

4. Results

We present the estimated orbital ephemerides in Table 2 and binary parameters in Table 3. We give one-sigma likelihood-estimation uncertainties for P and Ṗ, formal errors from PHOEBE for T₂,i,e and calculated formal errors² for M₁,M₂,R₁,R₂ and a.

Our approach in this work does not take into account any cyclic patterns of period changes that may come from the light travel time effect through the presence of a third body in the system or from stellar magnetic activity cycles (Liao & Qian 2010). The data for most of our selected stars have a relatively short time base compared to the typical periods of the third orbiting body or magnetic cycles, therefore we have opted to exclude a cyclic pattern as a free parameter because this would artificially introduce a new degeneracy to the fit. Any cyclic patterns in the period change are sources of possible errors in our results.

Neither did we take into account the apsidal motion in binaries with eccentric orbits (MT720 and MT059). The apsidal motion would visibly alter the shape of light and RV curves, but we did not see this in our data sets (see Appendix A).

Fig. 4 Parameter estimates for period and period change. Periods are given as a difference from the Pemp value (see Table 2). The coloured solid line represents the results for an empirical phase model, and the dashed line represents the results for the calculated PHOEBE phase model. The inner contour is the 68% and outer contour is the 95% confidence level.

Figure 4 shows the parameter P and Ṗ estimates for both the empirical and the PHOEBE model folded light curves. The shape of the resulting parameter estimation likelihood cloud shows the amount of degeneracy of the two parameters as well as the uncertainties in the parameter estimations. As a result of the degeneracy of P and Ṗ, the given uncertainties for these parameters are determined by marginalised distributions. If we had also fitted the HJD₀ value, the results would have been degenerated between all three of those parameters. Therefore, the value of HJD₀ was fixed to reduce the systematic errors of the two other fitted parameters that we are interested in. The distance between the locations of the parameter estimation likelihood clouds shows the systematic errors that arise from different model curve shapes, and the sizes of those clouds show the statistical errors. Clearly, one has to pay attention to the model folded light curves because they introduce some systematic errors compared to the true light curve.

Appendix A describes the model fits for each of the dataset separately. The visual comparison between the data and the models can be used to assess the applicability of the final resulting parameters.

By comparing the results to each other from the empirical and PHOEBE folded light curve models for each star in Fig. 4, the results look consistent for most cases. The results for Schulte 5 seem to disagree with each other, but because of small uncertainties in the parameter estimation, they agree well. The resulting uncertainties for Schulte 5 are very small because its time base is the longest of our selected objects. The results for MT720 show a very large systematic error between the two models. This is due to the relatively short time base as well as the difference between the shape of these two model folded light curves (see also Fig. A.4).

For four stars, we detect a non-zero period change values. The strongest evidence for period change is found for Schulte 5, while MT696 also shows period change at three-sigma confidence level and A36 and MT059 at two-sigma confidence level. MT059 and MT696 exhibit a negative period change, while A36 and Schulte 5 have a positive period change. Our calculations for B17, MT720, and Schulte 3 do not indicate any change in their periods.

The stellar parameters (Table 3) were obtained with the PHOEBE software where RV data were fitted alongside photometric data. The stellar parameter values were not fine-tuned for each individual star as they are a by-product of acquiring model phase curves. The parameter modelling is degenerate between many of the parameters, which means that adding high-resolution spectra for precise temperature and spectral-type measurements would constrain the models further.

There are no parameters for B17 because the only available RV data from Comerón et al. (2002) are published as a phased graph with no available ephemerides. Without RV data, no orbital parameters or masses can be modelled accurately. In addition, the RV data for MT059 are only available for the primary component, and thus the model stellar parameters are only one of the possible parameter fits because of the degeneracy of the mass ratio. All the fits of our data to the acquired binary models are presented in Figs. 5−10.

Fig. 5 Binary modelling for Schulte 5 with the PHOEBE software. The upper panel shows our photometric iTelescope data (data set 1) with the model curve (solid line), and the lower panel shows combined radial-velocity data from Bohannan & Conti (1976) and Rauw et al. (1999) with the model curves for both primary and secondary components.

Fig. 6 Binary modelling for A36 with the PHOEBE software. Designations of lines, labels, and photometric data set are the same as in Fig. 5. Radial-velocity data are taken from Kiminki et al. (2009).

5. Discussion

We used our estimated period change values to derive mass-loss rates for two distinctive cases. For an assumed conservative mass exchange between the two components with constant system orbital angular momentum and system mass, the mass transfer rate can be calculated as case I in Singh & Chaubey (1986): $\dot{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\frac{{\mathit{M}}_{\mathrm{1}}{\mathit{M}}_{\mathrm{2}}}{{\mathit{M}}_{\mathrm{1}}\mathrm{-}{\mathit{M}}_{\mathrm{2}}}\frac{\mathit{Ṗ}}{\mathit{P}}\mathrm{·}$(2)We can also assume that the mass is lost solely through stellar wind (case IV in Singh & Chaubey 1986): $\dot{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{M}}_{\mathrm{1}}\mathrm{+}{\mathit{M}}_{\mathrm{2}}\mathrm{)}\frac{\mathit{Ṗ}}{\mathit{P}}\mathrm{·}$(3)Table 4 summarises for all of our selected stars the primary component mass-loss rates calculated for both conservative mass exchange and mass loss through stellar wind using the two acquired period change values of Ṗ_emp and Ṗ_model. Equations (2) and (3) are for the idealised occurrences where no other effect interferes with the period-change value. The primary component mass-loss values in Table 4 represent isolated cases through either conservative mass exchange or stellar wind. They do not indicate the real values of the mass loss as in reality the mass loss is a combination of different factors. To constrain the mass-loss values further, some additional method to find the strength of stellar wind (radio observations or analysis of the stellar-wind-induced spectral line profiles) has to be included.

In the next sections we compare our results to the previous findings and present some of the problems that we came across for all stars individually.

Fig. 7 Binary modelling for MT059 with the PHOEBE software. Designations of lines, labels, and photometric data set are the same as in Fig. 5. Radial-velocity data are taken from Kiminki et al. (2008).

Fig. 8 Binary modelling for MT696 with the PHOEBE software. Designations of lines, labels, and photometric data set are the same as in Fig. 5. Radial-velocity data are taken from Kiminki et al. (2012).

Fig. 9 Binary modelling for MT720 with the PHOEBE software. Designations of lines, labels, and photometric data set are the same as in Fig. 5. Radial-velocity data are taken from Kiminki et al. (2008).

Fig. 10 Binary modelling for Schulte 3 with the PHOEBE software. Designations of lines, labels, and photometric data set are the same as in Fig. 5. Radial-velocity data are taken from Kiminki et al. (2008).

5.1. Schulte 5

Schulte 5 (V729 Cyg, BD +40° 4220) is an eclipsing binary of two O-type supergiants (Kobulnicky et al. 2014). There is also evidence for a third non-thermal radio source around the binary with a period of 6.7 years (Kennedy et al. 2010; Cazorla et al. 2014).

Linder et al. (2009) assumed the orbital solution from Rauw et al. (1999) to derive the absolute masses of the components (M₁ = 31.9 ± 3.2M_⊙ and M₂ = 9.6 ± 1.1M_⊙). Our modelling reveals similar masses for the components, 38.2 ± 5.2M_⊙ and 10.2 ± 1.4M_⊙, respectively.

Martins et al. (2005) have calculated a grid of models for O type stars using CMFGEN (Hillier & Miller 1998) atmosphere code, and they include theoretical spectroscopic masses for them. For the primary component of spectral type O7I, the spectroscopic mass would be M_spec = 40.91M_⊙. The authors note that one should take care in dealing with their spectroscopic modelled masses as the uncertainty may be as high as 35 to 50%. Nevertheless, our resulting mass for the primary is closer to the theoretical value than to the result from Linder et al. (2009).

We chose Schulte 5 as a testbed for our period change study because the eclipse-timing method for period change as well as the mass-loss rates have been studied on multiple occasions. Previous assessments for the mass-loss rates were derived by Persi et al. (1990) from infrared observations Ṁ = 2.5 × 10^-5M_⊙ yr^-1, by Conti & Howarth (1999) combined with multi-frequency radio observations Ṁ = (3.7 ± 1.3) × 10^-5M_⊙ yr^-1, and by Linder et al. (2009) using the period change rate of observed times of minimum Ṁ = (2.1 ± 0.6) × 10^-5M_⊙ yr^-1.

Following the calculation from Eq. (2), Linder et al. (2009) found a mass-loss rate $\dot{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\mathrm{(}\mathrm{4.7}\mathrm{\pm }\mathrm{1.5}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-6}{\mathit{M}}_{\mathrm{\odot }}\mathrm{y}{\mathrm{r}}^{-1}$, which they argued against because this value was lower than expected from the previous assessments derived from the infrared and radio data. Using the masses obtained from our PHOEBE model and the P and Ṗ from our Bayesian parameter estimation method, we also obtain similar mass-loss rates of $\dot{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\mathrm{(}\mathrm{3.0}\mathrm{\pm }\mathrm{0.2}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-6}{\mathit{M}}_{\mathrm{\odot }}\mathrm{y}{\mathrm{r}}^{-1}$ and $\dot{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\mathrm{(}\mathrm{3.4}\mathrm{\pm }\mathrm{0.2}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-6}{\mathit{M}}_{\mathrm{\odot }}\mathrm{y}{\mathrm{r}}^{-1}$ depending on the applied period change value. This difference between the Linder et al. (2009)$\dot{{\mathit{M}}_{\mathrm{1}}}$ value and our values is mostly due to the difference in the period change values. The value from Linder et al. (2009), Ṗ = (1.9 ± 0.5) × 10^-8, is higher than our obtained value of Ṗ_emp = (1.16 ± 0.03) × 10^-8 and thus their mass-loss rate is also slightly higher. Our period change value is probably more accurate because we have a time base that is longer by over 40 years.

Yaşarsoy & Yakut (2014) found a period change of P/Ṗ ~ 9.4 × 10⁵ from eclipse timings by acquiring new observations on 65 observing nights from 2010 to 2011 in Bessel UBVRI bands. With their given period of P = 6.597981 d, we can calculate Ṗ = 7.01913 × 10^-6, which is much higher for the value of Ṗ = (1.16 ± 0.03) × 10^-8 that we obtained. It may be possible that their given P/Ṗ value is in the units of years instead of days as one might assume from their period value. In this case, the value Ṗ = 1.92305 × 10^-8 agrees with the value of Ṗ from Linder et al. (2009). They also modelled the masses of the primary (31.6 ± 2.9 M_⊙) and secondary (8.8 ± 0.3 M_⊙) component. If we use their data for calculating mass loss through mass transfer, we obtain $\dot{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\mathrm{4.325}\mathrm{\times }{\mathrm{10}}^{-6}{\mathit{M}}_{\mathrm{\odot }}\mathrm{y}{\mathrm{r}}^{-1}$, which is very similar to the mass-loss rate that Linder et al. (2009) obtained, but differs from ours due to the value of period change and obtained component masses.

We find that the mass-loss rate through stellar wind (Eq. (3)) for Schulte 5 is $\dot{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1.56}\mathrm{\pm }\mathrm{0.07}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-5}{\mathit{M}}_{\mathrm{\odot }}\mathrm{y}{\mathrm{r}}^{-1}$ or $\dot{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1.74}\mathrm{\pm }\mathrm{0.08}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-5}{\mathit{M}}_{\mathrm{\odot }}\mathrm{y}{\mathrm{r}}^{-1}$ depending on the Ṗ value, which are comparable to the result of Linder et al. (2009; $\dot{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{2.1}\mathrm{\pm }\mathrm{0.6}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-5}{\mathit{M}}_{\mathrm{\odot }}\mathrm{y}{\mathrm{r}}^{-1}$). The previous mass-loss rate assessments of the stellar wind from radio and infrared observations have confirmed the value to be close to our finding. This indicates that the period change of the system is caused by the mass loss from stellar wind of the primary component instead of the mass transfer from the secondary component.

5.2. A36

A36 was first named in a paper by Comerón et al. (2002), who listed it as a possible binary. The eclipsing nature of A36 was found from NSVS data (Otero 2008; Hoffman et al. 2009) and it is classified as an Algol/β Lyrae candidate. Kiminki et al. (2009) determined the spectral type as a B0Ib+B0III system on the basis that the light curve and spectra show similar component temperatures and luminosities with differences in specific lines.

Our calculated mass ratio for the binary system is 0.66 and is within error limits of the Kiminki et al. (2009) estimation of 0.70 ± 0.06. They give a period of P = 4.6752 ± 0.0007 d that is slightly different from our result of P = 4.67604 ± 0.00002 d. This is mainly due to the lack of observations in the RVS for A36 to precisely determine the period. We find the period change from the combination of two models to be Ṗ = (3.85 ± 1.3) × 10^-8, thus the mass loss by direct mass transfer would be $\dot{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\mathrm{(}\mathrm{4.9}\mathrm{\pm }\mathrm{2.0}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-5}{\mathit{M}}_{\mathrm{\odot }}\mathrm{y}{\mathrm{r}}^{-1}$ and by stellar wind $\dot{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{5.7}\mathrm{\pm }\mathrm{1.9}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-5}{\mathit{M}}_{\mathrm{\odot }}\mathrm{y}{\mathrm{r}}^{-1}$.

The value of $\dot{{\mathit{M}}_{\mathrm{1}}}$ by stellar wind is too high for a B0Ib star because the maximum empirical estimates for mass-loss rates are typically lower than 10^-5M_⊙ yr^-1 (Vink et al. 2000; Puls et al. 2006). Therefore, the change in period must be caused by the mass transfer from the secondary to the primary. This also indicates that the secondary component has to be more evolved than the primary. Although the spectral classification by Kiminki et al. (2009) contradicts this conclusion, it may still be viable if the luminosity class of the primary component could be lowered. Our PHOEBE modelling results for the stellar radii and masses also support the argument that the primary component is closer to the luminosity class III.

5.3. B17

B17 was first discussed by Comerón et al. (2002), who discovered that it displays emission in the Brackett γ line (the B notation in front of the star number), which is an indicator of high mass-loss rate typical of the evolved stellar stages. It is a massive evolved binary system with a spectral classification estimation of O7Ia+O9I supergiants. Hoffman et al. (2009) classified it as an Algol/β Lyrae type using NSVS data.

For B17 we could not use any spectroscopic data. Stroud et al. (2010) have done an extensive photometric and spectroscopic campaign on B17, but their data are unfortunately not publicly available and they are only presented in a phased form on a graph. Their determined period of 4.02174 ± 0.00003 d coincides with ours, and because the time base between the NSVS photometry data and their observed spectra was only four years, they did not discuss any possibility of a period change. We did not find any indication of a period change within the binary system with values of Ṗ_emp = (1.1 ± 1.6) × 10^-8 and Ṗ_model = (−0.8 ± 1.7) × 10^-8, and therefore we find no mass loss in the system. No stellar parameters were modelled either because there were no available RV data.

5.4. MT059

MT059 (Schulte 1) is an O8V type star with a B, possibly later than B1V, type companion. Kiminki et al. (2008) have only measured radial velocities for the primary component and did not detected any spectral features of the companion, which indicates that the luminosity ratio has to be lower than L₂/L₁ ≲ 0.2. For their primary component mass assessment, they used the theoretical mass for a spectral type O8V from Martins et al. (2005) with a value of M_spec = 21.95M_⊙. They concluded that due to the absence of spectral signatures of the secondary component, the interpolation from its spectral type yields the mass of the secondary companion to be 5.1M_⊙ ≲ M₂ ≲ 13.8M_⊙ and the mass ratio 0.24 ± 0.02 ≲ q ≲ 0.64 ± 0.03. Our resulting mass ratio falls on the lower side of their prediction with a value of 0.29, and the masses we obtained agree with the previously mentioned predictions with values of M₁ = 21.23 ± 6.24M_⊙ and M₂ = 6.22 ± 1.83M_⊙. There are no RV data for the secondary companion, therefore we cannot independently assess the mass ratio, and thus the mass ratio is degenerate with the other stellar parameters. The semi-major axis of the system also yields large formal errors, which lead to large errors in the mass assessments. The modelled light curve (upper panel in Fig. 7) indicates the presence of a reflection effect (phase ~0.2) that is only loosely supported by the data. A more thorough investigation is needed to determine whether this effect is from the observational uncertainties or from the underlying physics of the system.

MT059 is indicated by Kiminki et al. (2008) as a SB1−type binary without mentioning its eclipsing nature. Our observations and NSVS photometric data show clear eclipses from the photometric folded light curve. The binary has an eccentricity of 0.14, which is also notable from the uneven timings in phase space between the primary and secondary eclipse.

We find that of all of our selected binary systems, MT059 has the highest period change value (Ṗ_emp = (−4.7 ± 1.7) × 10^-8 or Ṗ_model = (−3.6 ± 1.7) × 10^-8). This is due to the RV observations, which are also visibly shifted in our PHOEBE modelling (Fig. 7). The mass loss through conservative mass exchange with the results from our empirical model is $\dot{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{10.5}\mathrm{\pm }\mathrm{4}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-6}$M_⊙ yr^-1. As the period change is negative, the mass loss from stellar wind cannot be calculated using Eq. (3)because it would result in a positive $\dot{{\mathit{M}}_{\mathrm{1}}}$ value, and therefore the period change must be from mass transfer from the primary to the secondary component.

5.5. MT696

MT696 is a binary system of an O9.5V spectral type primary and an early-B type companion. Rios & DeGioia-Eastwood (2004) were the first to report the eclipsing binary status of MT696, indicating that it might be a W UMa type binary with a 1.46-day period. Kiminki et al. (2008) argued that the early spectral types of the components, the absence of emission features, and a period longer than one day suggest MT696 to be β Lyr type instead. They also improved the period estimate to 1.4692 ± 0.0005 days and estimated the inclination to be near 80° and the mass ratio q to be 0.85 ± 0.03, which are close to our estimates. Their minimum mass estimates are 15.1 ± 0.7 M_⊙ for the primary and 12.8 ± 0.5 M_⊙ for the secondary component, while our modelling gives total masses of 14.00 ± 0.81 M_⊙ and 12.67 ± 0.73M_⊙, respectively. According to Martins et al. (2005), the mass for spectral type O9.5V is M_spec = 16.46 M_⊙. Our mass estimates are within errors of the minimum masses from Kiminki et al. (2008), but lower than the theoretical primary mass, which means that our resulting masses might be underestimated. Another possibility is that the spectral types and thus the masses of the binary components might not be accurately determined by Kobulnicky et al. (2014).

Souza et al. (2014) found the period to be 1.46919 ± 0.00006 d, and they did not find any period change in their two-year-long data set. We, on the other hand, found the period change at three-sigma confidence level with the value of Ṗ_emp = (−0.37 ± 0.08) × 10^-8 and Ṗ_model = (−0.46 ± 0.08) × 10^-8. This is the lowest value of a period change amongst our studied binaries. As a result of the low value for the period change and the similar masses for the binary components, the errors for conservative mass-loss rates are very high because the total error depends on the difference between the companion masses. The acquired mass-loss rate for direct mass transfer is then $\dot{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{0.4}\mathrm{\pm }\mathrm{1.6}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-4}$M_⊙ yr^-1, and it is dominated by the measurement errors of the stellar masses.

5.6. MT720

Kiminki et al. (2012) have observed MT720 and measured its radial velocities on 32 occasions, noting that no known photometric variability have been studied. They indicated the period of the system to be 4.3622 ± 0.0003 days (compared to our P_emp = 4.36197 ± 0.00001 and P_model = 4.36177 ± 0.00001 days) and an eccentricity of 0.35 ± 0.02, which is evident from the RV curve. They also found minimum masses from the RV curve and suggested the spectral type of the binary system to be B0−B1V + B1−B2V. Their calculated minimum masses of M₁ ≥ 15.5M_⊙ and M₂ ≥ 11.1M_⊙ are in accordance with our modelling results for the component total masses of M₁ = 18.52 ± 1.51M_⊙ and M₂ = 13.19 ± 1.07M_⊙.

Salas et al. (2014) conducted a 1.5-year-long photometric survey in Cygnus OB2 and have included MT720 as a confirmed eclipsing binary with a period of 4.3619 ± 0.0001 days. Our photometric data also show clear evidence of eclipsing nature of the binary. At first we note that in phase space the primary eclipse coincides with the RV curve when the more massive star is eclipsing the secondary component, indicating that it must be less luminous. After modelling the light curve, it turns out that because of the high eccentricity of the system, the primary eclipse depth lessens into what we see as a secondary eclipse. This effect is due to the combination of the high eccentricity and a low inclination angle, resulting in a partial eclipsing of the primary component by the secondary star. The model in Fig. 9 is acquired by introducing a phase shift of 0.45 so the secondary eclipse would coincide with a phase zero point.

Modelled stellar parameters for MT720 are more prone to errors than the other binary systems because of the poor fit between the data and the model. As can be seen in Fig. 9, the secondary eclipse in the observations is narrower than the model curve. Our period change modelling for this system has a strong systematic effect that is due to the differences between the two phased light curve models. The above-mentioned modelling discord arises because MT720 has the largest photometric uncertainties because it is the faintest object in our selection. The systematic errors between the models prevent us from finding any reliable evidence for the period change in MT720.

5.7. Schulte 3

Schulte 3 is a massive binary system consisting of O6IV: and O9III stellar type components. Kiminki et al. (2008) found a mass ratio of q = 0.44 ± 0.08, which agrees within errors with our result of 0.53 ± 0.05, but differs from the theoretical mass ratio of q = 0.64 for O6IV: + O9III binary system (Martins et al. 2005). Our modelled masses (M₁ = 32.73 ± 7.04M_⊙ and M₂ = 17.21 ± 3.70M_⊙) agree with their minimum mass predictions of M₁ ≥ 17.2M_⊙ and M₂ ≥ 7.6M_⊙. Our acquired primary component mass is also close to the theoretical spectroscopic mass of O6V (M_spec = 31.73M_⊙), but the theoretical mass value for the secondary component (M_spec = 23.07M_⊙) is higher than our result.

Kiminki et al. (2008) argued that the difference between the mass ratios is due to the mass transfer between the components. The period we find coincides with their result (P = 4.7459 ± 0.0003 days) within the errors, but we did not detect any notable period change in the system. Our period change value for empirical model is Ṗ_emp = (−1.5 ± 1.0) × 10^-8, which does not imply a non-zero value at two-sigma confidence level, and for PHOEBE model the value is Ṗ_model = (0.0 ± 0.9) × 10^-8. These values show clear evidence for no period change within the system and indicate no mass transfer between the components.

6. Conclusions

We studied the period change of seven massive binaries in the Cygnus OB2 association. For this we used archival photometry and RV data together with our own photometric observations from 2011−2014. Our observations were essential in the period-change analysis for A36, B17, MT696, and MT720 because they increased the time base up to 13 years. We note for the first time the eclipsing nature of MT059.

Four stars show evidence of a non-zero period change at two-sigma confidence level. For three stars we reported the period change for the first time. Our calculations confirm the previous assessments for Ṗ of a star Schulte 5. We also show that the change in period in Schulte 5 originates from the stellar wind from the primary component. Two of the systems, B17 and Schulte 3, show period change values close to zero, indicating that there is no mass transfer between the binary components or that the period change values are too low for us to estimate them reliably with the available data. Our method had systematic errors in fitting the period and period change parameters for MT720, and thus we did not derive any reliable values for that system. Two of the systems, MT059 and MT696, have negative Ṗ values that show mass transfer from the primary to the secondary component. We argued that the change in period for A36 might stem from mass transfer from the secondary component to the primary, indicating a more evolved secondary stellar component.

We showed that by applying the Bayesian parameter estimation method the period change can be calculated using all the available data simultaneously for any binary star. Using all data points in contrast to using only the times of minimum improves the statistical significance of the calculations and opens up new data sets that can be included in the analysis. This is especially important with the trend to digitise archives of photographic plates because it adds new historical observations that were not optimised for the eclipses of binaries, which means that the times of minimum may be difficult to obtain.

Online material

Appendix A: Individual data sets and model curves

In this section we present the resulting phase diagrams for photometric and RV data sets for our studied stars. The diagrams are constructed by folding the data as Eq. (1) with the resulting orbital ephemerides (Table 2). Figures A.1−A.8 show individual data sets together with the phased model curves.

Fig. A.1 Phase diagrams for A36 using the Pemp and Ṗemp from Table 2. Blue dashed lines in the magnitude plots denote a model light curve from PHOEBE software, and the red solid lines denote an empirical light curve. The solid line in the RV plot represents the model for the primary component, and the red dashed line represents the model for the secondary component.

Fig. A.2 Phase diagrams for B17. Designations are the same as in Fig. A.1.

Fig. A.3 Phase diagrams for MT696. Designations are the same as in Fig. A.1.

Fig. A.4 Phase diagrams for MT720. Designations are the same as in Fig. A.1.

Fig. A.5 Phase diagrams for Schulte 5. Designations are the same as in Fig. A.1.

Fig. A.6 O−C diagram for Schulte 5. Points denote the time difference between the observed times of minimum and the calculated times of minimum using the Pemp and Ṗemp from Table 2. Dashed line shows the O−C value of zero.

Fig. A.7 Phase diagrams for Schulte 3. Designations are the same as in Fig. A.1.

Fig. A.8 Phase diagrams for MT059. Designations are the same as in Fig. A.1.

Acknowledgments

This work was supported by institutional research funding IUT40-1 and IUT40-2 of the Estonian Ministry of Education and Research, by the Estonian Research Council grants PUTJD5 and 8906. The authors thank Taavi Tuvi for his involvement in finding the literature data and for the preliminary analysis for finding the period change values. All the figures were made using the gnuplot plotting utility. In this work the SExtractor source extraction software was used (Bertin & Arnouts 1996).

References

All Tables

All Figures

	Fig. 1 Full frame of the Cygnus OB2 region from our first data set (iTelescope) with the selected binaries shown in red. The dashed line indicates the FOV of our second data set (RAITS).
In the text

Fig. 2 Available data shown in time domain with reduced Julian dates at the bottom and calendar years at the top. Point sizes are indicating the number of measurements in a one-day bin. We have changed the scale of the plot in the middle (dashed line) to better show the data available from 1990 onwards. Green points indicate our photometric data sets, red points are the available RV data, purple points show the times of minimum from Sazonov (1961), and black points are all other available photometric data listed in Table 1.

In the text

	Fig. 3 V-filter phase diagram of Schulte 5, folded with a constant period. The iTelescope data (black points) and Hall (1974) data (blue points) have been obtained 45 years apart. Solid lines are smoothed empirical light curves. The visible phase shift is due to the cumulative effect of period change (Ṗ).
In the text

	Fig. 4 Parameter estimates for period and period change. Periods are given as a difference from the Pemp value (see Table 2). The coloured solid line represents the results for an empirical phase model, and the dashed line represents the results for the calculated PHOEBE phase model. The inner contour is the 68% and outer contour is the 95% confidence level.
In the text

	Fig. 5 Binary modelling for Schulte 5 with the PHOEBE software. The upper panel shows our photometric iTelescope data (data set 1) with the model curve (solid line), and the lower panel shows combined radial-velocity data from Bohannan & Conti (1976) and Rauw et al. (1999) with the model curves for both primary and secondary components.
In the text

	Fig. 6 Binary modelling for A36 with the PHOEBE software. Designations of lines, labels, and photometric data set are the same as in Fig. 5. Radial-velocity data are taken from Kiminki et al. (2009).
In the text

	Fig. 7 Binary modelling for MT059 with the PHOEBE software. Designations of lines, labels, and photometric data set are the same as in Fig. 5. Radial-velocity data are taken from Kiminki et al. (2008).
In the text

	Fig. 8 Binary modelling for MT696 with the PHOEBE software. Designations of lines, labels, and photometric data set are the same as in Fig. 5. Radial-velocity data are taken from Kiminki et al. (2012).
In the text

	Fig. 9 Binary modelling for MT720 with the PHOEBE software. Designations of lines, labels, and photometric data set are the same as in Fig. 5. Radial-velocity data are taken from Kiminki et al. (2008).
In the text

	Fig. 10 Binary modelling for Schulte 3 with the PHOEBE software. Designations of lines, labels, and photometric data set are the same as in Fig. 5. Radial-velocity data are taken from Kiminki et al. (2008).
In the text

	Fig. A.1 Phase diagrams for A36 using the Pemp and Ṗemp from Table 2. Blue dashed lines in the magnitude plots denote a model light curve from PHOEBE software, and the red solid lines denote an empirical light curve. The solid line in the RV plot represents the model for the primary component, and the red dashed line represents the model for the secondary component.
In the text

	Fig. A.2 Phase diagrams for B17. Designations are the same as in Fig. A.1.
In the text

	Fig. A.3 Phase diagrams for MT696. Designations are the same as in Fig. A.1.
In the text

	Fig. A.4 Phase diagrams for MT720. Designations are the same as in Fig. A.1.
In the text

	Fig. A.5 Phase diagrams for Schulte 5. Designations are the same as in Fig. A.1.
In the text

	Fig. A.6 O−C diagram for Schulte 5. Points denote the time difference between the observed times of minimum and the calculated times of minimum using the Pemp and Ṗemp from Table 2. Dashed line shows the O−C value of zero.
In the text

	Fig. A.7 Phase diagrams for Schulte 3. Designations are the same as in Fig. A.1.
In the text

	Fig. A.8 Phase diagrams for MT059. Designations are the same as in Fig. A.1.
In the text