AN OPTIMIZED METHOD TO IDENTIFY RR Lyrae STARS IN THE SDSS×Pan-STARRS1 OVERLAPPING AREA USING A BAYESIAN GENERATIVE TECHNIQUE

The SDSS (Stoughton et al. 2002; Abazajian et al. 2009) is a deep spectroscopic and photometric survey (g < 23.3) that uses five filters (u, g, r, i, and z) to survey ∼12,000 deg² of the sky. Although most of the SDSS data are based on single-epoch observations, ∼270 deg² of the Southern Galactic hemisphere, the so-called Stripe 82, have been observed around 80 times.

The PS1 3π survey (Kaiser et al. 2002, 2010) is a ∼3.5 yr (2010 May–2014 March) multi-epoch photometric and astrometric survey that is being conducted in Hawaii. The PS1 telescope repeatedly observes the entire sky north of declination 30° (3π survey). It uses a 1.8 m telescope with a 7 deg² field of view. It is equipped with the largest digital camera in the world (1.4 gigapixels). One of its goals is to carry out a photometric and astrometric survey of stars in the Milky Way and the Local Group in five bandpasses (g_P1, r_P1, i_P1, z_P1, and y_P1) covering the spectral range of 4000 Å < λ < 10500 Å. More information about these filters can be found in Tonry et al. (2012). The PS1 obtains multiple images of three quarters of the celestial sphere in the optical and near-infrared (Kaiser et al. 2002) to ∼22 mag in g_P1 in individual exposures (Morganson et al. 2012). Specifically, it is designed to take four exposures per year and area with each of its filters (Morganson et al. 2012). By the end of the survey there should be ∼12 exposures per field and filter. Currently, the average number of observations in each of the five filters is 8 (Magnier et al. 2013).

The PS1 was mainly designed to detect potentially hazardous asteroids and near Earth objects (Kaiser et al. 2002). Because it is a deep survey that is repeatedly observing three quarters of the sky, its data are of interest for a wide range of different scientific topics. These topics cover different science areas, from solar system objects to cosmology. The PS1 data are of particular interest also for structural studies of the Milky Way affording a deeper and wider area coverage than previous surveys. When more than ∼4 epochs are available in at least two filters (∼4 epochs in each filter), the repeat observations of the PS1 allow one to identify variable stars such as RRL stars.

The SDSS colors of RRL stars have been studied and characterized using the 483 RRL stars detected in Stripe 82 (Sesar et al. 2010, and other studies). Since the SDSS (u − g) color serves as a surface gravity indicator for these stars, the range (∼0.3 mag) and the root-mean-square (rms) scatter (∼0.06 mag) are the smallest in this color (Ivezić et al. 2005).

The g, r, i, and z bands from the SDSS are similar to the g_P1, r_P1, i_P1, and z_P1 bands from the PS1, respectively. However, the u band is used only in the SDSS but not in the PS1, and the y_P1 band is found only in the PS1 but not in the SDSS. This is due to the difference in the surveys' major scientific goals and in the different sensitivities in the used cameras. The lack of the u filter in the PS1 is a disadvantage when it comes to finding RRL stars.

Additionally, the SDSS operates in a drift-scanning mode where the sky objects pass through its five different filters almost simultaneously. The correct colors of the observed sky objects can then be obtained unless they are variable on very short time scales (i.e., few minutes). Consequently, the SDSS drift-scanning technique gives the correct colors of RRL stars as these stars have periods in the ∼0.2–1 day range.

However, the correct colors of RRL stars are not provided with the PS1 photometric system because of the PS1 imaging technique. The PS1 images a selected patch of the sky with different filters at different times. Magnitudes in different filters correspond to different phases for short period variable objects like RRL stars.

Figure 1(a) illustrates the (g − r) versus (r − i) color–color diagram of stars in Stripe 82 from the seventh data release of the SDSS (SDSS DR7; Abazajian et al. 2009), and Figure 1(b) illustrates the PS1 (g_P1 − r_P1) versus (r_P1 − i_P1) color–color diagram for the same stars. Red dots represent a subsample of non-RRL stars while blue filled circles represent a subsample of the RRL stars detected in Stripe 82 (Sesar et al. 2010). While the RRL stars occupy a small and well-defined region in the SDSS color–color diagram (see Figure 1(a)), they are spread out over a large and wide region in the PS1 color–color diagram (see Figure 1(b)). This is a result of the different observing techniques used by the SDSS (near-simultaneous imaging using different filters) and PS1 (non-simultaneous imaging).

Zoom In Zoom Out Reset image size

We base our color cuts for selecting RRL candidates on colors from the SDSS DR7 photometric system and not on the colors from the PS1 photometric system due to the lack of the u band and of the true colors of RRL stars in the latter photometric system.

We use 636 pre-identified RRL stars selected from the catalogs of RRL stars in the CRTS (Drake et al. 2013) and LINEAR (Sesar et al. 2013) surveys for a better characterization of the SDSS colors and PS1 variability properties of RRL stars.

These 636 RRL stars are chosen based on their clean photometry in the SDSS DR7 and PS1 photometric systems. These stars have photometric errors of less than 0.2 in u and less than 0.1 in g, r, i, z, g_P1, and r_P1. These are primary objects that are not blended or saturated in both surveys and that have been observed more than twice by PS1 in both g_P1 (N$_{g_{P1} }\,{\ge}$ 3) and r_P1 (N$_{r_{P1} }\,{\ge}$ 3), respectively. N$_{g_{P1}}$ and N$_{r_{P1}}$ represent the number of PS1 observations in the g_P1 and r_P1 filters, respectively. The two last cuts were applied in order to study the variability of RRL stars in the PS1 multi-epoch data.

We corrected the magnitudes for extinction using the recalibration of Schlegel et al.'s (1998) dust map by Schlafly & Finkbeiner (2011). Since the RRL stars used here are located in areas where the extinction is small (i.e., at high Galactic latitudes), such color corrections can be used. The color densities of the 636 RRL stars in the SDSS photometric system are shown in Figure 2 where red and blue regions reflect large and small numbers of RRL stars, respectively. A sample of non-RRL stars are also plotted as small white dots to demonstrate the colors of these contaminant stars (i.e., main-sequence stars and stars in eclipsing systems). RRL stars occupy small areas in the color–color diagrams in Figure 2 and are concentrated in well-defined regions, especially in the (u − g) color, an advantage that helps in finding these stars.

Zoom In Zoom Out Reset image size

• 1.  
  We start by adopting initial rectangular color cuts from Sesar et al. (2010) to avoid downloading all the SDSS DR7 data in Stripe 82 (and later for the whole SDSS×PS1 footprint). Our RRL candidates must first pass the first four initial rectangular color cuts (Equations (6)–(9) in Sesar et al. 2010):

  $$\begin{equation} 0.75 < (u-g) < 1.45 \end{equation} \tag{ 1 }$$

  $$\begin{equation} {-}0.25 < (g-r) < 0.40 \end{equation} \tag{ 2 }$$

  $$\begin{equation} {-}0.20 < (r-i) < 0.20 \end{equation} \tag{ 3 }$$

  $$\begin{equation} {-}0.30 < (i-z) < 0.30. \end{equation} \tag{ 4 }$$

  These are single-epoch color ranges (Sesar et al. 2010) for RRab and RRc stars corrected for extinction using the Schlegel et al. (1998) dust map. The SDSS colors for RRL stars correspond to a random instant in their phase and depend on the time when the near-simultaneous SDSS photometry was obtained. It is thus safe to apply these color criteria to SDSS data but they are not suitable for PS1 data where the color range needs to be larger in order to account for the non-simultaneous observations in the PS1 filters. In order to avoid galaxies, these objects must be flagged as stars (${\tt type_{SDSS} = 6}$) in the SDSS. They should also be flagged as primary objects (${\tt mode_{SDSS} = 1}$) with clean photometry in the SDSS DR7 database (${\tt clean_{SDSS} = 1}$).Due to the noise and photometric errors resulting from the small number of PS1 epochs that we use in our method, some non-variable sources might appear as variables, especially faint sources with large photometric errors and bright sources that might saturate the CCD camera (see Section 4.3). To avoid this, we choose sources that are fainter than 14th and brighter than 20th magnitude in the PS1 g_P1 filter (Equation (5)). Although PS1 will eventually observe each object around 12 times in each filter, the survey is not finished yet and the average number of detections per star is ∼8 epochs in each of the PS1 filters. Some of these detections were not taken under good photometric conditions and therefore were flagged as bad sources by the PS1 pipeline. To ensure the reliability of our variability cuts, only clean PS1 detections that are not saturated or blended, and are not flagged as cosmic rays are used in our study (Morganson et al. 2012).Thus, we only choose stars that have been more than two clean detections in both the g_P1 and r_P1 filters (Equations (6) and (7)) in order to reliably distinguish variable from non-variable stars:

  $$\begin{equation} 14.0 < g_{P1} < 20.0 \end{equation} \tag{ 5 }$$

  $$\begin{equation} N_{g_{P1}} \ge 3 \end{equation} \tag{ 6 }$$

  $$\begin{equation} N_{r_{P1}} \ge 3. \end{equation} \tag{ 7 }$$

  Variability cuts in the i_P1, z_P1, and y_P1 filters are applied later.In the studied area of Stripe 82, we have ∼74,000 stars that passed the first four initial SDSS color cuts (Equations (1)–(4)), the PS1 magnitude cut (Equation (5)), and the PS1 threshold limit of the number of detections in both the g_P1 and r_P1 filters (Equations (6) and (7)).Although there are 374 RRL stars in the same area, we missed 85 of them. Around 92% of these 85 stars did not have more than two clean g_P1 or r_P1 PS1 detections (Equations (6) and (7)) while the rest 8% of the missed RRL stars did not pass all of the four SDSS color cuts (Equations (1)–(4)). This leaves us with 289 true RRL stars that we recovered in Stripe 82 (among the ∼74,000 stars that passed all the conditions in this step). The efficiency and completeness levels are then ∼0.39% (289/74, 000) and 77.3% (289/374), respectively.
• 2.  
  In order to optimize our color selection of RRL candidates, we define color selection boundaries using the 636 RRL stars (see Section 3.2) in the SDSS (u − g) versus (g − r) and (g − r) versus (r − i) color–color diagrams with GMM (VanderPlas et al. 2012). GMM is a Bayesian generative classification method that fits different classes with simple non-correlated Gaussians. These Gaussians are then used to compute the likelihood of a point to belong to each class. The class with the highest likelihood is the predicted result. In our case, GMM uses the colors of the 636 pre-identified RRL stars and compares them to the colors of non-RRL stars to find the GMM color selection boundaries. We choose this method instead of adopting sharp rectangular cuts (e.g., Vivas et al. 2001; Sesar et al. 2007, 2010) in order to optimize our efficiency and completeness levels when light curve analyses are not possible due to the small number of PS1 observations.In Figures 3 and 4, the GMM color selection boundaries are applied and plotted in green in the (u − g) versus (g − r) and (g − r) versus (r − i) color–color diagrams, respectively, for a subsample of stars in Stripe 82. The colors of the 636 pre-identified RRL stars used to find the GMM color selection boundaries are shown with black open circles. Stars that fall inside our GMM selection boundaries are shown as blue dots while stars that fall outside are plotted as red dots. Only stars that fall inside the GMM color selection boundaries in both color–color diagrams ((u − g) versus (g − r) and (g − r) versus (r − i)) are retained for further analysis.This step significantly reduces the number of stars in our sample from ∼74,000 to 1820 stars, out of which 260 are true RRL stars.Although the GMM color boundaries are computed using more than 600 well identified RRL stars distributed around the sky, 29 true RRL stars from Step (1) did not pass one or both of these GMM color boundary cuts. These stars either have relatively large SDSS magnitude uncertainties that are reflected in their colors or they fall close to, but outside of, our GMM color boundaries.Because 1820 stars passed all the cuts in this step (and the cuts in the previous step), and because we were able to recover 260 out of the 374 RRL stars found in Stripe 82, our efficiency level is ∼14.3% (260/1820) while the completeness level is 70% (260/374).
• 3.  
  After defining and applying the GMM selection boundaries for the SDSS colors in the previous section, we use the g_P1, r_P1, i_P1, z_P1, and y_P1 multi-epoch data from PS1 to distinguish a variable from a non-variable star.Since we cannot rely on our small number of PS1 detections to phase the light curves and find their periods, we calculate low-order statistics (e.g., standard deviation) and use them to define a GMM selection boundary cut for the g_P1 magnitudes as a function of the standard deviation in g_P1 ($\sigma _{g_{P1}}$) plus the standard deviation in r_P1 ($\sigma _{r_{P1}}$). In Figure 5, the GMM variability boundary plotted in green is computed by the GMM method (VanderPlas et al. 2012) which uses the ($\sigma _{g_{P1}}$ + $\sigma _{r_{P1}}$) values of the 636 pre-identified RRL stars compared to the ($\sigma _{g_{P1}}$ + $\sigma _{r_{P1}}$) values of non-variable stars to find the boundary of the variability cutoff. Although all of these 636 RRL stars are variable sources, only ∼90% of them fall above our variability boundary, while ∼10% show small or no variability due to the small number of epochs available from PS1. Only stars that fall above our GMM variability boundary are retained for further analysis. These stars have already passed the GMM selection boundaries for the SDSS colors discussed in the previous steps.In order to be considered as RRL candidates, stars that have more than two clean detections in the i_P1, z_P1, and y_P1 filters must pass the following additional variability criterion:

  $$\begin{equation} \sigma _{i_{P1}}+\sigma _{z_{P1}}+\sigma _{y_{P1}}\ge 0.1. \end{equation} \tag{ 8 }$$

  This threshold limit was adopted as more than 90% of the 636 pre-identified RRL stars (see Section 3.2) with more than two clean detections in the i_P1, z_P1, and y_P1 filters have $\sigma _{i_{P1}}+\sigma _{z_{P1}}+\sigma _{y_{P1}}\,{\ge}\, 0.1$. This criterion is applied to stars with $N_{i_{P1}}\,{\ge}\, 3$, $N_{z_{P1}}\,{\ge}\, 3$, and $N_{y_{P1}}\,{\ge}\, 3$ that have already passed all of our GMM color and variability selection boundaries. Stars that passed our GMM color and variability selection boundaries and that do not have more than two good detections in the i_P1, z_P1, and y_P1 filters are still considered RRL candidates.Applying the GMM variability selection cut (see Figure 5) for ($\sigma _{g_{P1}}$ + $\sigma _{r_{P1}}$) and the variability cut in the i_P1, z_P1, and y_P1 filters (see Equation (8)) to the 1820 RRL candidates from Step (2) reduces the number of RRL candidates to 255 stars, out of which 195 are true RRL stars and 60 are contaminant stars. We discuss the nature of the 60 contaminant stars in Section 4.3.At the same time, 65 RRL stars were lost when moving from Step (2) to Step (3). These stars did not show a significant amount of variability compared to other variable stars because their number of PS1 epochs is small (∼3) and their magnitudes in different detections are not significantly different as they have likely been multiply observed at a relatively close phase.In this final step, the efficiency significantly increases to ∼77% (195/255) and the completeness drops to ∼52% (195/374). This step greatly increases our efficiency level as it gets rid of a large fraction of non-variable stars with colors close to the colors of RRL stars (i.e., main-sequence stars with colors close to the colors of RRL stars).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We apply the regularly used color and variability rectangular cuts (Sesar et al. 2010) to the stars in Stripe 82 and compare the results using this technique with the results we achieved using the GMM technique to test whether the latter technique improves the recovery of RRL stars.

The first step in this technique is similar to Step (1) from the previous section where the number of RRL candidates is ∼74,000 stars of which 289 are known RRL stars. This step requires the SDSS rectangular color cuts, the magnitude cut, and the PS1 threshold limit of the number of detections in both, the g_P1 and r_P1 filters. The efficiency and completeness levels are then ∼0.39% (289/74, 000) and 77.3% (289/374), respectively.

Since we are not using the GMM technique in this section, we directly apply straight-line variability cuts in the PS1 filters. Stars with ($\sigma _{g_{P1}}+\sigma _{r_{P1}}\,{\ge}\,0.22$) that have passed the previous cuts in this section are retained for further analysis. The 636 pre-identified RRL stars were once again used to set the latter cut as more than 90% of these stars have $\sigma _{g_{P1}}+\sigma _{r_{P1}}\,{\ge}\,0.22$. Just like in Step (3), an additional cut ($\sigma _{i_{P1}}+\sigma _{z_{P1}}+\sigma _{y_{P1}}\,{\ge}\,0.1$) is applied for the retained stars with $N_{i_{P1}}\,{\ge}\,3$, $N_{z_{P1}}\,{\ge}\,3$, and $N_{y_{P1}}\,{\ge}\,3$. Retained stars that do not have more than two good detections in the i_P1, z_P1, and y_P1 filters are still considered RRL candidates.

There are ∼1600 stars that passed all of our cuts in this section of which 205 are known RRL stars. This yields efficiency and completeness levels of ∼13% (205/1600) and ∼54% (205/374), respectively.

The dependencies of the efficiency (dashed lines) and completeness (solid lines) levels in each step resulting from the GMM and rectangular cut techniques are plotted with red and blue lines in Figure 6, respectively. Although there was no significant change in the completeness level when using the rectangular cuts compared to the GMM technique, the efficiency level increased from ∼13% (205/1600, using rectangular cuts) to 77% ((195/255), using the GMM technique). Hence, we favor using the GMM technique in future studies.

Zoom In Zoom Out Reset image size

To understand the nature of the contaminant stars, we look for multi-epoch data in the CRTS database for the 60 contaminant stars we found in Stripe 82. Fifty six out of the 60 contaminant stars are found in the CRTS database and have been observed between ∼40 and 500 times.

Almost 40% of these stars showed no variability using the multi-epoch data from CRTS, which makes them non-variable stars that have passed our variability cuts. These stars were observed only three to four times with PS1 and have magnitudes close to our bright (g_P1 ∼ 14.0 mag) and faint (g_P1 ∼ 20.0 mag) magnitude cuts. Hence, it is not surprising that some non-variable sources passed our variability cuts as their variability statistics are based on a small number of observations where a single noisy epoch can bias the statistics and make a non-variable source appear as a variable one, and vice versa.

The remaining 60% of the contaminant stars in Stripe 82 appeared as non-RRL variable stars using the CRTS database. Their variability statistics reflected a change in their brightness over time but the shape of their light curve showed that most of them are W Ursae Majoris (W UMa), Algol binaries, δ Scuti, and SX Phe stars (Palaversa et al. 2013). Samples of the phased light curves for Algol binaries (P = 0.6684 day) and δ Scuti candidate (P = 0.11367021 day) stars that are contaminating our RRL stars are shown in panels (a) and (b) of Figure 7, respectively. We were able to recover the correct type and periods of these stars using the CRTS multi-epoch data. Using the current PS1 data available, there is no way of getting rid of all the contaminants.

Zoom In Zoom Out Reset image size

With the ∼77% (195/255) efficiency level computed in the previous section, we know that ∼23% of our RRL candidates are non-RRL stars (mainly non-variable stars and stars in eclipsing systems). However, we show in Section 6.1 that we are still able to detect halo substructures with such a contamination level. Additionally, the efficiency and completeness levels will be improved when more PS1 epochs are available in the near future. Our method can be useful in detecting RRL stars in surveys other than the PS1 where the number of detections per star is also small. Our efficiency and completeness levels as a function of g_P1 magnitudes are plotted in blue and red lines in Figure 8, respectively. The decrease in the efficiency and completeness levels as a function of magnitude reflects the increase in contamination for fainter stars.

Zoom In Zoom Out Reset image size

Using the 255 RRL candidates we detected in Stripe 82, we look for halo substructures in our covered distance range. We plot the number density distribution of these 255 RRL candidates in Figure 10. This plot includes our 60 contaminant stars in Stripe 82 (if we assume that Sesar et al.'s (2010) catalog of RRL stars is ≳ 99% complete).

Zoom In Zoom Out Reset image size

The density of the points that is accentuated by the white contours is shown in scaled density levels. The smoothed surface regions with a high number of stars are indicated in red while regions with low number of stars are indicated in dark blue. We recover the Hercules-Aquila cloud (Belokurov et al. 2007) at R.A.⁶ ∼−40° and d_h between ∼8 and ∼24 kpc. The trailing arm of the Sagittarius dwarf spheroidal's (dSph) tidal stream (Majewski et al. 2003; Law & Majewski 2010) is also recovered at R.A. ∼ 30° and d_h ∼ 23 kpc.

Both of our recovered substructures were seen using the ∼99% complete and efficient catalog of RRL stars in Stripe 82 (Sesar et al. 2010). Although our method is not as efficient and complete as the mentioned catalog, Figure 10 proves that the efficiency and completeness levels we achieved are good enough to select RRL candidates to trace stellar streams and substructures in spite of the inclusion of contaminant stars. Stripe 82 was visited ∼80 times by the SDSS, which made it relatively easy to find its RRL stars using light curve analysis (Sesar et al. 2010). In contrast, it was more difficult to find RRL stars in our study using only the SDSS colors and PS1 variability because of the small number of multi-epoch data available from PS1.

Nevertheless, we recovered ∼52% of the RRL stars (d_h < 70 kpc) not only in the Stripe 82 region, but in the whole SDSS×PS1 overlapping footprint. A detailed analysis of the distribution of the identified RRL candidates will be presented in a future paper.

Having additional PS1 epochs will improve the quality of our variability statistics which will improve the separation between variable and non-variable stars. Using the CRTS data, we showed in Section 4.3 that 40% of our contaminant stars are non-variable sources and have small number of PS1 epochs. We expect to get rid of at least 60% of these non-variable contaminant stars when more PS1 epochs (∼15 epochs in all filters) are available. However, it will not be possible to get rid of all of the contaminant stars as the number of PS1 epochs will not be sufficient to distinguish RRL from non-RRL variable stars using light curve and period analysis. Furthermore, having additional PS1 epochs will improve our completeness level as we missed many RRL stars due to the PS1 threshold limit of the number of detections in both the g_P1 and r_P1 filters (Equations (6) and (7)). We expect the efficiency and completeness level to increase to at least ∼83% and ∼65% when all of the PS1 epochs are available. Having additional epochs will also allow us to study stars in the halo that are further than 70 kpc (d_h).