FINDING, CHARACTERIZING, AND CLASSIFYING VARIABLE SOURCES IN MULTI-EPOCH SKY SURVEYS: QSOs AND RR LYRAE IN PS1 3π DATA

Pan-STARRS is a wide-field optical/near-IR survey telescope system located at Haleakala Observatory on the island of Maui in Hawaii. The PS1 survey (Kaiser et al. 2010) is collecting multi-epoch, multi-color observations undertaking a number of surveys, among which the PS1 3π survey (Chambers 2011) is the largest. It has observed the entire sky north of declination −30° in five filter bands (g_P1, r_P1, i_P1, z_P1, y_P1) with average wavelengths of 481, 617, 752, 866, and 962 nm, respectively (Stubbs et al. 2010; Tonry et al. 2012), with a 5σ single epoch depth of about 22.0, 22.0, 21.9, 21.0, and 19.8 mag in g_P1, r_P1, i_P1, z_P1, and y_P1, respectively. In contrast to the SDSS filters, the g_P1 filter extends 20 nm redward of g_SDSS, and the z_P1 filter reaches only to 920 nm. PS1 has no u band. In the near-IR, y_P1 covers the region from 920 nm to 1030 nm. More detailed descriptions of these filters and their calibration can be found in Stubbs et al. (2010) and Tonry et al. (2012).

Roughly 56% of the PS1 telescope observing time was dedicated to the PS1 3π survey, with an observing cadence optimized for the detection of near-Earth asteroids and slow-moving solar system bodies. The PS1 3π survey plan is to observe each position four times per filter per year, where the epochs are typically split into two pairs of exposures per year per band (transit-time-interval (TTI) pairs) taken ∼25 minutes apart in the same band. This allows for the discovery of moving or variable sources. The sky north of declination −30° was planned to be observed four times in each band pass per year (Chambers 2011). Through periods of bad weather and telescope downtime in practice, fewer epochs were observed.

Images are automatically processed using the survey pipeline (Magnier et al. 2008), performing bias subtraction, flat fielding, astrometry, photometry, as well as image stacking and differencing. The photometric calibration of the survey is better than one hundredth of a magnitude (Schlafly et al. 2011).

All data processing shown here is carried out under PS1 catalog processing version PV2, where the average number of total detections per source is 55 over 3.7 years, including some data taken in non-photometric conditions.

We perform a number of cuts on the PS1 data to remove outliers and unreliable data. These cuts fall into two categories: detection cuts that remove individual detections, and object cuts that remove all detections of a source from the analysis.

2.2.1. Detection Cuts

The most important detection cut we apply is to remove data taken in non-photometric conditions, according to Schlafly et al. (2011), and data from any Orthogonal Transfer Array (OTA) where the detections of bright stars on that chip are on average over 0.02 mag too faint. These cuts remove about 30% of detections.

The second most important detection cut we apply is to remove observations which land on bad parts of the detector, as indicated by having psf_qf_perfect < 0.95. This removes about 10% of detections. Similarly important, we exclude any observation where the point-spread function (PSF) magnitude is inconsistent with the aperture magnitude by more than 0.1 mag or four times the estimated uncertainties, removing 10% of detections.

We remove any detections with problematic conditions noted by the PS1 pipeline, according to the detections' flags. For the cleaning flags used, see Table 1 and also Magnier et al. (2013). This eliminates only about 2% of detections.

Table 1.  Bit-flags Used to Exclude Bad or Low-quality Detections

FLAG NAME	Hex Value	Description
PM_SOURCE_MODE_FAIL	0x00000008	Fit (nonlinear) failed (non-converge, off-edge, run to zero)
PM_SOURCE_MODE_POOR	0x00000010	Fit succeeds, but low-SN or high-Chisq
PM_SOURCE_MODE_SATSTAR	0x00000080	Source model peak is above saturation
PM_SOURCE_MODE_BLEND	0x00000100	Source is a blend with other sources
PM_SOURCE_MODE_BADPSF	0x00000400	Failed to get good estimate of object's PSF
PM_SOURCE_MODE_DEFECT	0x00000800	Source is thought to be a defect
PM_SOURCE_MODE_SATURATED	0x00001000	Source is thought to be saturated pixels (bleed trail)
PM_SOURCE_MODE_CR_LIMIT	0x00002000	Source has crNsigma above limit
PM_SOURCE_MODE_MOMENTS_FAILURE	0x00008000	Could not measure the moments
PM_SOURCE_MODE_SKY_FAILURE	0x00010000	Could not measure the local sky
PM_SOURCE_MODE_SKYVAR_FAILURE	0x00020000	Could not measure the local sky variance
PM_SOURCE_MODE_BIG_RADIUS	0x00100000	Poor moments for small radius, try large radius
PM_SOURCE_MODE_SIZE_SKIPPED	0x10000000	Size could not be determined
PM_SOURCE_MODE_ON_SPIKE	0x20000000	Peak lands on diffraction spike
PM_SOURCE_MODE_ON_GHOST	0x40000000	Peak lands on ghost or glint
PM_SOURCE_MODE_OFF_CHIP	0x80000000	Peak lands off edge of chip

Download table as:  ASCIITypeset image

Finally, we apply an outlier cleaning based on the z-score of the individual measurements ${z}_{i}=({m}_{i}-\mu ({b}_{i}))/{\sigma }_{i}$, where m_i is a given magnitude measurement, σ_i is its uncertainty, and μ(b_i) is the error-weighted mean magnitude of all measurements of that source in its band b_i. This eliminates 2% of detections, and we limit it to eliminate at most 10% of the detections of any individual source.

Figure 1 gives the number of epochs, as well as their cadence, in each band after all of these cuts have been applied. The average number of surviving epochs per source is 35 rather than the total 55 observations.

Zoom In Zoom Out Reset image size

The detection cuts we make are summarized in Table 2. We note that if a detection has one problematic condition, it is more likely than not to also be affected by other problematic conditions.

Table 2.  Cuts Used to Exclude Bad Detections

Condition	Fraction of Detections Removed
Photometric conditions	0.29
$| {\mathtt{ap}}\_{\mathtt{mag}}\quad -\quad {\mathtt{psf}}\_{\mathtt{inst}}\_{\mathtt{mag}}| \;\lt$ max(4 × σm, 0.1)	0.10
psf_qf_perfect > 0.95	0.11
Pipeline flags (Table 1)	0.017
$| {z}_{i}-{z}_{\mathrm{median}}| \lt 5\sigma$	0.02

Download table as:  ASCIITypeset image

2.2.2. Object Cuts

WISE is a NASA infrared-wavelength astronomical space telescope providing mid-infrared data with far greater sensitivity than any previous survey. It performed an all-sky survey with imaging in four photometric bands over 10 months (Wright et al. 2010). Nikutta et al. (2014) have shown that the color $W12=W1-W2\gt 0.5$ is an excellent criterion to isolate QSOs, because W12 is an indicator of the hot dust torus in AGNs. To aid in the QSO identification, we want to find objects with these unusual W12 colors, but need to make sure that these colors are not merely a consequence of poor WISE photometry. For objects with good measurements (σ_W1 < 0.3, σ_W2 < 0.3), we use W12 as a parameter for classification. 1.46 × 10⁸ of the 3.88 × 10⁸ selected objects from Section 2.1 have reliable W12 (σ_W1 < 0.3, σ_W2 < 0.3, where σ_W1, σ_W2 are the errors given on the WISE magnitudes).

The SDSS (York et al. 2000) is a major multi-filter imaging and spectroscopic survey using a dedicated 2.5-m wide-angle optical telescope at Apache Point Observatory in New Mexico, United States. The Sloan Legacy Survey covers about 7500° of the Northern Galactic Cap in optical ugriz filters, with average wavelengths of 355.1, 468.6, 616.5, 748.1, and 893.1 nm. In typical seeing, it has a 95% completeness down to magnitudes of 22.0, 22.2, 22.2, 21.3, and 20.5, for u, g, r, i, and z, respectively. Additionally, the Sloan Legacy Survey contains three stripes in the South Galactic Cap totaling 740 square degrees. The central stripe in the South Galactic Cap, Stripe 82 (S82), was scanned multiple times to enable a deep co-addition of the data and to enable discovery of variable objects. S82 has ∼60 epochs of imaging data in ugriz, taken over ∼5 years, where extensive spectroscopy provides a reference sample of nearly 10,000 spectroscopically confirmed quasars (Schneider et al. 2007; Schmidt et al. 2010). For S82, there is also a sample of 483 identified RR Lyrae available (Sesar et al. 2010).

The classification of QSOs and RR Lyrae in SDSS S82 will be used as a ground truth. This means, they will be used as training set for classification as well as for testing how well our classification method works (see Section 3).

Within −50° < α < 60°, −125 < δ < 125, there are 9073 QSO and 482 RR Lyrae from the samples mentioned above. Out of these, 7633 QSOs and 415 RR Lyrae are cross-matched to objects in our PS1 selection. We also select more than 1.85 × 10⁶ "other" objects from S82. Of these, 10% are missing because of the cuts of Section 2.2, and the remaining objects are outside our magnitude range of interest.

We start by laying out a very generic and non-parameteric measure for variability, simply to characterize the significance of variability by a scalar quantity. Specifically, we define

$$\begin{eqnarray}&&{\hat{\chi }}^{2}=\displaystyle \frac{{\chi }_{\mathrm{source}}^{2}-{N}_{\mathrm{dof}}}{\sqrt{2\;{N}_{\mathrm{dof}}}},\end{eqnarray} \tag{ 1 }$$

with

$$\begin{eqnarray}&&{\chi }_{\mathrm{source}}^{2}=\displaystyle \sum _{\lambda }\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{({m}_{\lambda ,i}-\langle {m}_{\lambda }\rangle )}^{2}}{{\sigma }_{\lambda ,i}^{2}},\end{eqnarray} \tag{ 2 }$$

where N is the total number of photometric points for one object across all n bands, the sum over λ is over the PS1 bands g_P1, r_P1, i_P1, z_P1, y_P1, and ${N}_{\mathrm{dof}}=N-n$ is the number of degrees of freedom.

Assuming that most of the sources are not variable, we expect the distribution of ${\hat{\chi }}^{2}$ to be a unit Gaussian distribution. In contrast, varying sources should form a "tail" of higher ${\hat{\chi }}^{2}$. Figure 3 shows the normalized distribution of ${\hat{\chi }}^{2}$, derived from the PS1 photometry of all selected objects in S82, with known QSOs (blue) and known RR Lyrae (red) shown in separate (normalized) distributions. The "other" objects have a ${\hat{\chi }}^{2}$-distribution close to that expected for non-varying sources (dashed line), confirming that most sources in the sky are non-varying (within a level of less than a few percent) and that the PS1 photometry is reliable. The QSOs and RR Lyrae appear well separated in the normalized distibutions. However, there are only 415 RR Lyrae and 7630 QSOs, compared to ∼1.85 × 10⁶ "other" objects in SDSS S82 cross-matched to PS1 and surviving the cuts of Section 2.2. Figure 3(b) shows how the distribution of "other" sources superimposes the distribution of QSOs and RR Lyrae due to the high number of "other" sources.

Zoom In Zoom Out Reset image size

Therefore, a simple criterion such as ${\hat{\chi }}^{2}$ is insufficient to identify QSOs or RR Lyrae. In the subsequent analysis, all objects are used, though for RR Lyrae only, one could in principle restrict oneself to objects with ${\hat{\chi }}^{2}\gt 10$ without losing completeness.

Beyond simply establishing variability, variable sources can and should be characterized by the amplitude of their variability and the timescales over which they vary. A useful and well-established tool for this is the structure function (Hughes et al. 1992; Collier & Peterson 2001; Kozłowski et al. 2009): it gives the mean squared magnitude difference (in a given band) between pairs of observations of some object's brightness (Δm) as a function of the time lag between the observations (Δt). There are various ways to parameterize such a structure function, and the damped random walk (DRW) is a useful function family. For a DRW, the structure function is specified by two parameters, τ and ${V}_{\infty },$ and is given by

$$\begin{eqnarray}&&V({\rm{\Delta }}t| \tau ,{V}_{\infty })={V}_{\infty }(1-{e}^{-| {\rm{\Delta }}t| /\tau }).\end{eqnarray} \tag{ 3 }$$

In this notation V(Δt) reflects the expectation value for the squared magnitude difference, Δm², among measurements separated in time by Δt; V(Δt) is simply $\sqrt{2}$ times the expected magnitude variance during Δt. ${V}_{\infty }$ is conventionally denoted as ω², and τ is called the decorrelation time of the DRW. The source variability is then characterized by two structure function parameters, ω and τ.

Objects of different classes typically occupy different regions in structure function parameter space. As we show below, the likelihood $p({\boldsymbol{m}}\;| \;\mathrm{SF}\ \ \mathrm{parameters})$ can be used to select remarkably pure and complete samples of QSOs as well as RR Lyrae, which makes selection by structure function parameters an efficient approach for both selecting stochastically varying and periodic variable objects (Schmidt et al. 2010).

The cadence of surveys like the SDSS provides data that allows application of the usual single-band formulation of structure functions. However, the cadence of PS1 3π data, which observes in different bands at different epochs (see Section 2), makes it necessary to extend this approach for multi-band fitting. We need a practical approach to turn the ∼6–9 epoch PS1 light curves in each band into a ∼35 epoch overall light curve. If objects were to vary the same way in all observed bands, implementing such an approach would simply amount to determining the (time-averaged) mean color of the object and shifting the light curves in the different bands to a common magnitude. However, in practice, most astrophysical objects vary more at shorter wavelengths. To account for this, the multi-band model we present here has, beyond ω and τ, a set of temporal mean magnitude parameters in each PS1 band, ${\boldsymbol{\mu }}$, and it links the variability amplitudes ω(b) in different bands b by a power law with exponent α. Specifically,

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{\mathrm{log}(\omega (b)/\omega (r))}{\mathrm{log}({\lambda }_{b}/{\lambda }_{r})},\end{eqnarray} \tag{ 4 }$$

where λ_b is the effective wavelength of the band b.

To assign a likelihood to an object's photometry, given a structure function model, we make use of a Gaussian process formulation for stochastic source variability. In contrast to single-band structure function models (e.g., Rybicki & Press 1992; Zu et al. 2011; Hernitschek et al. 2015), the Gaussian process is not applied to any particular band but instead to an arbitrarily constructed fiducial band which can be scaled and shifted onto the particular bands. This permits simultaneous treatment of multiple bands, even when the bands are not observed simultaneously at the same epochs. It is key in this context to realize that the fiducial band is a latent variable—it is never directly observed; only the scaled and shifted versions are observed, where substantial measurement noise is present.

The fiducial light curve can be described with a zero-mean and unit characteristic variance Gaussian process. That is, the prior probability distribution function (pdf) for a set of N fiducial "magnitudes" ${\boldsymbol{q}}$ that are instantiated at observed times t_n is a multivariate normal distribution:

$$\begin{eqnarray}&&p({\boldsymbol{q}})={ \mathcal N }({\boldsymbol{q}}\;| \;0,{C}^{q}),\end{eqnarray} \tag{ 5 }$$

where C^q is a N × N symmetric positive definite covariance matrix. In the case of a DRW model, C^q is given by

$$\begin{eqnarray}&&{C}_{{{nn}}^{\prime }}^{q}=\mathrm{exp}[-\displaystyle \frac{| {t}_{n}-{t}_{{n}^{\prime }}| }{\tau }].\end{eqnarray} \tag{ 6 }$$

This is identical to the usual single band DRW covariance matrix, except that we have dropped a scale factor ω² from Equation (6), because we have defined the fiducial band q to have unit variance. This factor reappears in our multi-band structure function through the scale factors that link the fiducial band to observed bands.

We consider now a given source for which we have N observations across N_band different bands. The data consist of the magnitude and uncertainty vectors ${\boldsymbol{m}}$ and ${\boldsymbol{\sigma }}$, the times of observation t_n, and the corresponding bands b_n. The source also has N_band temporal mean magnitudes ${\boldsymbol{\mu }}$. We define the N × N_band matrix ${\mathbb{M}}$ so that

$$\begin{eqnarray}&&{\mathbb{M}}{\boldsymbol{\mu }}=[\mu ({b}_{1}),\mu ({b}_{2}),\cdots ,\mu ({b}_{N})].\end{eqnarray} \tag{ 7 }$$

The likelihood of an individual measurement m_n, given its observational uncertainty σ_n and a value for the corresponding fiducial magnitude q_n, is found by shifting and scaling the fiducial magnitude and adding Gaussian noise. This makes the single-datum likelihood

$$\begin{eqnarray}&&p({m}_{n}\;| \;{q}_{n},{b}_{n},{\sigma }_{n}^{2})={ \mathcal N }({m}_{n}\;| \;\omega ({b}_{n}){q}_{n}+\mu ({b}_{n}),{\sigma }_{n}^{2}),\end{eqnarray} \tag{ 8 }$$

where ω(b_n) is the variability in bandpass b_n relative to the unit variability of the unobserved fiducial band.

Introducing the diagonal N × N matrix Ω, defined by Ω_ii = ω(b_i), the full likelihood is given by

$$\begin{eqnarray}&&p({\boldsymbol{m}}\;| \;{\boldsymbol{q}},{\rm{\Sigma }})={ \mathcal N }({\boldsymbol{m}}\;| \;{\rm{\Omega }}{\boldsymbol{q}}+{\mathbb{M}}{\boldsymbol{\mu }},{{\rm{\Sigma }}}^{2}),\end{eqnarray} \tag{ 9 }$$

where Σ is a diagonal matrix with Σ_ii = σ_i. Because everything is Gaussian, the latent fiducial magnitudes never have to be explicitly inferred; they can all be marginalized out analytically. This marginalization leads to the likelihood given the model, and the covariance matrix of the data:

$$\begin{eqnarray}&&p({\boldsymbol{m}}\;| \;\mathrm{SF}\ \ \mathrm{parameters},{\boldsymbol{\mu }})={ \mathcal N }({\boldsymbol{m}}\;| \;{\mathbb{M}}{\boldsymbol{\mu }},C)\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&C={\rm{\Omega }}{C}^{q}{\rm{\Omega }}+{{\rm{\Sigma }}}^{2}.\end{eqnarray} \tag{ 11 }$$

This is identical to the case of a single-band DRW model, except the rows and columns of C^q are scaled by amplitudes ω(b_n), $\omega ({b}_{{n}^{\prime }})$ for the bands b_n and ${b}_{{n}^{\prime }}$, and a contribution from the photometric uncertainties is added to the diagonal:

$$\begin{eqnarray}&&{C}_{{{nn}}^{\prime }}=\omega ({b}_{n})\omega ({b}_{{n}^{\prime }})\mathrm{exp}[-\displaystyle \frac{| {t}_{n}-{t}_{{n}^{\prime }}| }{\tau }]+{\sigma }_{n}^{2}\;{\delta }_{{{nn}}^{\prime }}.\end{eqnarray} \tag{ 12 }$$

Equations (10) through (12) provide a method for computing the probability of any set of observed magnitudes, given their metadata and the parameters ω(b), τ, and ${\boldsymbol{\mu }}$.

We are primarily interested in the structure function parameters and are relatively uninterested in the mean magnitudes ${\boldsymbol{\mu }}$. This is exactly the same situation as in Zu et al. (2011). Following that work, the likelihood of the structure function parameters, given the data, marginalized over ${\boldsymbol{\mu }}$, is given by:

$$\begin{eqnarray}&&p({\boldsymbol{m}}\;| \;\mathrm{SF}\ \ \mathrm{parameters})={ \mathcal L }\propto | C{| }^{-1/2}| {C}_{\mu }{| }^{1/2}\mathrm{exp}(-{\chi }^{2}/2)\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}{C}_{\mu } & = & {({{\mathbb{M}}}^{{\rm{T}}}{C}^{-1}{\mathbb{M}})}^{-1}\;\\ {\chi }^{2} & = & {({\boldsymbol{m}}-{\mathbb{M}}{\boldsymbol{\mu }})}^{{\mathsf{T}}}{C}^{-1}({\boldsymbol{m}}-{\mathbb{M}}{\boldsymbol{\mu }}).\end{eqnarray} \tag{ 14 }$$

We note that the factor of $| {C}_{\mu }{| }^{1/2}$ in Equation (13) comes from the marginalization over ${\boldsymbol{\mu }}$. We maximize $p({\boldsymbol{m}}\;| \;\mathrm{SF}\ \ \mathrm{parameters})$ to obtain best fit values of the structure function parameters. We then obtain ${\boldsymbol{\mu }}$ as the maximum likelihood values of ${\boldsymbol{\mu }}$ given the structure function parameters. That is, the mean magnitudes are given by

$$\begin{eqnarray*}&&{\boldsymbol{\mu }}={({{\mathbb{M}}}^{T}{C}^{-1}{\mathbb{M}})}^{-1}{{\mathbb{M}}}^{{\rm{T}}}{C}^{-1}{\boldsymbol{m}},\end{eqnarray*}$$

and have variance ${C}_{\mu }.$

One advantage of this approach is that it can be used to predict unobserved databased on observed data. Because both the process is Gaussian and the noise is assumed to be Gaussian, conditional predictions of the magnitudes can be made given the observed data and the structure function. The analysis is exactly the same as in Rybicki & Press (1992), with the exception that we adopt the multi-band structure function C of Equation (12). The magnitudes ${\tilde{m}}_{k}$ at K unmeasured times t_k, taken through bandpasses b_k, conditioned on the data in hand, are given by:

$$\begin{eqnarray}&&p(\tilde{m}| {\boldsymbol{m}})={ \mathcal N }(\tilde{m}| \tilde{\mu },\tilde{C})\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\tilde{\mu }={\boldsymbol{\nu }}+X\cdot {C}^{-1}\cdot [{\boldsymbol{m}}-{\mathbb{M}}{\boldsymbol{\mu }}]\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\tilde{C}=Y-X\cdot {C}^{-1}\cdot {X}^{{\mathsf{T}}}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\boldsymbol{\nu }}=[\mu ({b}_{1}),\mu ({b}_{2}),\cdots ,\mu ({b}_{K})].\end{eqnarray} \tag{ 18 }$$

In the case of a multi-band DRW model,

$$\begin{eqnarray}&&{X}_{{kn}}=\omega ({b}_{k})\omega ({b}_{n})\mathrm{exp}[-\displaystyle \frac{| {t}_{k}-{t}_{n}| }{\tau }]\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{Y}_{{{kk}}^{\prime }}=\omega ({b}_{k})\omega ({b}_{{k}^{\prime }})\mathrm{exp}[-\displaystyle \frac{| {t}_{k}-{t}_{{k}^{\prime }}| }{\tau }].\end{eqnarray} \tag{ 20 }$$

Here $\tilde{m}$ is the column vector of conditional predictions, $\tilde{\mu }$ and $\tilde{C}$ are a conditional mean vector and a conditional variance matrix, (temporary) mean vector ν is K-dimensional, and the matrices $\tilde{C}$, X, and Y are N × N, K × N, and K × K respectively. Vectors ${\boldsymbol{m}}$ and ${\boldsymbol{\mu }}$ and matrix C are defined in Section 3.2.

We have described a general technique for determining structure functions for multi-band, non-simultaneous data. The key ingredient is a description of the ratios of the variabilities in the different bands, which we characterize by a power law with exponent α (Equation (4)). The other elements of the structure function analysis—overall variability, timescale, and linear nuisance parameters (mean magnitudes)—are the same as in the single band case. For the case of the PS1 data at hand, it turns out that α is poorly constrained for individual objects, making it preferable to derive an external estimate of α from other data (SDSS S82), and fix it for the subsequent PS1 analysis. Assuming Equation (4), we used data from SDSS S82 to derive characteristic values for α. We found α ≈ −0.65 for QSOs and α ≈ −1.3 for RR Lyrae, both with an uncertainty of 0.01, which is in good agreement with Sesar (2012). We experimented fitting PS1 data with both choices of α, and obtained similarly good fits. Accordingly, we decided to adopt a single fixed α = −0.65 throughout this analysis. This choice of α corresponds to variability ratios ω(b)/ω(r) = 1.175, 1.00, 0.88, 0.80, 0.75, where b represents the PS1 bands g_P1, r_P1, i_P1, z_P1, and y_P1.

With α fixed, our fit to each source is described by the timescale τ, an overall variability scaled to the r band, ω_r, and the mean magnitudes ${\boldsymbol{\mu }}$. We calculate these from the PS1 photometry of all sources within the SDSS S82 area.

Figure 4 shows for example fits to the PS1 photometry of objects in SDSS S82: one QSO, one RR Lyrae, one "other" variable object and one seemingly non-varying object. For each object we show the light curve as observed in the five bands (top panel), and the combined light curve after shifting each band by the estimated μ(b); the structure function parameters ω_r and τ are listed for each case. Note that the QSO in Figure 4 (a) has τ of over a year, while the RR Lyrae in panel (b) has a τ of about a day. The figure also shows the interpolated light curves, given the observations and the structure function parameters, according to the technique of Rybicki & Press (1992).

Zoom In Zoom Out Reset image size

One could sensibly derive the pdf's for the parameters ω_r, τ, and ${\boldsymbol{\mu }}$ via MCMC; however, it proved computationally easier to calculate $p(m| {\omega }_{r},\tau )$ based on a reasonable parameter grid, and to do the linear optimization of the ${\boldsymbol{\mu }}$ for each grid-point. We use a log-spaced grid of $-2\lt \mathrm{log}{\omega }_{r}\lt 0.5$, 0.04 < τ < 5000 with 20 values evenly spaced in log ω_r and 30 in log τ to find the best-fit structure function parameters on the grid ω_r,grid and τ_grid. We have verified that this approximation agrees well with full MCMC runs. Figure 5 shows the gridded log-likelihood estimates for the same four sources as in Figure 4. The panels show the 68% CI of the $\mathrm{log}{ \mathcal L }$ distribution and the maximum likelihood values of the parameters.

Zoom In Zoom Out Reset image size

Figure 6 shows the distribution of variability parameters ω_r and τ, for all PS1 objects in the SDSS S82 area that survive the magnitude cut and which have significant variability, either satisfying ${\hat{\chi }}^{2}\gt 5$ or ${\hat{\chi }}^{2}\gt 30$ for objects within the stellar locus. This figure illustrates a number of points: first, and unrelated to variability, it shows the power of the WISE color $W1-W2$ to separate QSOs from other sources (Nikutta et al. 2014). Second, it shows that RR Lyrae and QSOs indeed populate different areas of (ω_r, τ) space. While they can only be roughly differentiated by their amplitudes ω_r, they have dramatically different timescales τ: RR Lyrae have typical τ ∼ 1 day and QSOs have τ ∼ 100–1000 days.

Zoom In Zoom Out Reset image size

We note that we also explored a power-law model for the structure function and found that it provided worse separation between QSO, RR Lyrae, and other objects in the structure function parameter plane. This can be explained by the cadence of the survey, as the definition of the power law makes the structure-function fitting more sensitive to the TTI pairs.

Figures 3 and 6 show that the light curve parameters will be very helpful in classifying variable sources. Yet, these figures also show that simple cuts on some parameters will not be optimal for differentiating object classes. A more sophisticated machine-learning method is needed here.

For classifying objects based on variability measures and mean magnitudes calculated before, we use the RFC (implemented in Python's scikit_learn package). Using a training set, it will give the classification probability of a target set's object being of a certain class, p_QSO and p_RRLyrae. However, we treat them as arbitrary numbers and not as probabilities, and instead calculate purity and completeness of the sample later on.

For using a RFC, a training set is needed, with observed object parameter values as well as classification labels. We use the classification of QSOs (Schneider et al. 2007; Schmidt et al. 2010) and RR Lyrae (Sesar et al. 2010) in SDSS S82 as a ground truth. Positions from SDSS S82 are cross-matched to PS1 positions within 1 arcsec. The object selection and outlier cleaning of Section 2.2 are applied. This results in a set of 7633 QSO, 415 RR Lyrae and more than 1.87 × 10⁶ other objects. This "training set template" is then extended to deal both with magnitude uncertainties and different amounts of reddening.

As an RFC cannot deal with measurement uncertainties by default, we deal with measurement uncertainties by extending the "training set template" by copies of itself, sampled within the assumed errors of the PS1 and WISE data. We take five samples for each object in the training set, in addition to the original one.

We additionally extend the training set to account for uncertainties in reddening. Dereddening is done using the reddening-based $E(B-V)$ dust map from Schlafly et al. (2014). We extend the training set by presenting additional QSOs, RR Lyrae, and other objects to the classifier, where we have artificially introduced a small dereddening error. We do this in the following way:

• (i)  
  make $E{(B-V)}_{\mathrm{sample}}$ drawn from Gaussian $G{(E{(B-V)}_{\mathrm{catalog}},\delta E(B-V)=0.1E(B-V)}_{\mathrm{catalog}})$ at the position of the training set source;
• (ii)  
  5% chance that $E{(B-V)}_{\mathrm{sample}}=0$, irrespective of catalog entry;
• (iii)  
  sample new mean magnitudes in bands g_P1, r_P1, i_P1, z_P1, y_P1 for PS1, and W1, W2 from WISE within their errors;
• (iv)  
  deredden them by $E{(B-V)}_{\mathrm{sample}}$;
• (v)  
  brighten magnitudes so that r_PS1 after dereddening by $E{(B-V)}_{\mathrm{sample}}$ is the same as after dereddening by $E{(B-V)}_{\mathrm{catalog}}$.

Each of the sources within the "training set template" is re-sampled five times to make the training set. This results in a training set having 38165 sources labeled as "QSO," and 2075 sources labeled as "RR Lyrae" at different reddenings, and a few million the "other" objects.

In principle, this technique could be used to train a classifier that could robustly classify sources even at large $E(B-V)$. In practice, however, our current classifier does not operate reliably at large $E(B-V)$. This is because in our training set, we use only objects on Stripe 82, where the reddening is small, and so no objects with large reddenings, and, correspondingly, large reddening errors, exist in the training set. We defer to later work accurate treatment of highly reddened objects in the plane, and focus on the lightly reddened high-latitude regions in this paper.

When using an RFC, missing values have to be replaced by some dummy values ("imputation") in the training and target sets. A common solution is replacing missing values by the mean of the available ones. This can be done not only for missing values, but also for unreliable values. As imputation of the median is impractical for the way we process the data, we had tested if an imputation of −9999.99 instead behaves comparably. We found that using −9999.99 versus median had no effects on our results.

In addition to imputing missing values, we use imputation when values are considered as unreliable. Accordingly, we impute a value of −9999.99 also in cases where σ_W1 > 0.3, σ_W2 > 0.3, or when magnitude errors are not available.

The Table 3 summarizes the parameter set being used for the RFC.

Table 3.  Parameter Set for the Random Forest Classifier

Parameter	Description
ωr,grid,τgrid	best fit structure function parameter on log-spaced grid
${\hat{\chi }}^{2}$	normalized χ2 statistic, see Equation (1)
${(g-r)}_{{\rm{P1}}}$, ${(r-i)}_{{\rm{P1}}}$, ${(i-z)}_{{\rm{P1}}}$, ${(z-y)}_{{\rm{P1}}}$	colors from dereddened PS1 mean magnitudes
$\langle r{\rangle }_{{\rm{P}}1,\mathrm{deredd}}$	dereddened PS1 mean rP1 magnitude, only used for calculation of pQSO
W12	$W1-W2$, helps with QSO identification
${i}_{{\rm{P1}}}-W1$	separates RR Lyrae from QSO

Download table as:  ASCIITypeset image

Though the mean r band magnitude is helpful in detecting RR Lyrae in general, we do not use it here as it introduces too strong a bias in distance, as the training set covers only the range 14.5 ≲ r_P1 ≲ 21.5 and we want to identify candidates fainter than 20.25 mag. Among the colors, the dereddened ${(i-z)}_{{\rm{P1}}}$ is a helpful gravity indicator that helps to reduce contamination (Vickers et al. 2012).

Table 4.  The Catalog of Variable Sources in PS1 3π

Column	FITS Format Code	Description
1	E	R.A. in degrees
2	E	decl. in degrees
3	E	scalar variability quantity ${\hat{\chi }}^{2}$, Equation (1)
4	E	best fit structure function parameter ωr (r band variability amplitude) on log-spaced grid
5	E	best fit structure function parameter τ (timescale) on log-spaced grid
6	E	error-weighted mean gP1 band magnitude $\langle {g}_{{\rm{P1}}}\rangle$
7	E	error-weighted mean rP1 band magnitude $\langle {r}_{{\rm{P1}}}\rangle$
8	E	error-weighted mean iP1 band magnitude $\langle {i}_{{\rm{P1}}}\rangle$
9	E	error-weighted mean zP1 band magnitude $\langle {z}_{{\rm{P1}}}\rangle$
10	E	error-weighted mean yP1 band magnitude $\langle {y}_{{\rm{P1}}}\rangle$
11	E	W1 – W2 color from WISE
12	E	pQSO
13	E	pRRLyrae

Note. Structure of the Catalog of Variable sources in PS1 3π. The Catalog of Variable sources is available in its entirety in machine-readable format. A FITS file is available.

(This table is available in its entirety in FITS format.)

Download table as:  ASCIITypeset image

In order to test the efficacy of the selection and classification method, we carried out detailed testing on the S82 area, using PS1 light curves, with the object classifications from S82 (Schneider et al. 2007; Sesar et al. 2010) as the "ground truth." To quantify purity and completeness of our classifications, we use S82 both as the training and validation set. A randomly selected 50% of the ≳1.85 × 10⁶ cleaned S82 objects is used for creating the training set, with the other half as the validation set.

For any one of the two categories, say RR Lyrae, we can define a candidate sample ${ \mathcal S }$ by the choice of a minimum p_RRLyrae. We can then calculate on the basis of the S82 ground truth the completeness and the purity of this sample. Here, purity is defined as the fraction of all RR Lyrae stars in ${ \mathcal S }$, and the "completeness" is the fraction of actual RR Lyrae stars contained in ${ \mathcal S }$. In both instances, we would expect completeness to be monotonic and purity to be nearly monotonic in p_RRLyrae.

For the QSOs, analogous definitions apply. Depending on context, we describe a sample ${ \mathcal S }$ either by a cut on p_RRLyrae/QSO, or by the corresponding purity and completeness of this sample as determined on Stripe 82.

Figure 7 shows precision-recall curves (Powers 2011) for the trade-off between purity and completeness with respect to the total cross-matched sources. These values are calculated for sources fulfilling the criteria of having all $\langle {g}_{{\rm{P}}1}\rangle$, $\langle {r}_{{\rm{P}}1}\rangle$, $\langle {i}_{{\rm{P}}1}\rangle$ available and between 15 and 20. We also show the case of all sources ($15\lt \langle {g}_{{\rm{P}}1}\rangle$, $\langle {r}_{{\rm{P}}1}\rangle$, $\langle {i}_{{\rm{P}}1}\rangle \lt 21.5$) using all available information as dashed blue lines.

Zoom In Zoom Out Reset image size

The left column refers to QSO classification, the right one to RR Lyrae classification. This figure shows that, as expected, for small completeness the purity is maximal, while the completeness is maximized with severe expense to the purity. What compromise needs to be made between completeness and purity in sample selection depends in detail on the science question, but the top panels of figure 7 suggest that the purity increases only a little at the expense of completeness, less than 80%. This may be a sensible threshold for an inclusive sample, whenever PS1 light curves and mean colors, as well as WISE colors are available. At the top of the horizontal axis we have indicated the relation between completeness and p_RRLyrae, p_QSO. Using probability thresholds of ∼0.05, we get purity both for QSO and RR Lyrae at a level of 70%, and completeness at a level of 98%. Using probability thresholds of 0.2, we get purity at a level of 76%, completeness of 94% for both QSO and RR Lyrae.

The different lines in the upper panels of Figure 7 illustrate the relative importance of the different pieces of data that may enter the classification; we have not only carried out classification with the full parameter set from Table 3, but also tested the cases where only color-related or variability-related information were used.

For RR Lyrae, Figure 7 shows that the variability information is absolutely indispensible to define a sample with a interesting combination of purity and completeness. For QSOs, (time-averaged) PS1 color together with WISE already do a very good job in selecting QSOs. The PS1 variability provides a significant, but not decisive improvement of purity and completeness. These different precision-recall curves also indicate what one might expect for purity and completeness, when a particular source lacks some information, for instance, a detection in WISE or particular PS1 bands.

Given that our training sets are finite in size, the purity and completeness will depend in detail on the chosen training sample. The individual lines in the lower panels of Figure 7 reflect different samplings of the training set. For a training set of the size available in S82, the effect is noticeable, but small.

Completeness and purity may depend on the brightness of objects under consideration. The difference between the blue solid and dashed curves in Figure 7 shows how purity and completeness change when using a bright sample versus using the entire sample. For QSOs, purity and completeness are significantly reduced, though only a little effect is evident for RR Lyrae. The validity of these conclusions, however, depends on the completeness of our training set at faint magnitudes; see Section 3.7.

To test the performance of the classifier with regard to the signal to noise of the light curves, and for RR Lyrae implicitly their distance, we considered the purity and completeness for objects classified through p_RRLyrae ≥ 0.2, as a function of their apparent magnitude (Figure 8).

Zoom In Zoom Out Reset image size

For the S82 RR Lyrae sample, using a threshold of p_RRLyrae ≥ 0.2, we get a purity = 75%, completeness = 92% within S82. The upper panel gives the dependence of purity and completeness on the apparent r_P1 band magnitude. The training set consists of the RR Lyrae in the lower panel, showing large variation in the number of sources depending on the mean r_P1 band magnitude. To account for the different number of sources in different $\langle {r}_{{\rm{P1}}}\rangle$ bins, we calculate error bars on the purity and completeness in the upper panel from the 68% confidence interval of a Poisson distribution.

According to Figure 8, we assume the purity and completeness for p_RRLyrae ≥ 0.2 to be constant within $15\leqslant \langle {r}_{{\rm{P1}}}\rangle \leqslant 20$, but might decrease or be affected by shot noise from the low number of SDSS S82 RR Lyrae beyond this interval. For the detection of faint RR Lyrae, we have to account for a loss in candidates with increasing distance, as discussed in Section 4.2.2.

Our method of automatic source classification is subject to several limitations. The most important of these are:

• (i)  
  mismatch between our ground-truth training set and other regions of sky,
• (ii)  
  incompleteness of the training set, and
• (iii)  
  the inhomogeneity of the available data over the sky.

We address these limitations in the following subsection.

We train our classifier using data on SDSS S82, where existing large catalogs of RR Lyrae and QSOs provide an almost complete sample of objects in the region. We wish to train the classifier on this region, but apply it to other regions, where no similar classifications already exist. In general, however, the application of the classifier to regions other than S82 is only justified when the region has distributions of RR Lyrae, QSOs, and potential contaminants similar to that in S82. Over most of the high latitude sky, this is the case, and we can apply our technique without difficulty.

However, at low latitudes the number of contaminants is relatively much larger than the number of RR Lyrae and QSOs in S82, since in these regions the data include very large numbers of metal-rich disk stars. Additionally, the data is itself qualitatively different: the presence of reddening changes the colors of sources, and variation in reddening as a function of distance means that even with a perfect 2D reddening map, dereddened colors may no longer match the true colors of objects. Accordingly, at low latitudes we do not expect our classifier to perform with the same purity and completeness as at high latitudes, and our S82-based estimates of purity and completeness will no longer apply.

The second problem with our technique is that even in high latitude regions, our adopted training set is imperfect. This is especially the case for our adopted QSO training set. We use spectroscopically selected QSOs from Schmidt et al. (2010), which are complete only down to roughly an i_P1-band magnitude limit of 21.25. Therefore, in our training set, fainter objects are marked as non-QSOs, and our classifier learns to discard these objects—even when they are, in fact, QSOs, as indicated by their WISE $W1-W2$ color and variability. This results in a quasar sample from our technique whose purity and completeness is really only relative to S82 spectroscopic quasars, rather than the underlying population of QSOs falling in our magnitude range.

We expect that our training set is more complete for RR Lyrae, though the very small number of distant RR Lyrae (Figure 8) means that we would run the risk there of discarding all distant objects as well, were it not for the fact that we do not include the r_P1-band magnitude as a parameter when classifying RR Lyrae (Table 3).

A final concern with our technique is that our ability to determine if an object is in fact a QSO or RR Lyrae depends on what information is available for it. The purities and completenesses we compute are properties of the entire sample of selected objects. The assignment to classes of individual objects within that sample may be relatively uncertain, if, for instance, those objects lack specific PS1 colors or detections in WISE. Figure 7 serves to show what may happen to the purity and completeness of subsamples of objects, for which only limited information is available.

QSOs should be distributed isotropically across the sky, with a mean number density of candidates, of about 20 objects per deg² in the magnitude range 15 < mag < 21.5 (Hartwick & Schade 1990; Schneider et al. 2007; Schmidt et al. 2010). This allows us to test the large scale homogeneity of our classification in areas of high Galactic latitude, and it allows us to look at the changing completeness and purity toward low latitudes. As we expect contaminants to increase at low latitudes, we expect many more candidates with low p_QSO. Until dust extinction and disk star contamination become severe, we may still get an approximately uniform density of objects with high p_QSO.

Some of these expectations are borne out in the candidate selection near the Galactic north pole: as shown in Figure 9 the selection of candidates with p_QSO ≥ 0.6, accounting for a purity in S82 of 82% and a completeness of 75%, is uniform to a high degree.

In regions away from the Galactic plane, we get a homogeneous distribution of the QSO candidates. We see this homogeneity in Figure 9 as well as in Figure 18 down to $| b| \sim 10^\circ$. When comparing the increase in the cumulative source density between p_QSO = 0.6 and p_QSO = 1, we find an increase of about 20 sources per deg² (see Figure 9(b)). The number of sources per deg² at a given minimum p_QSO is compareable for all $| b| \gt 20^\circ$, and compareable to S82. At high latitudes, the increase of candidates with p_QSO is similar on and off S82, as illustrated in Figure 9. A sample selected using a lower threshold of p_QSO shows inhomogeneities caused by contamination at almost all Galactic latitudes.

Around the Galactic anticenter (see Figure 10), the number of sources with low p_QSO per deg² is much higher than around the Galactic north pole, by a factor of ∼5. This higher overall source density does not lead to an (presumably erroneous) increase of the number of candidate objects with a high p_QSO. Indeed, the number of candidates decreases, caused by dust or varying WISE depth, to less than 2 objects per deg² (see Figure 10 (b)).

Across PS1's entire 3π area, we find 399,132 likely QSO candidates with p_QSO ≥ 0.6 (with an expected high-latitude purity = 82%, completeness = 75%), 892,131 candidates with p_QSO ≥ 0.2 (purity = 77%, completeness = 95%), and 1,596,319 possible candidates (purity = 0.72%, completeness = 98%) with p_QSO ≥ 0.05.

In this section, we present the properties of the resulting RR Lyrae candidate sample. In particular, we test whether the completeness and purity of our selection is good enough to recover known halo substructure, as well as whether it can compete with the classification from other surveys the method is not trained on.

Figures 11, 12, and 19 present the diagnostics of our RR Lyrae candidate identification, analogous to the figures for the QSO candidates. Because we expect the angular and 3D distribution of RR Lyrae to be highly structured, diagnosing the quality of our candidate identification across PS1 3π is more complex than for the QSOs. Even Figure 11, showing the RR Lyrae distribution around the Galactic north pole in its top-panel, shows gradients and structure; the overdensity seen between l = 220° and 315° is the Sagittarius (Sgr) tidal stream. The bottom panel of Figure 11 shows that we have one very likely RR Lyrae (p_RRLyrae ≥ 0.6, purity = 88%, completeness = 66%) per deg² and two per deg² with p_RRLyrae ≥ 0.05 (purity = 71%, completeness = 98%). This fits with the expectation of about 1–2 RR Lyrae per deg² from SDSS S82 (Sesar et al. 2010).

Zoom In Zoom Out Reset image size

At low latitudes, around the Galactic anticenter (see Figure 12) where the total source density is five times higher, the number of RR Lyrae candidates with p_RRLyrae ≥ 0.6 is comparable, only 1–2/deg². This may reflect the combination of higher contamination (reducing the number of p_RRLyrae ≥ 0.6 candidates), with actual RR Lyrae in the Galactic disk. The number of possible RR Lyrae (candidates), with p_RRLyrae ≥ 0.05 is much higher than around the Galactic north pole, by a factor of ∼5–10, which must reflect, foremost, increased contamination. Compared to the QSO's, we have chosen a more inclusive criterion for further consideration of RR Lyrae candidates, i.e., p_RRLyrae ≥ 0.05, because subsequent period-fitting (B. Sesar et al. 2015, in preparation) can dramatically reduce the contamination, while preserving high completeness.

The panoptic view of the PS1-selected RR Lyrae candidates (Figure 19) is quite striking, as it reveals how prominent the Galactic disk and bulge are in the map of likely RR Lyrae candidates. Note that this is in stark contrast to the large-scale distribution of probable QSOs, whose density drops toward the Galactic plane. Therefore, these data may, in addition to contaminants, be revealing enormous numbers of RR Lyrae candidates throughout the disk and the bulge. Bulge RR Lyrae have been surveyed extensively, e.g., by OGLE (Udalski 2003); yet, to date there have been very few studies of RR Lyrae throughout the Galactic disk (Mateu et al. 2012). This survey therefore represents the largest sample of Galactic disk RR Lyrae candidates, by a wide margin. Of course, they require extensive verification and follow-up.

We find 59,888 possible candidates with p_RRLyrae ≥ 0.05 (purity = 71%, completeness = 97%) and 42,674 highly likely halo RR Lyrae candidates at Galactic latitudes of $| b| \gt 20^\circ$ outside of the bulge, having p_RRLyrae ≥ 0.2 (purity = 75%, completeness = 92%).

Within $| b| \lt 20^\circ$, where reddening and contamination mean our method is less likely to be reliable (Section 3), we find 187,393 possible RR Lyrae candidates with p_RRLyrae ≥ 0.05 and 110,474 highly likely RR Lyrae candidates with p_RRLyrae ≥ 0.2. Out of them, 19,958 with p_RRLyrae ≥ 0.05 and 12,967 with p_RRLyrae ≥ 0.2 belong to the bulge as being in a radius of 20° around the Galactic center. From this, we get 167,435 possible and 97,510 highly likely disk RR Lyrae candidates outside of the bulge. Within the complete area covered by PS1 3π, we find 247,281 possible RR Lyrae candidates and 153,151 highly likely RR Lyrae candidates.

At higher Galactic latitudes, the PS1 3π includes sky regions with known halo substructures or satellite galaxies that contain RR Lyrae, and we can make use of this to verify our candidate selection. Known structures, clusters, and satellite galaxies are labeled¹² in Figure 19. Many of them show up in the our map of likely RR Lyrae. Note that M31 and M33 appear in these maps, presumably because their (unreddened) Cepheids get misintepreted as RR Lyrae by our classifier.

4.2.1. The Sagittarius Stream

The dominant substructure in the Galactic halo (aside from the Magellanic clouds) is the Sagittarius stellar stream, with its leading and trailing arms (Majewski et al. 2003). Already in Figure 19, the Sagittarius stellar stream can be seen as an overdensity crossing l = 0° and l = 180°. It is useful to show the geometry of the Sagittarius stream by selecting stars near its presumed orbital plane, and then showing a projected view of this orbital plane. Specifically, we plot the angular and distance distribution for RR Lyrae candidates with p_RRLyrae ≥ 0.2 (formal purity = 76%, completeness = 94%) in Figure 13 using the heliocentric Sagittarius (orbital plane) coordinates $({\tilde{{\rm{\Lambda }}}}_{\odot },{\tilde{B}}_{\odot })$ defined by Belokurov et al. (2014) and a distance modulus D from the mean magnitude $\langle r\rangle$. In this coordinate system, the equator is aligned with the plane of the Sagittarius trailing tail, and ${\tilde{{\rm{\Lambda }}}}_{\odot }$ increases in the direction of Sagittarius motion. The latitude axis ${\tilde{B}}_{\odot }$ points to the Galactic north pole.

Zoom In Zoom Out Reset image size

Distances D in parsec were taken from

$$\begin{eqnarray}&&D={10}^{((\langle r{\rangle }_{\mathrm{deredd}}-{M}_{r}+5)/5)}.\end{eqnarray} \tag{ 21 }$$

where $\langle r{\rangle }_{\mathrm{deredd}}$ is the dereddened r mean magnitude. We adopt M_r ∼ M_V = 0.60 mag from Sesar et al. (2010), who used the Chaboyer (1999) ${M}_{V}-[\mathrm{Fe}/{\rm{H}}]$ relation under the assumption that the mean metallicity of RRab stars is equal to the median metallicity of halo stars ([Fe/H] = −1.5, Ivezić et al. 2008). As we do not distinguish between RRab and RRc stars from our analysis, and RRab stars are most common, making up 91% of the observed RR Lyrae, we use M_r ∼ M_V = 0.60 mag for all RR Lyrae candidates.

Figure 13, showing the RR Lyrae candidates in the Sagittarius plane, provides a striking view of the stream, with its trailing and leading arm to distances of about 100 kpc. We can compare the structure in Figure 13 to Figure 6 in Belokurov et al. (2014) as well as to Figures 6 and 17 in Law & Majewski (2010a), which show the best-fit N-body debris model in a triaxial halo and observational constraints from 2MASS+SDSS for the leading and trailing arm.

We compare our results to Ruhland et al. (2011), who traced the Sagittarius stellar stream using BHB stars and compared it to Law & Majewski (2005). From our results, we can confirm that there is an extension of the trailing arm at distances of 60–80 kpc from the Sun as given, e.g., by Ruhland et al. (2011). Furthermore, we find a cloud-like overdensity at ${\tilde{{\rm{\Lambda }}}}_{\odot }\sim 110^\circ$, 5 ≲ D ≲ 25 kpc, that can be identified with the Virgo overdensity. This overdensity can be seen in a number of works (Newberg et al. 2007; Cole et al. 2008; Ruhland et al. 2011), but our RR Lyrae candidates show the three-dimensional structure especially clearly.

4.2.2. Distance Accuracy from the Draco dSph

The Draco dwarf galaxy, at known distance and known to contain many RR Lyrae, provides us with an opportunity to estimate the distance precision of the RR Lyrae candidates, using their inferred mean magnitude in the r band. As we expect many RR Lyrae in this direction, we consider an inclusive set of candidates with p_RRLyrae ≥ 0.05. Draco is entirely dominated by old stars, and is affected by near-negligible reddening, which increases the likelihood of dealing with true RR Lyrae stars as compared to the candidates seen in the region of the Galactic disk. Out of the 272 RR Lyrae listed by Kinemuchi et al. (2008), in 248 cases we found a cross-matching PS1 source within our cuts, which compares well with our 10% selection loss (Section 2.2).

Out of these, 164 have a p_RRLyrae ≥ 0.05. The first panel of Figure 14 shows their angular distribution (black points); the second panel shows their distribution in distance D. Our result of 75.3 ± 4.0 kpc is in good agreement with Kinemuchi et al. (2008), who found a distance of 82.4 ± 5.8 kpc, and Bonanos et al. (2004), who found a distance of 75.8 ± 5.4 kpc. Remarkably, the variance in our estimated distances is only ∼4 kpc, or 5% in distance. This provides us with an excellent empirical estimate of the distance precision of RR Lyrae candidates, before period-fitting (B. Sesar et al. 2015, in preparation). Note that many other satellites within ∼100 kpc also show clusters of RR Lyrae candidates (see Figure 19).

Zoom In Zoom Out Reset image size

In Figure 15 we show the heliocentric distance distribution of RR Lyrae candidates with p_RRLyrae ≥ 0.05 (purity = 71%, completeness = 97%) at Galactic latitudes $| b| \geqslant 20^\circ .$ Half of them are within 20 kpc. The most distant candidates with p_RRLyrae ≥ 0.05 are ∼150 kpc away. A profile ∼D^−1.5 related to a galactocentric halo density profile ρ ∼ D^−3.5 is overplotted for illustrative purposes. Such a halo profile is in the ball-park of recent determinations (Sesar et al. 2010; Deason et al. 2011; Xue et al. 2015). Comparing this profile to the distance distribution of our RR Lyrae candidates, we find this fits well up to ∼80 kpc. However, a rigorous derivation of the RR Lyrae density profile is beyond the scope of this paper.

Zoom In Zoom Out Reset image size

4.2.3. Comparison to the Catalina Survey

Of course, PS1 is not the first large-area RR Lyrae survey at high Galactic latitudes; so in selected areas, we can also compare with previous surveys, e.g., SDSS (York et al. 2000), Catalina (Drake et al. 2009), QUEST (Mateu et al. 2012), and PTF (Rau et al. 2000). Having used SDSS S82 in the training of the classifier, we focus here on the CSS (Drake et al. 2009), which has covered the region around the Galactic north pole down to b = 30°, but only for bright sources ≤19 mag. CSS is a survey program for finding new near-Earth objects, composed of the original CSS, the Siding Spring Survey (SSS) and the Mt. Lemmon Survey (MLS). Catalina photometry covers objects in the range −75° < δ < 70° and $| b| \gtrsim 15^\circ$. In addition to an asteroid search, the complete Catalina data is analyzed for transient sources by the Catalina Real-time Transient Survey (CRTS), resulting in catalogs of RR Lyrae (Drake et al. 2009, 2013a, 2013b). We use this to verify and check our RR Lyrae candidate identification, by cross-matching in this region with respect to CSS and SSS. We focus on the magnitude range in common between both surveys ∼15–18.5 mag.

In Figure 16, we compare our RR Lyrae candidate list with the RR Lyrae identified by CSS. For this, we cross-match sources from our RR Lyrae candidate list to those of CSS with a matching distance of 5 arcsec and keep only the source with the closest match, following position errors reported in Casetti-Dinescu et al. (2015) which we also found for CSS. Additionally, we cut to a magnitude range covered by both CSS and our analysis, 15 < V < 18.5. In this magnitude range, the total number of CSS RR Lyrae with b > 30° is 5108. For 4879 of them, we find a PS1 source within 5 arcsec. Out of these, 4686 have p_RRLyrae ≥ 0.05, 193 have p_RRLyrae < 0.05, and 229 never enter our analysis, as the PS1 source does not fulfill the selection criteria of Section 2.2.

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

With respect to CSS, we get a completeness of 92% (i.e., we find 92% of their RR Lyrae), and a cross-identified fraction of 30% (i.e., they find 30% of our sample), if we adopt the above magnitude cuts and p_RRLyrae threshold. The completeness of 92% can be explained by the 10% selection loss (Section 2.2).

When comparing to the SSS RR Lyrae, we also chose a matching tolerance of 5 arcsec and kept only the nearest match. Restricting to 15 < V < 18.5, there are 3148 RR Lyrae in the region covered both by PS1 and SSS with −30° < δ < −15°. Out of these, 2785 have p_RRLyrae ≥ 0.05, 233 have p_RRLyrae < 0.05, and 130 never enter our analysis. To assess the cross-identified fraction, we have to consider $| b| \gt 15^\circ$, as SSS roughly misses $| b| \lt 15^\circ$. The number of PS1 RR Lyrae candidates in the overlapping region and magnitude range and p_RRLyrae ≥ 0.05 not cross-matched to SSS is 11,336. The number of SSS RR Lyrae within these boundaries is 2725.

In total, this results into a completeness with regard to SSS of 88%, and a cross-identified fraction of 20%.

We find 2–3 times more RR Lyrae candidates with p_RRLyrae ≥ 0.05 than the pure samples of CSS and SSS RR Lyrae. This is in line with our assessment of purity at such lenient candidate criteria p_RRLyrae ≥ 0.05. Taken together, Considering CSS and SSS's claim of 70% completeness (Torrealba et al. 2009), our 10% selection loss as of Section 2.2, and our purity of ∼71% at p_RRLyrae ≥ 0.05, we expect about 44% of our candidates to be cross-identified in CSS or SSS; this is close to the actual fraction of 30% for CSS within 5 arcsec. In the SSS, we obtain a lower cross-identified fraction of 20%; this suggests that the completeness of the SSS is in fact lower than that of the CSS.

We have processed 3.88 × 10⁸ PS1 3π sources that fulfill the cuts described in Section 2.1. From these, we provide a catalog of all likely variable point sources in PS1 and of all likely QSOs, a total of 25,790,103 sources. We include all sources fulfilling the criterion of $\mathrm{log}{\hat{\chi }}^{2}\gt 0.5$ (see Figure 3) or W12 > 0.5. The latter criterion is intended to ensure that we provide variability statistics for almost all QSOs.

The Catalog of Variable Sources is available in its entirety in machine-readable format in the online journal. A table structure is shown in Table 4 for guidance regarding its form and content.