AUTOMATED VARIABLE STAR CLASSIFICATION USING THE NORTHERN SKY VARIABILITY SURVEY

The Robotic Optical Transient Search Experiment (ROTSE-I; Akerlof et al. 2000) was a project whose main goal was to detect the so-called orphan gamma-ray bursts (GRBs) and provide follow-up for high-energy missions with rapid optical observations. With nearly a year in baseline, observations nearly nightly, and all-sky coverage, the project is an excellent resource for finding variable stars. The data were compiled as the SkyDOT database and released as The Northern Sky Variability Survey (NSVS; Woźniak et al. 2004). The database contains about 14 million unique objects, and we have identified and classified 4659 objects into several variable star classes. This was achieved with period, light curve information, Fourier coefficients, and some manual inspection, to be discussed later.

The identification of pulsating stars is important for the calibration and determination of distance indicators. Early in the 20th century, Cepheids were discovered to have a period–luminosity (and later refined to a period–luminosity–color) relation that makes them useful as distance indicators (see Tammann et al. 2003, 2008 for examples of modern studies). Cepheids typically have periods longer than 1 day and light curves ranging from nearly sinusoidal to strongly asymmetric. Having a larger sample of Cepheids allows for refinement of the relations. However, in old stellar systems, there are no Cepheids for distance estimates, but RR Lyrae variables are present. RR Lyr variables are radial pulsators with periods usually between 0.2 and 1.0 day with light curves similar to Cepheids. RR Lyr variables can also be used as distance indicators as they show a metallicity–luminosity relation (Sandage & Tammann 2006). These stars are essential for use as standard candles in nearby galaxies as well as globular clusters in our own galaxy and also for understanding stellar pulsation and evolution.

The δ Scuti variables are another type of pulsating star in the instability strip. These are typically main-sequence or slightly evolved stars of spectral types of A or F with periods between roughly 0.5 and 7 hr. The δ Scuti stars are divided into two main types: the high amplitude delta Scutis (HADS) and the small amplitude δ Scutis. HADS have photometric variations of more than 0.25 mag while the small amplitude types can have variations on the order of millimagnitudes. The small amplitude class usually pulsates in multiple modes and has nonradial modes, creating a very chaotic light curve as beating occurs. The HADS are typically monotone radial pulsators and have more stable light curves (Laney et al. 2003). As discussed by Breger (2000) and others, mode identification of multimode δ Scuti stars is very difficult, but is very important to constrain stellar pulsation models. Additionally, HADS with stable light curves can be used as a short-period extension of the Cepheid period–luminosity–color relation (Laney et al. 2003).

W Ursae Majoris-type stars (W UMa) are close contact binary systems. They are among the most numerous variable stars in the sky. The stars in W UMa systems are main-sequence stars of spectral types from early A to early K and typically have periods between 0.25 and 1.2 days (Rucinski 2004). Recent theories suggest that W UMa binaries start out as detached systems with P ≈ 2 days and very different masses. The binary loses angular momentum due to magnetized stellar winds and eventually forms a common envelope with mass transfer (Stępień 2006; Eker et al. 2007). Despite having very different mass ratios ranging from q = 0.97 to q = 0.066 (Rucinski et al. 2001), the surface brightness and effective temperature are nearly identical over the entire contact envelope surface. Due to their unique properties, W UMa systems can be used similarly to pulsating stars by using a period–luminosity–color relation as a distance indicator among solar-type stars (Rucinski 2004).

Eclipsing binaries of Algol and β Lyrae type offer a unique opportunity to determine stellar parameters with a high degree of accuracy, including a direct measurement of the radius of each star if the period, inclination, and radial velocity of each star is known. Algol and β Lyrae eclipsing binaries are differentiated by constant variation in the light curve outside of eclipses for the β Lyrae systems to the nearly constant brightness outside of eclipses of the Algol-type systems. Refer to Hoffman et al. (2008) for more information on these eclipsing binaries.

The initial requirements for a variable star candidate in the NSVS were the same as those of Hoffman et al. (2008), requiring at least a 0.1 mag variation about the median magnitude and at least 30 observations. This immediately reduced the number of variable star candidates from 14 million in the NSVS to around 100,000 objects.

These candidates are then put through Fourier analysis, magnitude ratio, and phase dispersion minimization routines to determine and identify periodicities in the data. This step is again identical to that of Hoffman et al. However, since we will later be fitting a Fourier series to the data, we invoke an additional constraint on the data of having at least one data point in each 0.05 phase bin. This ensures that the light curve is well sampled at all phases. We compute the median magnitude value in each bin to produce an averaged smoothed light curve. Combined with the magnitude ratio used in Hoffman et al., we determine whether the candidate is a pulsating (one minimum per cycle) or eclipsing/ellipsoidal variable (two minima per cycle). The process involves finding the global minimum and then searching 0.5 in phase from the minimum for a local minimum. If the light curve has another local minimum at this point or is relatively flat (as an Algol with no secondary eclipse might be), it is flagged as an eclipsing/ellipsoidal variable. Otherwise, the candidate is flagged as a pulsating star and the period is halved for the final tabulation. For the Fourier decomposition to follow, candidates flagged as pulsating stars are phased to twice their real period so that all candidates can be directly compared.

The next step in the process is Fourier decomposition (Rucinski 1993). We use the Levenberg–Marquardt nonlinear least-squares fitting routine to fit the following equation (Rucinski 1993) to each light curve:

$$\begin{equation} l(\theta) = \sum _{n=0}^6(a_n \cos (2n\pi (\theta - \phi))), \end{equation} \tag{ 1 }$$

where θ is the phase and ϕ is the phase offset to compensate for the unknown time of minimum/maximum. While six orders are used for the fit, only the first four are used for classification. If the fit does not converge (i.e., low-order coefficients near zero or unrealistically large values), the routine is run again with slightly different initial conditions. If it fails to converge again, the light curve is phased at integer harmonics of the period. If the fit still does not coverage, the candidate light curve is discarded as this is indicative of a nonperiodic variable or a light curve with insufficient data.

Light curves with periods very close (within 0.001 days) to an integer harmonic of 1 day are discarded. This is usually an indication that the fast Fourier transform performed earlier failed to find significant periodicities in the data. However, there are indeed some real variable objects with such periods that are discarded, but due to the rarity of such objects and the large number of rejected light curves, manual inspection was not performed for these sources.

If a₄ has a positive value, the candidate is assigned an initial classification of RR Lyr if the period is less than 1.2 days and as a long-period/Cepheid candidate otherwise. Algol/β Lyrae candidates identified previously in Hoffman et al. (2008) are removed from the candidate list. RR Lyr variables identified in the NSVS by Kinemuchi et al. (2006) are flagged. If a₄ is negative, the light curve is analyzed to determine the depths of the primary and secondary minima and combined with period and Two Micron All Sky Survey (2MASS) color information for further classification.

The objects are filtered through the SIMBAD database to find any previous identifications. Since all the objects in the NSVS are relatively bright, almost all of them are expected to have at least been identified by their coordinates (e.g., the Tycho-2 catalog; Høg et al. 2000), though only a small fraction has been previously identified as variables or has been classified.

Objects with |a₄| < 0.02 are flagged for close manual inspection. These objects are in the Fourier coefficient region where they could be either RR Lyr or W UMa variables with periods under 1.2 days. The light curves of the 728 objects that met those criteria were manually inspected and classified by their light-curve shape, as the classification by Fourier coefficients is not reliable in this region. However, even after manual inspection, there are likely to be many RR Lyr type variables of subclass RRc misidentified as W UMa variables and vice versa. This is because their light curves can appear nearly identical.

The new objects classified as Algol/β Lyr variable are listed in Table 1. This table does not include 548 objects classified as such previously in Hoffman et al. (2008). Table 1 columns include the R.A., decl., NSVS object identification number, period, the J − H, H − K, and K magnitudes extracted from the 2MASS database, the mean ROTSE unfiltered magnitude, the amplitude of oscillations determined from the Fourier fit, any previous identification, any previous classification, and notes. The candidates that have a₄ values greater than the line defined by a₄ = a₂(0.125 − a₂) (Rucinski 1997) can be considered β Lyr candidates, as this is the equation defining the envelope for inner contact. Objects that have a₄ values less than the line are likely detached eclipsing binaries or Algols. Figure 1 is a plot of the a₂ and a₄ Fourier coefficients of all 4569 light curves. The curved line is the theoretical limit for inner contact. The detached/semidetached variable star candidates (red crosses) can be considered β Lyr candidates if they are below the curved line. Figure 2 shows some example light curves of objects in this category with two Algol and two β Lyr light curves. The main source of contamination in Table 1 are W UMa type variables. Sometimes the Fourier fit does not reproduce the bottom of the eclipses very well and the algorithm mistakenly identifies the eclipse depths as very different, and thus an Algol system. While all of the objects were manually inspected for obvious period phasing errors, there could possibly be some Algol systems which have little or no secondary eclipse that will be misphased to twice the true period. After manual inspection it was found that of the objects initially identified as Algol/β Lyr candidates, ∼5% were later classified as another class.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Objects classified as RR Lyr variables are listed in Table 2. These are the open diamonds in Figure 1 and generally have positive a₄ values. The vast majority of these objects are of the RRab subclass (asymmetric light curves). The columns of Table 2 are the same as that of Table 1. RR Lyr variables identified by Kinemuchi et al. (2006) are denoted in the last column. Figure 3 shows some example light curves of objects in Table 2. Bumps on the descending branch can sometimes be seen depending on the amplitude and photometric errors. Since there is little overlap in period between the major classes of pulsating stars (δ Scuti, δ Ceph, and RR Lyr), it is fairly easy to classify these objects. There is contamination in the shorter-period candidates, with the most likely contaminant being W UMa stars. There could also be some γ Dor type variables in the sample. As seen in Figure 4, top panel, there is some contamination peaking around 0.35 days, which is due to the W UMa sample (bottom panel). The solid line is the distribution from Samus & Durlevich (2009) and the gray histogram is from the present study. Kinemuchi et al. also had this problem (dotted line). Since W UMa and RRc objects have nearly identical light curves, it is difficult to distinguish them based on single-color photometry alone. Thus, objects in Table 2 with periods less than ∼0.4 days should be treated with suspicion. After manual inspection it was found that of the objects initially identified as RR Lyr candidates, <1% were later classified in another class, though it is clear that there is a significant percentage of the sample that are likely W UMa contaminants. Additionally, ∼56% of our RR Lyr candidates were previously identified by Kinemuchi et al. (2006).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3 lists the objects classified as Cepheids and long-period variables. These are denoted as filled blue triangles in Figure 1, and have little correlation in their a₂ and a₄ values. Due to photometric constraints, it is generally not possible to distinguish between classical and W Vir Cepheids, although it is certainly possible to make such a classification with high-signal/noise photometry and well-sampled light curves. The columns of Table 3 are the same as that of Table 1, except for the final column. The period determined by the Fourier analysis is unreliable when it approaches the baseline of the data. Thus, objects with listed periods of greater than ∼40 days should be treated with suspicion because of the observational window. Figure 5 demonstrates this problem at longer periods. Our sample has significantly more longer-period Cepheids than the GCVS (Samus & Durlevich 2009) would predict. The final column in Table 3 denotes such issues that arise for each object. Figure 6 shows four example light curves of Cepheid candidates from Table 3, and shows the great variety of Cepheid light curves in the database. Unfortunately, it is difficult to distinguish long-period Cepheid and Mira type variables even with 2MASS colors as the color difference is not always statistically significant. Thus, there are some Mira variables that contaminate the Cepheid sample. There is also the possibility of long-period RR Lyr variables in the sample. This possible RR Lyr contamination can be seen by the slight deviation from the GCVS sample near periods of 1 day in Figure 5, though the rest of the period distribution matches quite nicely. As with Vilardell et al. (2007), there is a hint of a secondary peak around 15 days. Less than 1% of these candidates were later classified as another variable star type, but a much higher fraction have questionable periods, as discussed above.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Objects classified as W UMa type variables are listed in Table 4. This is the most common classification of the variables, which is expected since ∼95% of all variable stars in the solar neighborhood are W UMa type (Eggen 1967). The W UMa candidates are denoted as filled black circles in Figure 1 and generally fall below the curved line of theoretical inner contact. Figure 7 shows some example W UMa light curves for objects listed in Table 4. The most likely contaminant is RR Lyr variables of the RRc subtype. Without color information and higher precision photometry, W UMa and RRc variables cannot be easily differentiated. There can also be some β Lyr objects in this table if the Fourier fit failed to reproduce the eclipse depths accurately. There can also be ellipsoidal variables such as RS CVn stars in the table that are difficult to distinguish with the present data set. The bottom panel of Figure 4 shows a comparison of the period distribution between the NSVS sample and the W UMa objects in the GCVS (Samus & Durlevich 2009). The distributions match quite well considering our sample consisted of almost nine times more candidates than listed in the GCVS. The slight oversampling at the peak is likely due to RRc contamination. Approximately 5% of these W UMa candidates were later reclassified as Algol/β Lyr candidates. An unknown, but relatively significant fraction is likely RRc variables.

Zoom In Zoom Out Reset image size

The final table is Table 5, which lists the short-period variables. These are objects with periods less than 0.2 days and are denoted as green stars in Figure 1. The objects in this list are likely dominated by δ Scuti variables, though there could be a few β Ceph variables (pulsating early B-type giants/subgiants) in the list. The δ Scuti stars in this list are necessarily monoperiodic as any significant beating of modes would have blurred the light curve over the long time period of the NSVS and thus would have been rejected by our selection criteria. Figure 8 shows example NSVS δ Scuti light curves. Multiperiodic pulsating stars could be detected by multiple Fourier transform iterations, but because of the low-amplitude nature of these objects and the relatively poor photometric precision of the NSVS data, this was not attempted. Less than 1% of the candidates were later reclassified after manual inspection.

Zoom In Zoom Out Reset image size

During the manual inspection process, some light curves were noted as unusual in their characteristics. Figures 9–12 are examples of such objects. Figure 9 is EQ CMa, previously identified as an Algol binary. The light curve is quite unusual for an Algol system, with unequal light at phases 0.25 and 0.75. This could be due to a spotted surface, or that the system is actually a cataclysmic variable (CV) and we are seeing the hotspot being eclipsed. NSVS 5455133, whose light curve is plotted in Figure 10, was previously identified as an eclipsing binary in SIMBAD. The light curve shows a well-defined primary eclipse, but the secondary eclipse, which appears to have nearly the same depth as the primary is very poorly defined. Phasing errors are unlikely as the cause because only the secondary is not well defined and the data are spread over many months. A variable spotted surface on the secondary star could be the cause. Figure 11 shows another unusual light curve, NSVS 4491084, an object with no prior identification. There is a hint of a secondary eclipse at phase 0.5. Like EQ CMa, this could be a spotted system or a CV, though the odd saw-tooth shaped light curve outside of eclipse is unusual. Finally, Figure 12 shows a previously unclassified object, NSVS 4027822, with unequal light at phases 0.25 and 0.75. This is most likely a heavily spotted W UMa or spotted ellipsoidal RS CVn system because of the equal minima. All these systems require follow-up observations to determine their true nature.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Some objects with unusual light curves can be classified incorrectly with the techniques we have used. An example of such an object that was rejected after manual inspection is NSVS 6300521 (RX J0019.8+2156, QR And). Figure 13 shows the NSVS light curve of the object. Our routine classified the asymmetrical saw-toothed light curve as an RR Lyr star due to the ∼15 hr period of the object, but the object is actually a supersoft X-ray binary system with an unusual light curve. While such a situation is rare, there will be objects in our tables whose true nature is nothing like the variable class that we have assigned.

Zoom In Zoom Out Reset image size

Using a combination of Fourier coefficients, light-curve properties, period, and color information, 4659 variable objects in the NSVS were identified and classified into one of the five categories. The method greatly improves upon Hoffman et al. (2008) as evidenced by nearly twice as many Algol/β Lyr candidates found by this method. The main improvement is the fact that Fourier decomposition allows for light curves with only partial coverage of an eclipse to be useful. The method also almost doubled the number of identified RR Lyr variables in the NSVS database as found by Kinemuchi et al. (2006). If the classification of the W UMa and Cepheid variables are confirmed, they will also greatly increase the number of such objects known. With these variable stars, we have now exhausted the variable objects in the NSVS using our techniques. While more could potentially be found by relaxing the initial selection criteria, the large increase of false-positives would mean that more manual inspection would be required. A more productive use of the method is to apply it to other all-sky survey projects currently coming online (e.g., Pan-STARRS, LSST, etc.). The higher precision photometry and multicolor data of such surveys make them ideal for our pipeline to be applied as the misidentifications and false-positives would be reduced dramatically over that of the NSVS.

The total misidentification rate, where a candidate was later moved to another variable star class after manual inspection, is ∼4%. However, the true failure rate is higher, mainly due to the W UMa/RRc degeneracy problem, errors in the manual inspection, and the quality of the data. Approximately 15% of the total candidates were discarded after manual inspection due to poor data quality, or the obvious nonperiodicity of the data resulting from Fourier fitting errors. We expect that this fraction can be dramatically reduced by further refinement of the Fourier peak selection criteria outlined in Hoffman et al. (2008). Lower photometric errors would also reduce the discarded and misidentification rates, and lower-amplitude multimode systems could be identified.

This publication makes use of the data from the Northern Sky Variability Survey created jointly by the Los Alamos National Laboratory and University of Michigan. The NSVS was funded by the Department of Energy, the National Aeronautics and Space Administration, and the National Science Foundation.