© ESO, 2017

1. Introduction

The detailed stellar evolution theory of massive stars of over 8 M_⊙ is currently an open question, as physical phenomena like magnetic fields, core convective overshooting, internal stellar rotation and angular momentum distribution inside stars are poorly understood. There have been parametrised descriptions of those phenomena in corresponding models but they remain uncalibrated with high-precision observational data (Aerts et al. 2013). Concerning photometric monitoring, the space missions MOST, CoRoT and Kepler have observed hundreds of OB dwarfs during the last decade and the pulsational modelling of a few tens of objects has greatly contributed to the asteroseismic description of their internal structure. But this is still an exploratory phase, and observations for ensembles of massive stars with known metallicity and age are needed (Aerts 2015). The case of evolved OB stars is even more challenging, mainly due to the very small number of objects monitored from space up to now. An additional obstacle is the relatively short campaign lengths ranging from weeks (MOST) to months (CoRoT; Kepler 2; BRITE) which is not sufficient to investigate longer oscillation periods reaching up to several months in the case of OB supergiants (e.g. Moravveji et al. 2012). The long-term high-duty-cycle monitoring is also useful for main-sequence variables (β Cep and SPB stars) to check the possible long-term frequency or amplitude variability.

As a rule, ground-based photometric observations can not achieve the required precision of 0.1 to 0.5 mmag that is needed for proper asteroseismic analysis (Aerts et al. 2013). However, ground-based multicolour light curves (giving a precision of a few mmag) can be used to detect the pulsation modes with suitably high amplitudes (Saesen et al. 2013). For the less-studied supergiants, this approach would be especially profitable.

The spatial concentration of massive stars is highest in young open clusters and OB associations. For example, ~70% of Galactic O-type stars are found in young clusters or loose OB associations (Portegies Zwart et al. 2010), making these systems ideal targets for massive-star studies. The scarcity of open clusters and associations older than 1 Gyr is attributed to encounters with giant molecular clouds, supported by the observed median cluster age in the solar neighbourhood of around 250 Myr (Kharchenko et al. 2005). This cluster age segregation means that different clusters contain massive stars in various evolutionary phases. In this way, the targeted photometric monitoring of clusters and associations to get a high yield of massive stars will also provide an opportunity to study the possible relationship between the age and variability characteristics of selected objects. In addition, the fraction of massive binary stars in bound systems can reach as high as ~ 90±10% when corrected for observational bias (Kiminki & Kobulnicky 2012). Although photometric monitoring campaigns typically only detect binarity for system geometries that are close to edge on, each new eclipsing binary discovered among massive stars helps to constrain the binarity fraction.

It is not feasible to observe bright objects in large sky areas (e.g. clusters and associations) with high cadence using large-aperture telescopes because of their low availability and typically limited field of view. But a number of small telescopes dedicated to astrophotographers exist that are suitable for bright-star observations with relatively large field of view and photometric filters. Commercially available telescopes (e.g. iTelescope.net or slooh.com) give us access to astronomical sites with better astroclimates and minimal effort of maintaining the equipment.

The use of automated small telescopes has resulted in many variable-star all-sky surveys such as the All Sky Automated Survey (ASAS, Pojmanski 1997), the Northern Sky Variability Survey (NSVS, Woźniak et al. 2004), and WASP (Butters et al. 2010) amongst others. These surveys have revealed numerous new discoveries of variable stars (Pojmanski et al. 2005; Norton et al. 2007; Nedoroščík et al. 2015) from eclipsing binaries to stars with as yet undefined variability types. But the need to cover nearly all of the observable sky takes its toll on the time resolution of the light curves for a given star, and therefore smaller-scope campaigns are needed to improve the observational time coverage for specific fields. In addition, the brightest stars in the sky are often overexposed in those surveys and consequently understudied (bright limit for ASAS and NSVS is V ≈ 8 mag).

In order to study massive stars in different evolutionary phases, we photometrically observed a number of open clusters and OB associations in the northern sky over the 2011–2013 period with high cadence. We used a small-aperture telescope with exposure times optimised for the brightest objects in those fields. Our basic goal is to investigate selected variable stars case-by-case and to incorporate additional spectroscopic data with appropriate spectral and temporal resolution (e.g. Laur et al. 2015; Kolka 2013). The goal of this paper is to expand the variable-star knowledge in our selected fields by extracting the variability types from the light curve information. Using a discrete Fourier transform method and time series feature extraction tool, we have analysed over 20 000 stars in our fields, and found 354 variable stars. We describe the released catalogue and give an overview of the variable stars found in the data.

2. Observations and data reduction

We observed 22 northern open clusters and associations in four different constellations across 15 fields. The fields were chosen for their number of bright OB stars and the age difference of the clusters. The observations were not centred on the clusters but instead the field of view was set to maximise the number of massive stars in each field based on the catalogue of supergiants and O stars in associations and clusters (Humphreys 1978).

The data set was obtained between 2011 and 2013 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 from Custom Scientific was used, providing a field of view of 60.5×40.8 arcmin with pixel scale of 1.64 arcsec.

The time resolution was one of the priorities of the campaign and therefore observations were carried out as frequently as possible. During the campaign, the telescope was used on 488 nights out of 996, having an uptime of 49%. We observed every available field once per night with occasional nights dedicated to specific fields with multiple observations.

A single observation consisted of three consecutive exposures in a set of VI or BVI passbands with exposure times optimised for the brightest sources in each field. For some fields, we opted for two different exposure times in order not to saturate the brightest stars and, at the same time, to monitor some of the fainter massive stars in that field. Analysis of our data set was mostly done using the V passband. B passband exposures were only obtained for two of the fields, so we can not use it for consistency reasons, whereas I passband exposures were optimised for the red supergiants, leaving most of the hotter sources underexposed and thus having larger scatter than in V passband.

Table 1 lists all of the observed fields, the number of nights they were monitored, and open clusters within those fields. The ages of the clusters have been gathered from literature with the references given in parentheses. In addition to smaller works, two larger catalogues were used for the cluster age determination. The first is by Kharchenko et al. (2013) who compiled a catalogue of over 470 million stars and use a homogeneous method for determining cluster parameters for 3006 open clusters which are published as the Milky Way Star Clusters Catalogue (MWSC). The other compilation of cluster data is the DAML02 (Dias et al. 2002) database, containing cluster parameters for more than 2167 open clusters. On a side note, the age discrepancy of NGC 7654 is thought to be caused by multiple and independent star-formation periods as there seems to be an older cluster behind a group of younger stars (Pandey et al. 2001).

Data reduction was carried out using the PHOTWORK software (Tuvikene 2012) based on the IDL Astronomy User’s Library¹ and the magnitudes were extracted by means of aperture photometry. We used an aperture size that was scaled by the full width at half maximum (FWHM) of stars in any given frame with the typical aperture radius of 1.8×FWHM. For further analysis, we only used instrumental magnitudes with uncertainties less than 0.1 mag.

Additionally, SExtractor source extraction software (Bertin & Arnouts 1996) was used for PSF photometry. This was necessary to disentangle the light curves of stars where the signal of a source in a given aperture was contaminated by neighbouring stars. We determined this contamination from FWHM-magnitude plots by a detection of a strong correlation between the FWHM and magnitude values. This correlation indicates a nearby star that affects the total extracted magnitude when the aperture size changes with the FWHM value. To disentangle these contaminated sources, PSF photometry was used.

We used Bouguer’s extinction plots to flag data points affected by bad weather (Sterken & Manfroid 1992). We flagged data points for which the brightness drop due to extra atmospheric extinction exceeded 0.8 mag from the fitted magnitude to airmass slope. This indicates that the bigger than average extinction in the atmosphere is due to clouds.

Flat-field frames were supplied by the service provider for between one and ten nights per month, with three to six flat frames per night. Master flats were constructed individually for each observing night using raw flats from the ten closest nights. We excluded frames that exhibited systematic differences over 5% from the median flat of neighbouring nights. Typically between six and twelve frames from between two and four nights were combined to an individual master flat. During periods with very few provided flats, we constructed master flats from frames of a single night with minimal time separation from the observations.

Scattered light reaching the focal plane of the camera causes illumination gradients in flat-field frames and using these flat fields in the data reduction may yield position-dependent systematic uncertainties in the resulting photometry (Manfroid 1995). The T4 telescope had a German equatorial mount and due to the aforementioned scattered light problem, the data show systematic effects depending on the telescope being on the east or west side of the pier. To correct for these effects we constructed correction surfaces for the full frame. This was done by dividing the data of constant stars into distinct time intervals based on the systematics in the light curves and then measuring the offset of mean magnitudes relative to the first time interval. In this way we found a global difference surface over the whole frame for each of these time intervals in every filter. In total, we constructed correction surfaces for 14 time intervals that were separated based on the combination from the pierside position and other systematics in the data. We generated two different types of correction surfaces using either a median smoothing or a polynomial fitting. The analysis of light curves were done using a combination of both corrections.

The amount of flat fielding error increases with the distance from the centre of the frame. So in order to minimise correction uncertainties from comparison stars, they were selected as close to the frame centre as possible. The correction was applied to both the programme star and the comparison stars; the mean correction value is 0.017 mag. The nightly mean precision of brightest stars is better than the mean correction value (cf. Fig. 5) and as a result, the total signal-to-noise ratio of the data can be improved by implementing this correction.

A separate correction was carried out between Julian dates (JD) 2 456 189 and 2 456 243 as the whole field was rotated about 30 degrees. Also, the period between JD 2 456 350 and 2 456 413 needed to be corrected as there were erroneous flat-field frames supplied.

An example correction of a light curve is given in Fig. 1. The raw data (upper panel) show offsets in the light curve (inset panels) that are corrected (lower panel) based on the surfaces combined from constant stars. The need for the correction is more evident in the phased light curves (big panels) where the scatter of a phase curve can be seen reduced after applying the correction.

Due to telescope control software issues, telescope centring on the target was also not consistent. These issues caused the frames to be at offset with respect to the targeted field. In the worst case, that offset was up to 1/4 of frame size, excluding about one third of the target stars. However, this problem affected only a limited number of exposures.

To analyse the data, we used multiple comparison stars that were chosen to be, preferably, close to the centre of the frame to minimise systematic errors across the frame. They were also chosen to be usable in both theV and I passbands. The final comparison stars were selected by studying star pairs that matched the aforementioned criteria and provided minimal scatter in the combined light curve. We used between two and five comparison stars per field. This combination of comparison stars gives us a relative photometric precision of better than 0.01 mag.

In addition to the photometric observations, the classification of selected brighter variable stars without known spectral types was carried out using the 1.5 m AZT-12 telescope at Tartu Observatory, Estonia. Spectral data were collected using a long slit spectrograph ASP-32 with 600 lines/mm grating, yielding R ≈ 1550 and covering wavelength range λ ≈ 3680 to 5830 Å.

All of the collected spectra were reduced using IRAF⁶ tasks. First, a zero-level correction was applied along with the pixel-to-pixel non-uniformity corrections using incandescent lamp flat-field data. Then, for the wavelength calibration, a dispersion relation was established using a ThAr spectral lamp. Finally, to correct for the instrumental sensitivity effects, target spectra were corrected using one or more spectrophotometric standard stars from CALSPEC⁷ database (Bohlin 2014) that were observed on the same nights as the targets. For classification purposes, a digital spectral classification atlas⁸ by R. O. Gray for late-type stars and synthetic spectra for B-type stars⁹ (Gummersbach & Kaufer 1996) were used.

3. Calibrated data

In total, 28 082 CCD frames were obtained and calibrated, of which 1942 were obtained in B, 13 428 in the V passband and 12 712 in the I passband.

We give the airmass histogram in the upper panel of Fig. 2. As the observation plans for the campaign were sequenced to observe every field preferably at lower airmasses, the airmass values are mainly in the range of 1 to 1.5. The two peaks of airmass values are due to two distinct maximum altitude ranges for our selected fields as seen from the New Mexico site.

The FWHM values on the lower panel in Fig. 2 were obtained from the CCD frames by analysing the brightest stars. These values varied from 3 to 14 arcsec with the median value of 5.2 arcsec, affected by telescope optics. Our obtained FWHM was good for observing bright stars, as we could use longer exposure times without overexposing our target stars. Although larger FWHM lowers the signal-to-noise ratio of faint stars and affects our limiting magnitude, the focus was on the brightest cluster members and the trade-off was deemed acceptable. We limited our analysis to using FWHM values of less than 10 arcsec, as the data become much worse at high values.

We obtained light curves for 25 395 sources. Figure 3 shows the V magnitude distribution of our observed stars. The measured stellar brightnesses were not transformed into standard photometric system. Instead we used an instrumental system that was shifted close to Johnson-Cousins system by using known standard magnitudes from the APASS¹⁰ DR9 catalogue. The fall-off in the number of stars with the apparent magnitudes fainter than 15 mag is due to the short exposure times.

The V passband data are summarised in Table 2. We list the V passband exposure times in Col. 2 (Exp.) and the total number of detected objects in Col. 3 (N). The given V magnitude range in Col. 4 is where the detected objects fall into. The limiting V magnitude (Lim.) values in Col. 5 are given as the limit where fainter stars are detected on less than 90% of the frames taken for an individual observed field (cf. Fig. 4). Stars that are detected in over 90% of the frames are defined as stars with good photometric quality (N₉₀). These are stars for which the detection on our frames was not affected by the airmass or seeing values. The total number of N₉₀ stars is around 2900. The quick dropoff in detection for fainter stars is to be expected due to variable FWHM values and strict photometric detection limit.

Column 6 in Table 2 lists the typical photometric internal root-mean-square (rms) noise (σ) of the constant stars to quantify our photometric precision at the limiting V magnitude (Lim.). The rms noise is defined by Kjeldsen & Frandsen (1992) and is calculated from the consecutive magnitudes measured during the same night. We give the precision as an intersection on a σ-magnitude diagram between the internal noise and the limiting V magnitude value.

A representative σ-magnitude diagram of Cygnus OB2 field is shown in Fig. 5. The resulting precision at limiting magnitude was obtained as a median internal noise from stars in a 0.4 mag bin around the limiting V magnitude value (between the two dashed vertical lines).

4. Detection of variability

Variable stars were initially selected on the basis of a variability index. We define the variability index as a ratio between the standard deviation of the full light curve and the internal rms noise (Kjeldsen & Frandsen 1992). For a constant star, the value would be close to unity whereas stars with values above 1.5 are counted as photometrically variable objects. The variability index is useful for a wide variety of observations, as the internal rms noise takes into account all instrumental uncertainties and the resulting value is therefore independent of the total scatter of the data, effectively measuring only the signal-to-noise ratio of variability.

Throughout the variability determination, we use measurements in the V passband because of the better signal-to-noise ratio for most of the sources. As the variability index is a measure of the signal-to-noise ratio of variability, we calculated the value of the variability index for most of our objects. We only omitted stars that had fewer than 100 measurements from the selection to exclude very faint stars and spurious detections. For every star, we calculated the variability index based on three different data sets. Two of these data sets were acquired by applying the two flat-field correction surfaces described in Sect. 2 and the third was the raw light curve. The final selection was based on the lowest variability index value of those data sets to minimise erroneous variability detections.

We then flagged our selected variable stars, based on the presence or absence of systematic effects, under different flags. As the aforementioned problems requiring the flat-field correction were too severe in some cases, we had to flag a number of stars unusable when a clear separation between the data from the east and west side of the pier was seen. Additional flagging was done based on the photometry extraction method used. We used aperture photometry for brighter sources and PSF photometry for fainter and crowded stars.

To analyse the light curve shapes of different type of variable stars, we used a Python library called Feature Analysis for Time Series (FATS). It extracts up to 64 different features from an unequally sampled time series that are described in Nun et al. (2015). We used some of these features as supplementary information to describe our data set.

5. Results

In total we found 354 variable stars in the observed fields. We list objects in tables separated by variability type in Appendix A with the full catalogue available at the CDS.

Table 3 is an excerpt from the full table (Table A.1 at CDS) of all the eclipsing binaries we found. The table lists our field names and internal ID-s of the stars. This was added to help identify the stars in our finding charts that are added in Appendix B. The designations are based on the SIMBAD database and, in case the star was absent there, we used The Two Micron All-Sky Survey (2MASS) designations. Given coordinates are based on 2MASS All-Sky Point Source Catalog¹¹.

The spectral types listed in Table 3 are combined from SIMBAD database, the Tycho-2 spectral type catalogues (Fabricius et al. 2002; Wright et al. 2003), observations by Currie et al. (2010), and spectral classification compilations by Skiff (2014) and Wright et al. (2015). Spectral types with lower case letters are obtained from photometric observations in the Vilnius 7-colour system (Milašius et al. 2013).

Based on our spectral observations from Tartu Observatory, we obtained spectral classifications for ten variable stars of which six had no previous classifications. These classified variable stars are listed in Table 4 and marked with an asterisk (*) behind the classification in Table 3.

From our spectra, BD+60 358 is a B1Ve type star as both Hβ and Fe lines are in emission. V472 Cas is an M4 or later type star and based on the Ca I 4227 line, its luminosity class is from V to III. V2538 Cyg is a previously known eclipsing binary star (Nicholson & Scottish 2006) and the lines from its companion are also visible in our spectra, making it an SB2-type spectroscopic binary with a B1III-V type primary and a B-type secondary component (cf. Fig. 6). The spectra for V2538 Cyg was taken at a phase of 0.225 and the component radial velocities of −207 and 226 km s^-1 were estimated by Gaussian fitting the line components.

To calculate the cluster membership fraction among the variable stars, we used the MWSC project (Kharchenko et al. 2013) which incorporates cluster membership probabilities around open clusters. They give three different probabilities based on proper motion and photometric J–K and J–H colour. These probabilities were determined from the location of the stars with respect to the reference sequences of either the average cluster proper motion in kinematic diagrams or the isochrones in photometric diagrams. They find about 400 000 stars to be cluster members having kinematic and photometric membership probabilities higher than 60%. We adopt their membership probabilities (kinematic probability (Pkin) >50% and photometric probabilities (PJKs+PJH) >100%) to determine the cluster membership of our variable stars.

In Table 3 we compare our variables with previously known variability information on these objects. For this comparison we use the International Variable Star Index (VSX)¹² database run by the American Association of Variable Star Observers (AAVSO) to identify possible first discoveries of variability. Due to the extensive compilation of different sources, the VSX database includes information on 398 030 (as of July 2016) variable stars. In the VSX “info” and “period” columns, we label stars that have no entry in the VSX database with “–” and those matches in the VSX catalogue with no variability classification with “?”. Table A.11 is a compilation of all the variable stars that are not listed in the VSX database.

In Table 3, we give an amplitude value as a modified FATS amplitude. The amplitude value in the table is defined as the difference between the median of the maximum 5% and the median of the minimum 5% magnitudes.

To search for periodic variables and eclipsing binaries from our data set, we used the discrete Fourier transform code SigSpec (Reegen 2007). It calculates a significance value which describes the probability that an amplitude is not caused by noise. SigSpec derives the resulting frequencies by iteratively prewhitening the frequency power spectrum using the most significant signal components from each run. The dominant frequencies were confirmed by visually identifying an amplitude of at least twice of the folded phase curve scatter.

We used the significance of any given frequency to determine the probability of the signal being either from a pulsating star, an eclipsing binary or a false signal due to the aliasing from our window function. We distinguish between a pulsational star and an eclipsing binary by visually inspecting the phase diagram. In the case of suspected eclipsing binaries, we also used the phase dispersion minimisation method (Stellingwerf 1978) to improve the period estimation as it accounts better for the different shapes of binary folded light curves than the Fourier method does.

In Table 3, we report only the periods corresponding to the dominant frequencies. We list β Cep and slowly pulsating B (SPB) stars with multiple frequencies in Table 6. A frequency power spectrum peak is commonly called significant when the signal-to-noise ratio of an amplitude is S/N> 4 (Breger et al. 1993) which corresponds to a significance value of 5.46 in SigSpec. Following the more conservative approach by Chapellier et al. (2011), we list only the frequencies with significance value over 20 because the periodograms below this began to show a forest of peaks.

Ellipsoidal variables (Morris 1985) are close binary stars that are non-spherical due to tidal interactions and exhibit sinusoidal phase curves due to the changes in the cross-sectional area that is visible to the observer. In addition, some cases of overcontact binaries exhibit similar sinusoidal phase curves. The phase curves of these binary stars are identical to pulsating variables and we only classified stars with light curves that show unequal minima or maxima as eclipsing binaries as this indicates tidal distortion effects or an uneven surface distribution of chemical elements (Kourniotis et al. 2014). Variables with no identified periods and unknown spectral types were left unclassified.

6. Classification of variable stars

Table 5 summarises our observed variable stars. The classification of variability was based on the period of variability as well as position on the HR diagram (Aerts et al. 2010). As we observed only the brighter objects in each field, the variability types we found are also from the upper part of the HR diagram, containing mostly α Cygni, β Cephei, Be, SPB, Cepheid, semi-regular (SR) and Mira type variables. The full tables listing all the members of each type are available in Appendix A.

6.1. Eclipsing binaries

We classified eclipsing binaries based on their previous classifications from VSX or, in case of new discoveries, based on their non-sinusoidal phase curves. In total, we found 80 eclipsing binary stars (Table A.1 at CDS) in all of the fields, for 69 of which we identified their orbital periods. The remaining 11 were named binary candidates as their light curves include some eclipse-like drops, but no period could be extracted due to the limited data. According to VSX, 47 eclipsing binaries are new discoveries with no known period and variability designation. Using our data set, we have previously analysed seven of the known eclipsing binaries in Cygnus OB2 association and found that A36, Schulte 5, [MT91] 059, and [MT91] 696 show a non-zero period change value at two sigma confidence (Laur et al. 2015).

6.2. β Cephei stars

βCephei stars are pulsating early-to-mid B-type dwarfs to giants with periods ranging between two and eight hours. The total number of known β Cephei stars is around 200 (Pigulski & Pojmański 2008).

We find 13 β Cephei candidates (Table A.2 at CDS) of which two have been previously listed in the VSX database and eight have a known spectral type. We report four β Cephei stars that have multiple frequencies over the significance value of 20 (Table 6).

The only β Cephei star in our selection with a previously determined period is V665 Per. Gomez-Forrellad (2000) reported two periods of 0.19493 and 0.24233 days. Their first period can be found in our data with significance value of eight. Their second period is closest to our dominant period of 0.242277 d but we would like to point out that our data do not support their period of 0.24233 d. This difference in periods can be interpreted by a period decrease of around five seconds in 15 yr (Ṗ/P = −0.24Myr^-1). This period change rate could be used to constrain evolutionary models for β Cephei stars by assuming different amounts of initial rotation and convective core overshoot and comparing the theoretical rate of period change to the measured one as done in Neilson & Ignace (2015).

6.3. SPB stars

SPBs consist of mid-to-late B-type dwarf to giant stars with periods roughly between 0.8 and 3 days (Aerts et al. 2010).

We find 16 possible SPBs (Table A.3 at CDS) based on their periods, of which 12 also have a spectral classification in the literature. We report three SPB stars with two frequencies over significance value of 20 (Table 6).

Our only matching SPB in the VSX database, V359 Per, is listed there as a β Cep star. The previous determination is done by the IOMC pipeline (Alfonso-Garzón et al. 2012) and they find a period of 0.275 d¹³. This period is not supported by visual inspection of the IOMC phase curve nor by our data. We find a different period of 1.2309 d, which is more in line with the star being an SPB type (cf. Fig. 7).

6.4. Be stars

While some Be stars show oscillations that are in accordance with either β Cep or SPB stars (Walker et al. 2005; Balona 1995), Be stars infrequently exhibit photometric outburst behaviour linked to the dispersion of the circumstellar disk (Rivinius et al. 2013).

As Be stars are classified from their spectra, we have compiled here every variable B star that has shown some kind of emission in the past. In total, we find 62 Be stars (Table A.8 at CDS) in our selection to be variable. Two of them, TYC 4279-2067-1 and V356 Per, exhibit additional β Cep-like pulsations with periods of 0.582 d and 0.342 d, respectively. V356 Per is previously classified by IOMC pipeline as a β Cep candidate¹⁴. They also find a period of 0.30595 d which is not supported by our data. Instead, we find a period of 0.342 d with a phase curve that is more sinusoidal than the IOMC data.

Be stars with photometric outbursts are indicated as “Y” in the outburst column in Table A.8 at CDS. Outbursts can be caused by a heating front that propagates across the disk and, after reaching the critical density, pushes the disk into a hot state with high mass transfer and high accretion luminosity. (Mennickent et al. 2002) distinguish four types of photometric behaviour with type 1 being caused by accretion of matter with lower mass transfer rates and type 2 by a sustained period of a high mass-transfer rate. 27 of our Be stars exhibit outbursts with sudden rises and declines (type 1) during the three seasons. The highest amplitude outburst of 0.87 mag was seen in EO Per during the 2012 season. In addition, VES 203 exhibits behaviour of sudden brightness jumps (type 2) where the brightness dropped 0.2 mag for 100 days in the beginning of 2012 season and then stayed at a lower level for the remainder of the campaign.

6.5. α Cygni stars

αCygni variables are non-radially pulsating supergiants of B–A spectral types with irregular variability caused by the superposition of many oscillations with close periods (Aerts et al. 2010; Saio et al. 2013). The typical observed timescales range from 1 to 50 days (Eyer & Mowlavi 2008).

We classified all of our variable B-to-early-A type supergiants (luminosity class I and II) as α Cygni variables. In total, we find 30 α Cygni variables (Table A.7 at CDS) and for three of them, HD 229059, HD 14143 and V757 Per, we find pulsation periods. HD 14143 and V757 Per have previous period determinations of 0.0892 (Koen & Eyer 2002) and 0.1643 (De Cat et al. 2007) days, accordingly, although we found different period values with 0.5321 and 0.2453 days. Altogether, VSX lists periods for six of our determined α Cygni variables although none of those periods can be seen in our data. We may not have observed the same periods in our data because the variability of α Cygni type stars may be cyclic in nature, caused by interior mixing processes and atmospheric instabilities (Aerts et al. 2010).

6.6. Wolf-Rayet stars

Almost all Wolf-Rayet (WR) stars show no strict periodicities and the origin of their variability remains poorly understood (Aerts et al. 2010).

According to the spectral classification, we find three variable WR stars (Table A.6 at CDS): WR 141, WR 158 and V444 Cyg. WR 141 and WR 158 exhibit irregular variability with amplitudes of 0.108 and 0.073 mag, respectively. WR 158 is previously known spectroscopically variable star (St-Louis et al. 2009). WR 141 is previously classified as an SB2 binary (Marchenko et al. 1998) with an O5V–III companion and an orbital period of 21.698 d. Only from the first season of our observations do we find some signal with a period of 21.7 d, whereas there is no signal from the other seasons. This might be due to the flickering nature of WR 141 (Khalack 1996). V444 Cyg is a previously known eclipsing binary (Cherepashchuk & Khaliullin 1973) with a WR component. Its orbital period from IOMC data is 4.213 d¹⁵ which can be confirmed with our data.

6.7. Late-type stars

Variable stars situated on the red side of the classical instability strip and with long periods (P > 80 days) are classified as Mira or SRa type variables if they are periodic, and SRb type if they possess alternating intervals of periodic and slow irregular changes in their light curves (Aerts et al. 2010; Eyer et al. 2012). F–K spectral type weak-lined giants and supergiants are classified as SRd type variables. As we monitored our fields for three seasons our time series are not long enough to comfortably determine very long periods and therefore we classified all of the late-type stars into one group.

There are 46 variable stars in our fields with late-type spectral classification. Based on our measured V–I colours and the shape of the light curves, we find additional 17 late-type candidates (Table A.9 at CDS) with V–I over 1.8 based on the effective temperature and colour tables from Currie et al. (2010). Due to the varying degree of reddening in young clusters, some of the candidates may be included because of the reddening difference between the comparison stars and target stars. We tried to account for that by calculating the differential reddening of the comparison stars by colour transforming the known APASS and 2MASS colours in the field.

From the total of 63 late-type variables, 30 are listed in the VSX database with 16 of them having a determined period. We find semi-regular periods for 19 stars only three of which are close to the previously determined periods. Some periods were obtained after de-trending the light curve with a wide median-filter to smooth out the slow irregular changes in semi-regular stars. As the periods are semi-regular and vary in time, the discrepancy between our newly determined and the previously determined periods is understandable.

6.8. Cepheids

Cepheids are radially pulsating F–K type giants and supergiants and their luminosity classes change roughly from III to I during the pulsation cycle (Aerts et al. 2010; Turner 2012). A typical light curve of a Cepheid is skewed and periodic, ranging from 1 to 50 days.

We find seven Cepheid variables (Table A.4 at CDS) of which five were previously known. V1319 Cyg has been previously classified as an SRa type, probably due to its longer period. On the other hand, the folded light curve shows a classical skewed Cepheid shape and the spectral classification of K0I, determined from our observations at Tartu Observatory, is in accordance with a Cepheid rather than SR type of star. Of the two previously unknown Cepheid candidates, NGC 6913-36 exhibits classical Cepheid phase curve shape and has a spectral class estimation of “f” (photometrically determined F type) based on Vilnius photometric system (Milašius et al. 2013). The second candidate, 2MASS J20173903+3804143, is added as a candidate because of its brightness, period length and phase curve shape, although the photometric scatter is too high to distinguish between a noisy sinusoidal shape and a classical Cepheid shape.

We note that for the Cepheid NSVS 5731720, only the last two seasons can be used to fold the light curve nicely with the period (Fig. 8). The first season has approximately half the amplitude of the latter two, and when folded with the same period of 3.20612 days, the V passband light curve exhibits a binary-like shape with two unequal maxima. There is also a phase shift of 0.6 between the first two seasons. We are unable to explain this behaviour with systematic effects.

6.9. Other variables

There is one previously known γ Doradus star V1132 Cas (Table A.5 at CDS) for which we confirm the period of 0.626 days (Wyrzykowski et al. 2002). As γ Dor stars have masses of 1.5−1.8 M_⊙ and are of luminosity type IV–V (Aerts et al. 2010), the apparent magnitude of 8.16 mag for V1132 Cas suggests that it is a nearby field star.

There are an additional 80 variable stars that we could not classify due to the lack of spectral or pulsational information (Table A.10 at CDS). V1156 Cas and TYC 3695-992-1 have previously determined periods but as the stars are too faint in our data we were unable to confirm them. We find periods for 13 of these unclassified objects (Table A.10 at CDS), ranging from 0.5 to 102 days.

6.10. Examples

We present some light curve examples in Fig. 7. Phase curves in panels a and b show the difference in scatter for one bright and one faint star for which we managed to find a period. V359 Per is a SPB star with V magnitude of 9.06 and shows a scatter of around 0.02 mag. 2MASS J23163318+6031477 is a fainter β Cep star with V magnitude of 15.188 and shows a scatter of around 0.15 mag. While it is below our limiting magnitude for the NGC 7510 field, the high amplitude of 0.489 mag is enough to detect the period. We did not find any previous variability information on this star from any database, probably due to its low brightness

Panels c and d in Fig. 7 are examples of irregular or semi-regular behaviour. V504 Per is a Be star with one of the most numerous outburst counts, showing five major outbursts during the observing campaign. OQ Cep is an M7 type star exhibiting SR variation of over four magnitudes during the second observing season, the highest amplitude variation across all our data. We also find a semi-regular period of 385 days, based on the first two seasons, which is comparable to the previous period determination of 370 days (Alfonso-Garzón et al. 2012).

We present the amplitude distribution of variable stars in Fig. 9. The effective temperatures corresponding to specific spectral types found in the literature are acquired from Tables 7–9 in Currie et al. (2010). In case there was no luminosity class determined, their temperature is given as their spectral subgroup mean temperature. The figure shows that the late-type stars have a larger variability amplitude compared to the rest of the selection. Be stars skew the amplitude of hotter stars towards higher values due to their outbursts but as their variability is caused by the circumstellar disk instead of their intrinsic pulsations, they are marked separately with open circles.

7. Discussion

To test the suitability of the matches to eclipsing binary and pulsating variable categories, we use the FATS skewness and kurtosis statistics. Eclipses in eclipsing binaries cause the magnitude distribution to have a tail towards the fainter magnitudes and therefore positive skewness values, whereas pulsating stars have a close to zero or negative skewness values due to a more even magnitude distribution. This separates the two groups of periodic variability, eclipsing binaries and pulsating stars, on the skewness-kurtosis plane.

Figure 10 shows the skewness-kurtosis plane for the eclipsing binaries and pulsating stars we detected in our fields. The pulsating stars clump up on the left hand of the diagram and the eclipsing binaries are in a long tail in the right-upper part of the diagram, with the two types having a separation around a skewness value of 0.5 as indicated by Graczyk & Eyer (2010). The collection of binary stars up to a skewness value of 1 consists of ellipsoidal variables that show sinusoidal light curves similar to other pulsating stars but are a little more skewed due to irregularities in the uneven maxima of two consecutive periods. The coloured circles represent our binary candidates (where no orbital periods were found) and they follow the skewness-kurtosis trend of the other binaries.

Despite finding variability in previously known objects, our photometric depth does not always allow us to find periodicities to compare with already published values. Saesen et al. (2013) monitored NGC 884 over three observing seasons and found periods for 131 stars in their field. Eight of them match our listed variables and of these, we determined a period for only one object, due to the limiting magnitude in the field of our data.

Omitting eclipsing binaries, 35 out of the 275 intrinsically variable stars were found to be probable cluster members, based on the cluster membership probabilities of individual objects from MWSC (Kharchenko et al. 2013). This makes the total cluster membership fraction for variable stars less than 13% (membership fractions for individual types are given in Table 5). Using a limiting magnitude of 14 mag, that is, close to the limiting magnitude of our observations, we find a total of 21 792 stars around our observed clusters in the MWSC database. Applying the same membership criteria as we did with our variable stars (kinematic probability (Pkin) >50% and photometric probabilities (PJKs+PJH) >100%), we find that 11% of all those stars are cluster members. Thus the cluster membership fraction among variable stars is comparable to all of the stars around clusters according to the MWSC.

One of our goals in monitoring stellar clusters of different ages was to study massive stars in various evolutionary phases. To find independent stellar ages, we need accurate membership information for these stars. To this end, we compared MWSC catalogue to other works to test the reliability of their cluster membership estimations. Wright et al. (2015) studied 167 probable cluster members in the Cygnus OB2 association. Matching their 167 cluster members with MWSC database, we find that out of 155 matched stars, only 51 are probable cluster members according to MWSC.

The DAML02 database (Dias et al. 2002) also includes membership probabilities for stars around clusters although their sample size is much smaller than in MWSC. For example, in Berkeley 87 open cluster DAML02 has 157 nearby stars whereas MWSC has 23567. Matching the two data sets, we find 141 common stars. For this sample, DAML02 gives a membership fraction of 72% (membership probability over 50%) whereas MWSC gives a fraction of 29%. This shows that MWSC underestimates the cluster membership fraction compared to other works. But as it is the only database that covers most of our observed stars, we were obliged to use it for its completeness.

In order to find an upper limit to the cluster membership fraction, we studied the spatial distribution of variable stars in the plane of the sky and measured how variable stars are positioned around stellar clusters. For each star, we calculated the relative distance to the centre of its closest cluster d′ = d/r_C, where d is the angular distance and r_C is the angular cluster radius from MWSC. This cluster radius is defined as the distance from the cluster centre to where the surface density of stars becomes equal to the average density of the surrounding field. Based on these distances, we calculated the probability densities of finding a star at a certain relative distance from the centre of any cluster. We used the bootstrapping technique to get the probability density of all stars with 99.7% (3σ) confidence limits.

The resulting probability density functions (PDF) can be seen in Fig. 11. Panel a represents the combined distribution of α Cyg, β Cep, Be, SPB, Cepheid, γ Doradus, Wolf-Rayet and late-type variables. This distribution shows that there are relatively more variable stars inside one unit of cluster radius than there are non-variable stars (represented by the coloured area). Panels b and c show the PDF for α Cygni and Be type stars, respectively. They both display the same signal as the combined distribution, albeit the 3σ confidence area is larger due to smaller sample size. The PDF on panel d, on the other hand, consists of eclipsing binaries and shows no deviation from the 3σ level. This means that eclipsing binaries have no positional preferences in relation to the clusters. As the variability of binaries is not intrinsic, we omitted them together with unclassified variables from panel a.

The distribution of all stars shows a maximum at around 0.5 r_C with the width of the excess around 1 r_C which is to be expected from the definition of the cluster radius. The maximum of our observed variable stars in the plane of the sky is at around 0.3 r_C, meaning that the stars in the core of the clusters are more likely variable than the stars at the edges of the clusters.

Integrating the PDF of variable stars to one cluster radius, we find that 35% of variable stars are located inside the clusters. As we studied the spatial distribution in the plane of the sky, this can only be the upper limit of the cluster membership fraction due to fore- and background object contamination.

These examples indicate that the cluster membership probabilities in MWSC might be too strict and thus the membership fraction of 13% amongst our variable stars is underestimated and can go as high as 35%. We did not use DAML02 database for cluster memberships, as their sample size is much smaller and did not include many of our variable stars. But luckily, the reliability of membership probabilities will improve with the arrival of Gaia’s more accurate distance, proper motion and radial velocity data as it measured even the brightest of sources in our fields (Gaia Collaboration 2016).

8. Conclusions

We monitored 22 northern open clusters and associations in 2011–2013 across 15 fields. Using a small-aperture telescope, we optimised our observations for the brightest stars in each field and obtained light curves for ~3000 stars which were detected on more than 90% of the frames. We used a discrete Fourier transform code called SigSpec to obtain periods for pulsating stars and phase dispersion minimum algorithm to fine-tune the periods for eclipsing binaries.

Based on the variability index method we found 354 variable stars in the fields. In addition to the pulsational information, we used a pulsation HR diagram to classify our variable stars. We classified 80 eclipsing binaries, 31 α Cyg, 13 β Cep, 62 Be, 16 SPB, 7 Cepheid, 1 γ Doradus, 3 WR and 63 late-type stars. We compared our results with the Variable Star Index (VSX) to identify possible first discoveries of variability amongst the detected variable stars. We found 47 eclipsing binaries, 10 α Cyg, 11 β Cep, 29 Be, 15 SPB, 2 Cepheid, 1 WR and 33 late-type stars that are absent in the VSX database. An additional 70 variable stars could not be classified due to the lack of spectral or pulsational information.

We find that 13% of the variable stars are cluster members based on the membership information in the MWSC. The MWSC membership estimations are lower compared to other membership estimations, but we had to use it as it is the only database that covers most of our observed stars. We analysed the probability density function of variable stars and found the upper limit of membership fraction to be 35%. We list all 354 variable stars in CDS divided into different variability types. In addition, the light curve data will also be released in the CDS.

---

Acknowledgments

This work was supported by institutional research funding of the Estonian Ministry of Education and Research (IUT26-2, IUT40-1, IUT40-2), by the Estonian Research Council grants 8906 and PUTJD5, and by the European Regional Development Fund (TK133). This research has made use of the International Variable Star Index (VSX) database, operated at AAVSO, Cambridge, Massachusetts, USA. Figures were made using the Gnuplot plotting utility and APLpy, an open-source plotting package for Python hosted at http://aplpy.github.com. The use of SIMBAD database and AAVSO Photometric All-Sky Survey (APASS), funded by the Robert Martin Ayers Sciences Fund. are acknowledged. We thank an anonymous referee for constructive suggestions that have helped to improve our paper.

References

Appendix A: Variable star catalogue

The full variable star catalogue (Tables A.1–A.10) is available at the CDS along with light curves. Here we only list a subset of the catalogue for stars that were not found in the VSX database.

Appendix B: Individual fields

We present finding charts for our variable star candidates across the 15 fields in Figs. B.1–B.15. These finding charts were obtained by median smoothing our observed V passband frames taken with the 0.25-m Takahashi Epsilon telescope, located in Mayhill, New Mexico, USA. All of the frames were with a fixed field of view of 60.5×40.8 arcmin with pixel scale of 1.64 arcsec. Variable stars are indicated with a red circle and our internal ID number above them. We have only marked the stars listed in our variable star catalogues in Appendix A and the IDs can be cross-referenced from the second column.

All Tables

All Figures