Variable classification in the LSST era: Exploring a model for quasi-periodic light curves
J. C. Zinn, C. S. Kochanek, S. Koz?owski, A. Udalski, M. K. Szyma?ski, I. Soszy?ski, ?. Wyrzykowski, K. Ulaczyk, R. Poleski, P. Pietrukowicz, J. Skowron, P. Mróz, M. Pawlak
aa r X i v : . [ a s t r o - ph . I M ] D ec MNRAS , 1–19 (2016) Preprint 20 July 2018 Compiled using MNRAS L A TEX style file v3.0
Variable classification in the LSST era: Exploring a modelfor quasi-periodic light curves
J. C. Zinn, ⋆ C. S. Kochanek, , S. Koz lowski, A. Udalski, M. K. Szyma´nski, I. Soszy´nski, L. Wyrzykowski, K. Ulaczyk, , R. Poleski, , P. Pietrukowicz, J. Skowron, P. Mr´oz, and M. Pawlak Department of Astronomy, The Ohio State University, 140 West 18th Avenue, Columbus OH 43210, USA Center for Cosmology and AstroParticle Physics, The Ohio State University, 191 W. Woodruff Avenue, Columbus OH 43210, USA Warsaw University Observatory, Al. Ujazdowskie 4, 00-478 Warszawa, Poland Department of Physics, University of Warwick, Gibbet Hill Road, Coventry, CV4 7AL, UK
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
LSST is expected to yield ∼ light curves over the course of its mission, which willrequire a concerted effort in automated classification. Stochastic processes provideone means of quantitatively describing variability with the potential advantage oversimple light curve statistics that the parameters may be physically meaningful. Here,we survey a large sample of periodic, quasi-periodic, and stochastic OGLE-III variablesusing the damped random walk (DRW; CARMA(1,0)) and quasi-periodic oscillation(QPO; CARMA(2,1)) stochastic process models. The QPO model is described by anamplitude, a period, and a coherence time-scale, while the DRW has only an amplitudeand a time-scale. We find that the periodic and quasi-periodic stellar variables aregenerally better described by a QPO than a DRW, while quasars are better describedby the DRW model. There are ambiguities in interpreting the QPO coherence timedue to non-sinusoidal light curve shapes, signal-to-noise, error mischaracterizations,and cadence. Higher-order implementations of the QPO model that better capturelight curve shapes are necessary for the coherence time to have its implied physicalmeaning. Independent of physical meaning, the extra parameter of the QPO modelsuccessfully distinguishes most of the classes of periodic and quasi-periodic variableswe consider. Key words: methods: data analysis – stars: oscillations – stars: variables: general –stars: variables: Cepheids – stars: variables: RR Lyrae – Magellanic Clouds
Time-domain astronomy has resulted in a range of observedastrophysical variability regimes that may be broadly cate-gorized as periodic, quasi-periodic, stochastic, and transient.Transient events are short time-scale changes in flux suchas gamma ray bursts (Klebesadel et al. 1973; Cano et al.2016), supernovae (for a review see, e.g., Woosley & Weaver1995), cataclysmic variables (Robinson 1976), and stellarflares (Variable Star Network; Kato et al. 2004). These canfrequently be modeled with a template light curve whosestructure and parameters can be used to classify the tran-sient (e.g., photometric SN classification, Pskovskii 1977;Hamuy et al. 1996; Sako et al. 2008, 2011). Periodic vari-ables such as RR Lyrae, Cepheids, and eclipsing bina- ⋆ E-mail: [email protected] (JCZ) ries are classified based on their period and the struc-ture of their phased light curves (e.g., Pojma´nski 2002;Debosscher et al. 2007; Soszy´nski et al. 2008; Sarro et al.2009; Soszy´nski et al. 2009a; Graczyk et al. 2011). Quasi-periodic sources such as Miras or spotted stars have oneor more dominant frequencies describing their variability,but do not maintain consistent phase and/or amplitude asa function of time (e.g. Howarth 1991; Bedding et al. 2005).In this regime, sources are classified by a period and aver-age light curve structure (usually Fourier component ratios,e.g., Soszy´nski et al. 2008), but there is no measure of the co-herence of the variability. Stochastic variables such as AGNvary without any obvious pattern and may be modeled bystochastic processes (e.g., Kelly et al. 2009; Koz lowski et al.2010; MacLeod et al. 2010; Zu et al. 2013).While studies of periodic sources are ubiquitous in as-tronomy, there are far fewer quantitative studies of quasi- c (cid:13) Zinn et al. periodic systems. Presently, the primary large-scale ap-plication is the determination of stellar rotation ratesfor
Kepler stars (Balaji et al. 2015; McQuillan et al. 2013;Reinhold et al. 2013). In this case, the variability is quasi-periodic due to the finite lifetime of star spots and differen-tial rotation. These periods have been used to study rota-tion rates in stars (e.g., Frasca et al. 2011), and are the basisfor gyrochronology, where stellar ages are estimated basedon models for spin-down rates (Garc´ıa et al. 2014). Quasi-periodic variables such as Long-Period variables (LPVs) arealso well known. Qualitative measures of irregularity in lightcurves have been used for classification, most notably in theGeneral Catalogue of Variable stars (Kholopov et al. 1985)and the Optical Gravitational Lensing Experiment Cata-logue of Variable Stars (Soszy´nski et al. 2008), but quan-tification of the degree of quasi-periodicity for these sourcesremains little-explored.The study of variable sources benefits from a broadrange of surveys. There are many ongoing surveys for tran-sients such as the Optical Gravitational Lensing Experi-ment (OGLE; e.g., Udalski et al. 2008), the Catalina Real-Time Transient Survey (Drake et al. 2009), the PalomarTransient Factory (Law et al. 2009), QUEST (Snyder 1998),the All-Sky Automated Survey (Pojma´nski 1997), and theAll-Sky Automated Survey for SuperNovae (Shappee et al.2014). Apart from these variability surveys, searches fortransiting exoplanets such as the Hungarian AutomatedTelescope (HAT, Bakos et al. 2002),
Kepler (Borucki et al.2010), CoRoT (Auvergne et al. 2009), KELT (Pepper et al.2004), and SuperWASP (Pollacco et al. 2006; Smith et al.2014) have also produced a wealth of light curves forother variable phenomena. More general surveys such asthe Sloan Digital Sky Survey (York et al. 2000), Pan-STARRS (Kaiser et al. 2002), and the Dark Energy Sur-vey (The Dark Energy Survey Collaboration 2005) also in-clude time domain components. Future surveys such as TheLarge Synaptic Survey Telescope (LSST; Ivezic et al. 2008;Tyson 2002; LSST Science Collaboration et al. 2009), TESS(Ricker et al. 2014), PLATO 2.0, and the Zwicky TransientFacility (Smith et al. 2014) promise to expand time domainstudies still further.Two challenges to science with variability surveys areto first classify the variable sources and then to extractphysical information about the variability from the lightcurves. A promising avenue for parsing the huge numbers oflight curves is machine learning (e.g., Wo´zniak et al. 2004;Debosscher et al. 2007; Richards et al. 2012; Bloom et al.2012; see Kim et al. 2014; Ball & Brunner 2010 for reviewsof machine learning in astronomy). In most machine learn-ing approaches, a large number of statistics are calculatedfor each light curve and an algorithm is trained to classifyvariables based on the properties of a sample of previouslyclassified variables (but see, e.g., Mackenzie et al. 2016, foran example of unsupervised classification of variable stars).While these methods work well for identifying transientssuch as supernovae and highly periodic sources like Cepheidsand RR Lyrae, stars with irregular periodicity and multi-periodicity have proven more difficult to accurately classify.One problem may be that the statistics employed in many ofthe current approaches have no physical significance. Param-eters with physical meaning such as periods may be betterfor classifying light curves than ad hoc statistical measures, in addition to providing insights into the nature of the vari-able source. The challenge is the definition and extractionof new physically meaningful parameters from quasi-periodicor stochastic light curves.The introduction of stochastic process models forquasar variability (Kelly et al. 2009; Koz lowski et al. 2010)is a good recent demonstration of these principles. Adamped random walk (DRW) stochastic process charac-terized by an amplitude and a coherence time-scale pro-vides a good, compact statistical model of quasar vari-ables on time-scales of days to decades (MacLeod et al.2012). The parameters of the model are then correlatedwith the observed wavelength and the intrinsic proper-ties of the quasars (e.g., MacLeod et al. 2010), potentiallyproviding insight into the physical origin of the variabil-ity (Kelly et al. 2009; Dexter & Agol 2011; Andrae et al.2013). There is evidence for deviations from the DRW, par-ticularly on short time-scales (e.g., Mushotzky et al. 2011;Zu et al. 2013; Graham et al. 2014; Andrae et al. 2013;Kasliwal et al. 2015), and this may allow the extraction ofadditional physical scales from AGN light curves. The DRWmodel has also seen application in He et al. (2016) to modelMiras, in which they find that adding a DRW componentin addition to a purely periodic component can improve thefidelity of period estimates by up to 20 per cent.The DRW model is the first ( p = 1) in a series of con-tinuous time autoregressive moving average (CARMA( p, r ))processes, which are the solutions to stochastic linear differ-ential equations of order p driven by white noise and up to r < p derivatives of white noise processes. The higher or-der processes in this sequence have occasionally been usedin astronomy (e.g., Koen 2005, 2012; Andrae et al. 2013;Kelly et al. 2014; Aigrain et al. 2016; Kasliwal et al. 2015).In particular, Kelly et al. (2014) offer a package implement-ing fits of these models to light curves based on meth-ods from the forecasting literature (e.g., Brockwell & Davis2002).Physically, the p = 2 process models a quasi-periodicoscillation (QPO) characterized by an amplitude, a period,and a coherence time. Koen (2005, 2012), Kelly et al. (2014),and Aigrain et al. (2016) give examples of applying themodel to periodic and quasi-periodic variables, but it hasnever been employed for a large-scale variability survey. Herewe carry out such a survey, applying the QPO model to largenumbers of variable stars and quasars in the OGLE-III fieldsfor the Large Magellanic Cloud (LMC). We also compare theQPO results to those from modeling the same sources usinga DRW, and explore to what extent certain variables pre-fer a QPO to traditional purely periodic models. The goalsare to explore whether the two models successfully distin-guish between variable stars and quasars and whether theQPO coherence time provides a useful means of classify-ing and understanding variable stars. In § § § § The Optical Gravitational Lensing Experiment (OGLE) ob-serves the Magellanic Clouds, and portions of the Galac-
MNRAS , 1–19 (2016) xploring a model for quasi-periodic light curves tic bulge and disk in the I- and V-bands using the 1.3-mWarsaw telescope at the Las Campanas Observatory. Weuse light curves of sources in or behind the LMC fromthe OGLE-III survey (Udalski et al. 2008), which operatedfrom 2001 to 2009. We examine 3361 classical Cepheids(Soszy´nski et al. 2008), a subset of 50000 Long PeriodVariables (Soszy´nski et al. 2009b), and 24906 RR Lyrae(Soszy´nski et al. 2009a), as well as well as OGLE-III/IVlight curves of 753 QSOs identified in the Magellanic QuasarsSurvey (MQS, Koz lowski et al. 2013). We used only the I-band data as they have a better mean cadence ( ∼ ∼ I = 19 . σ ′ i = p ( γσ i ) + ǫ , to match the observed dis-persion in standard stars as a function of magnitude, where γ and ǫ are computed on a field-to-field basis as describedin (Wyrzykowski et al. 2009).The OGLE LPVs have been broadly categorized basedon light curve properties into Semi-Regular Variables, Mi-ras, and OGLE Small-Amplitude Semi-Regular Variables(OSARGs) by Soszy´nski et al. (2009b). All three types areunderstood to be giant branch or asymptotic giant branchstars. OSARGs have variability amplitudes ranging from0.005 to 0.13 mag with periods ranging from 10 to 100 days(Wray et al. 2004), and the pulsation mechanisms responsi-ble for their light curves are not well understood (but seeSoszy´nski et al. 2004; Christensen-Dalsgaard et al. 2001).Miras are generally believed to be variable as a result of in-herently unstable radial modes due to opacity effects in theatmosphere (e.g., Wood 1974), with periods ranging from100 to 1000 days and amplitudes of up to 10 mag. Semi-Regular Variables (SRVs) have similar periods to Miras, butare distinguished from Miras by their lower amplitudes andbecause they pulsate in both the fundamental and the over-tone radial modes, instead of only the fundamental mode.Although their light curves exhibit qualitatively differentlight curves, with SRVs generally having less coherent peri-odicity than Miras, Soszy´nski et al. (2013) hypothesize thatas red giants evolve, they first pulsate as OSARGs, thenas SRVs, and finally as Miras (see also Bedding & Zijlstra1998).Among the best-studied variable stars are Cepheidvariables—giant stars with pulsation periods mostly from 2to 20 days and amplitudes from 0.2 to 1 mag in V-band dueto the κ -mechanism (Zhevakin 1953). The OGLE Cepheidsare divided into sub-classes based on their pulsation modes,from fundamental mode (“F”) and first overtone (“1”) pul-sators up to third overtone pulsators, with a small numberof multi-mode pulsators (Soszy´nski et al. 2008). RR Lyraestars are horizontal branch stars in the instability strip that also pulsate due to the κ -mechanism (Christy 1966).Based on their periods and light curve shapes, RR Lyrae aresubdivided into RRab (fundamental mode pulsators), RRc(first overtone pulsators), RRe (second overtone pulsators)and RRd (simultaneous fundamental and first overtone pul-sators). Typical variability amplitudes for RR Lyrae rangefrom 0.2 mag to 2 mag in V-band, and the pulsation periodsrange from of 0.2 to 1 days.We model light curves using a family of stochastic pro-cess models known as Continuous time Autoregressive Mov-ing Averages (CARMA). These models statistically describethe solution to stochastically-driven ODEs by way of an au-tocovariance function, S ≡ h ~s | ~s i , whose structure dependson the order of the ODE, p , and the weights given to the r < p noise processes driving the system. For the purposesof this paper, we call the p = 1 case a damped randomwalk (DRW) and the p = 2 case a quasi-periodic oscilla-tion (QPO). The autocovariance functions of the CARMAprocesses are the sum of p exponentials e ω i t , where the ω i must be real or appear in complex conjugate pairs (seeBrockwell & Davis 1996). The DRW model is the lowest or-der example, and its covariance function between times t i and t j is S ij = σ DRW exp ( −| t j − t i | /τ ) , (1)where σ DRW describes the variance of the light curve on longtime-scales and τ is the coherence time. We use the DRW forcomparison to the QPO model in order to explore whether aquasi-periodic model is a better model for the stellar variablesources (or the reverse for quasars). The QPO model is thecase p = 2, with a complex conjugate pair of w i , whichwe can divide into an oscillation period, P = 2 π/ω , and acoherence time, τ , to give the autocovariance function S ij = σ −| ( t j − t i ) | /τ ] cos [ ω ( t j − t i )] . (2)We could also add a phase, but we omit a phase factor infavour of assuming that all real variable sources have a max-imum in their autocovariance functions as | t j − t i | → τ parameter indicates the timescaleover which variability is correlated, the QPO τ parameter isthe time-scale over which the process appears as a coherentoscillation. Thus, the ratio of the coherence time to the pe-riod, τ /P , roughly corresponds to the number of oscillationsover which the light curve will maintain phase coherence.Figure 1 illustrates this meaning of the coherence time witha series of phased light curves ranging from an essentiallyperiodic Cepheid at the top to a non-periodic MQS quasarat the bottom. We adopt τ /P > as our operational defi-nition of periodic, although RR Lyrae and Cepheid variablescan be even more strictly periodic. In practice, we cannotdistinguish among models with τ > days due to theOGLE time baseline (cf., the discussion of Figure 2, below;see also Koz lowski 2016). When fitting models, we impose τ < days. In the limit that τ /P ≪
1, illustrated by thequasar, the QPO model becomes the DRW model.In addition to the process model for the variability, weassume the data have uncorrelated Gaussian errors (noise)with an autocovariance function N ≡ h ~n~n i = δ ij σ i . Thedata are then described by the total covariance matrix, C = S + N , which combines the covariance due to the signal S andthe covariance due to the noise, N . The model also includes MNRAS000
1, illustrated by thequasar, the QPO model becomes the DRW model.In addition to the process model for the variability, weassume the data have uncorrelated Gaussian errors (noise)with an autocovariance function N ≡ h ~n~n i = δ ij σ i . Thedata are then described by the total covariance matrix, C = S + N , which combines the covariance due to the signal S andthe covariance due to the noise, N . The model also includes MNRAS000 , 1–19 (2016)
Zinn et al. −0.10.00.10.20.30.4
OGLE-LMC-CEP-0069 τ/P = 5532.01 −0.50.00.51.0
OGLE-LMC-RRLYR-00001 τ/P = 38.62 −1.0−0.50.00.51.01.52.02.5
OGLE-LMC-LPV-00082 τ/P = 5.18 −20246
OGLE-LMC-LPV-06171 τ/P = 2.09 −2−1012345
OGLE-LMC-LPV-18816 τ/P = 1.29 phas −0.50.00.51.0
MQS J043801.81-692014.1 τ/P = 0.17 I [ m a g ] Figure 1.
Folded light curves illustrating the meaning of the coherence time-scale in the QPO model. Each panel shows a folded lightcurve of a source ranging from coherent ( τ/P ≫
1) sources at the top to incoherent ( τ/P ≪
1) sources at the bottom. The most coherentsource shown is the Cepheid OGLE-LMC-CEP-0069, followed by an RR Lyrae, a series of LPVs, and ending with the quasar MQSJ043801.81-692014.1. That τ/P is not larger for OGLE-LMC-RRLYR-00001 is discussed in § a set of linear parameters to deal with the light curve meansand long time scale trends (see the Appendix). In all caseswe use the linear parameters to remove the light curve mean.As we explore in §
3, long term trends can affect the resultsfor the LPVs. We experimented with removing polynomialtrends up to third order, ultimately settling on using a cubictrend model for all stellar sources (LPVs, RR Lyrae, andCepheids), and two trends for QSOs corresponding to theOGLE III and IV portions of the light curves because theycan be offset. We only used the OGLE III light curves forthe variable stars so there was no need to allow for meanoffset other than for the quasars. Our method automaticallyincludes the uncertainties in the linear parameters into theuncertainties in the process parameters.The final likelihood function in this framework is then L = exp (cid:0) − ( ~y − L ˆ q ) T C − ( ~y − L ˆ q ) / (cid:1) (2 π ) P/ | C | / | L T C − L | / , (3)where L ˆ q is the best-fitting trend model described by thematrix, L , and the best-fitting linear coefficients, ˆ q (seethe Appendix), ~y is the observed light curve, and P is the number of parameters in the stochastic process model. Thecomputational barrier to performing such an analysis lies inthe inversion of the covariance matrix, C . Rybicki & Press(1992) found a method to reduce the scaling from O ( N )to O ( N ) for the cases p = 1 and p = 2, and this was themethod used by Koz lowski et al. (2010) and MacLeod et al.(2010) to rapidly analyse OGLE-III and SDSS light curveswith the DRW model. Ambikasaran et al. (2014) has sub-sequently found an algorithm which is O ( Np ), for all p ,and we use this sparse matrix method. The forecasting al-gorithms used by Kelly et al. (2009, 2014) are also O ( N ),but in Koz lowski et al. (2010) we found that evaluating thefull likelihood given by Equation 3 had significantly morestatistical power than the forecasting methods.We model the light curves by first finding best-fittingparameters that maximize the likelihood (Equation 3) usinga simulated annealing method, as implemented in PyGMO .Our problem necessitates a MLE approach for optimization https://github.com/esa/pagmo/tree/master/PyGMO MNRAS , 1–19 (2016) xploring a model for quasi-periodic light curves as opposed to a grid search because the length of the OGLElight curves requires large numbers of period samples evenbefore adding extra parameters. For a fixed phase error, δφ ,the period has to be sampled by ∆ P/P ∼ ( δφ/ π )( P/T ),where T ∼ P ∼ P ∼ . δφ ∼
1, this requires ∆
P/P ∼ / / emcee (Foreman-Mackey et al. 2013). We use standard logarithmicpriors for the process parameters. Our chain length is signifi-cantly larger than the autocorrelation time (e.g., Sokal 1996)of the estimated parameters calculated by emcee for mostlight curves. In this limit, the mean marginalized parameterestimated is representative of the true mean of the parame-ter. Parameter errors are quoted at 1- σ , and the best-fittingparameter is taken to be the mean of the MCMC chains.To compare models with differing numbers of param-eters, we use the Akaike Information Criterion (AIC) andthe Bayesian Information Criterion (BIC). For a light curvewith likelihood L , N data points, and k model parameters,the Bayesian likelihoods (Equation 3) are modified to be AIC ≡ − L + 2 k (4)and BIC ≡ − L + k ln N. (5)The BIC penalizes the addition of new parameters moreheavily than does the AIC. These allow us to compare mod-els with differing numbers of parameters, and, in particular,address the questions of which of the models (DRW or QPO)better fits the data, and whether a quasi-periodic model isa better fit than a periodic model. We present the AIC re-sults. Our initial criterion for rejecting sources better fit bywhite noise (the limit τ → AIC < − f ( x ) = sinh − (cid:20) AIC − AIC (cid:21) = sinh − x, (6)This transformation provides a linear scaling for | x | < ∼ | x | for | x | >>
1. Thenormalization is chosen such that the likelihood ratio is 95per cent at f = ± χ = 4 for χ statis-tics). For f < − f > | f | <
1, the maximum likelihoodprobability for model 1 (model 2) is greater than 5 per cent.The probability of the less favoured model plummets outsidethe | f | < .
07 per cent for | f | > We next discuss a series of experiments to help understandthe global results. We start by simply confirming that theprocess is internally consistent by generating QPO modellight curves with cadences, amplitudes and uncertaintiestypical of the OGLE data and showing that we can re-cover the input parameters reasonably well. We also show that we can recover the periods of the LPVs with some de-pendence on the treatment of long time scale trends in thelight curves. For shorter period variables, there is essentiallynever a problem in estimating the period. When we thenapplied the method to the actual light curves, we foundthat the results showed a significant dependence on the lightcurve shape, behaving as expected for relatively sinusoidallight curves but with counter-intuitive results for very non-sinusoidal variables like long period Cepheids. In particular,it was surprising to find that QPO models of extremely pe-riodic RRab and long period Cepheids frequently preferredlow coherence QPO models or even the DRW model. Thismotivated a series of experiments modeling purely periodicsine and triangle waves as a function of the noise relativeto the amplitude that demonstrate the effect of light curveshape on QPO parameter estimation. Finally, we introducea higher order, p = 4, QPO model to confirm that a higherorder model that can better mimic the structure of less sinu-soidal light curves provides a means of addressing the prob-lem. To provide a guide to interpreting the results in §
4, wefirst tested how well we can recover QPO parameters simplydue to the modeling process itself. We generated a sampleof 12000 artificial QPO light curves using cadences typicalof the OGLE light curves and with amplitudes, periods, andcoherence time scales covering the range of the variables westudy in §
4, with periods varying from 0 . days, τ varying from 0.1 to 10 days, and amplitudes varying from0 .
01 to 1 mag. We used a fixed photometric uncertainty of0 .
02 mag, which is roughly twice the 95 th percentile of thephotometric error distribution of the OGLE light curves. Wegenerated the light curves using the Cholesky decompositionmethod outlined in Zu et al. (2013). If the covariance ma-trix, C , can be Cholesky decomposed into C = MM T , and R is a vector of Gaussian random deviates of unit amplitude,a light curve realization is simply y = MR .Figure 2 illustrates the results. Figure 2a shows thatthe fractional errors in the periods are well modeled by | ∆ P | /P = 1 / [ B + ( τ /P ) C ] with B ≃
14 and C ≃ . τ /T . .
1, then the median fractional uncertaintyin τ is ≈
10 per cent, with a scatter of ± . τ becomes large as it approachesthe light curve timespan, T , generally in the sense that itcan be greatly over-estimated. That τ is poorly recoveredfor τ /T & τ = 10 days ( τ ≈ T ) is the best model.We also demonstrate in Figure 2b that τ is poorly recoveredfor τ /T . − , which is the typical sampling time-scale( ∼ τ is comparableto T (MacLeod et al. 2010; Koz lowski 2016).We also tested the ability of the QPO models to recoverthe published OGLE-III periods of the variable stars. Wefind that we can always recover the period of short-periodvariables, with periods .
100 days. Care must be taken,however, to adequately resolve period space when searchingfor periods much less than the baseline of the light curve( T ∼ MNRAS , 1–19 (2016)
Zinn et al. -3 -2 -1 coherence τ in /P in -7 -6 -5 -4 -3 -2 -1 | P o u t − P i n | / P i n -4 -3 -2 -1 τ in /T -4 -3 -2 -1 τ o u t / τ i n Figure 2.
Fractional error in period, | P out − P in | /P in as a function of the coherence of the oscillation, τ in /P in , for artificial light curves(left). A simple model for the fractional error in period as a function of the coherence is shown as a solid line to guide the eye. Recoveredversus input coherence time, τ out /τ in , as a function of the coherence time in units of the light curve baseline, τ in /T , where T ≈ τ in /T = 1 × − ) areincreasingly difficult to recover. lihood peaks become increasingly narrow due to the largenumber of periods spanned by the data.As shown in Figure 3 for the Miras, the situation is quitedifferent for the LPVs. The period agreement is very goodfor shorter periods ( P <
300 days), but there are significantnumbers of cases where the QPO period estimate is signif-icantly longer than the OGLE period. It is not surprisingthat the discrepancies are concentrated at longer periods,since even for truly periodic sources the fractional perioduncertainty will be on the order of the period divided by thetime span of the data. However, the discrepancies are largerthan we would expect based on our simulations of model-ing randomly generated QPO light curves (see Figure 2b).Figure 4 shows the light curves and QPO models for twoexamples, LPV-00144 where the period discrepancy is large( P OGLE = 364 versus P QPO = 577 days) and LPV-00743where the period discrepancy is smaller ( P OGLE = 216 ver-sus P QPO = 209 days). The obvious distinction is that LPV-00743 lacks the large amplitude changes of LPV-00144. Suchlarge scale features are common among the period outliers.This suggested that the period differences might in partbe driven by the cubic model for long term trends that isincluded in our standard analysis, so we repeated the cal-culations removing constant (light curve mean), linear andquadratic trends. Figure 3 compares the results for simplyremoving a constant to the default cubic model, and we seethat the number of outliers increases significantly. We alsoshow in Figure 3 where outliers for the two trend models liein the distribution found for the other trend model. Giventhe results of these experiments, we kept the cubic modelas our default for stellar sources. The trend model choice is also important for the recovery of DRW parameters, par-ticularly for QSOs. Unsurprisingly, fitting trend models ofhigher order than a simple mean results in removing genuinestochastic behavior and biases the estimates of τ toward ar-tificially small or high values (not shown). For this reason,we opted to only remove a mean each for the OGLE-III andOGLE-IV sections of the QSO light curves for the DRWmodel. A cubic was removed for each of the two light curvesections when fitting QPO models to QSOs.One technical point to emphasize about Figure 4 is thatthe QPO model curve is not a particular “best-fitting” lightcurve, but the statistical average of light curves consistentwith the data (see the discussion in Rybicki & Press 1992).The “error snake” surrounding this curve is the 1 σ dispersionof these light curves around their mean. A particular statisti-cal realization of a light curve consistent with the data wouldtrack the mean with deviations statistically bounded bythe “error snake”, and the construction of such constrainedlight curve realizations is also described in Rybicki & Press(1992). These constrained realizations will show more smallscale structure than the mean light curve, and this is a largepart of the variance captured by the “error snake”.As we developed our analysis, we found (see §
4) thatperiodic variables with light curves that strongly deviatedfrom sinusoids (e.g., fundamental mode Cepheids and ab-type RR Lyrae) were frequently either favouring QPO mod-els with short coherence times or even the DRW model overthe coherent, periodic limit of the QPO model. Figure 1 illus-trates the issue for the RRab variable OGLE-LMC-RRLYR-00001, which is clearly a coherent periodic variable over theOGLE baseline, but has a best-fitting QPO coherence time
MNRAS , 1–19 (2016) xploring a model for quasi-periodic light curves l o g O G L E P e r i o d [ d a y ] l o g O G L E P e r i o d [ d a y ] Figure 3.
A comparison of QPO and OGLE period estimates for Miras using either cubic (left) or constant (right) temporal trendmodels. The diagonal line represents a one-to-one correspondence and the error bars along the bottom indicate the median QPO perioduncertainties for a series of logarithmic bins in period. The red points in each panel are the 3 σ period outliers in the other panel. Ingeneral the cubic trend model provides better agreement and is adopted as our standard model. HJD-2450000 [day] I [ m a g ] P OGLE = 364.20P
QPO = 577.25
OGLE-LMC-LPV-00144
HJD-2450000 [day] I [ m a g ] P OGLE = 216.10P
QPO = 208.69
OGLE-LMC-LPV-00743
Figure 4.
Miras with large (left) and small (right) period differences between the OGLE and the QPO models. The points are theOGLE data and the curve is the average of QPO models at the indicated, fixed P QPO that fit the data well. The green “error snake”represents the 1 σ dispersion of these light curve models around the mean. Median photometric errors are shown as a single vertical bar.The periods are listed in each panel along with lines to graphically show the scale of the periods and their differences relative to the lightcurves. The divergence of the light curve model in the left panel is a consequence of the cubic trend model.MNRAS , 1–19 (2016) Zinn et al. of τ ∼
30 days for a period of P ≃ .
64 days so that τ /P ∼
40 instead of τ /P ∼ ∞ .To test the hypothesis that light curve shape could biascoherence time estimates and model preference, we mod-eled purely periodic sinusoid and triangle wave variability.Artificial light curves were generated using the same I-bandmagnitude (18 . σ i = 0 . τ of rescaling the reported noiseby a factor of f err , meaning that we analyse light curveswith the error defined to be f err σ i instead of the true error, σ i . The light curve itself is not modified. Under-estimatingthe true error should drive the solution towards τ → f err = 1 . . f err ≈ . χ /dof (left axis) as a function of f err . Model preference isdetermined by the AIC criterion and the scaling defined inEquation 6, for which a positive (negative) value indicatesa preference for the periodic (QPO) model. The sine, saw-tooth, and OGLE-LMC-RRLYR-00001 light curves all con-verge to a periodic model for a sufficiently large f err , essen-tially tracking the τ parameter estimates. The goodness of fitshows two branches. If the errors are greatly under-estimated( f err → τ → χ /dof <
1) andthe process can be periodic ( τ → ∞ days). For the sinu-soid, we see that for f err = 1, we have χ /dof = 1 and τ → days, as we would expect. For the sawtooth, the fitat f err = 1 is poor ( χ /dof ≃ . τ ≈
170 days), and the QPOmodel is favoured even though the light curve is periodic.If we increase f err to f err ≈ .
4, we can make χ /dof ≃ τ → days, and have a preference for the periodic model.The RR Lyra behaves similarly to the sawtooth, but wouldneed f err ≃ χ /dof = 1. The non-sinusoidallight curve shape requires the assumption of additional whitenoise (higher values of f err ) for the preferred solution to beperiodic.In the second experiment, we generated light curveswith different ratios σ/ A between the errors and the peak-to-peak variability amplitude, 2 A . Figures 5c and 5d showthe best-fitting coherence ratio, τ /P , the goodness of fit χ /dof , and the degree of preference for the QPO or pe-riodic models for this case. The sinusoidal model behavesexactly as expected, always finding the limit τ → days,with χ /dof ≃
1, and an overwhelming preference for theperiodic solution. The sawtooth model, on the other hand, prefers low coherence QPO models at high signal-to-noiseand then switches to preferring the periodic solution at lowsignal-to-noise. In the transition between these two limits atintermediate signal-to-noise, the fits to the data are poor,with χ /dof significantly greater than unity. OGLE-LMC-RRLYR-00001, whose parameters are indicated with an ‘x’,behaves like a sawtooth light curve.We can explain these behaviors by thinking aboutthe problem in terms of the power spectrum. The QPOmodel produces a Lorentzian power spectrum P ( ω ) = σ / √ π P ± (cid:2) τ − + ( ω ± ω ) (cid:3) − for a period of P = 2 π/ω and a damping time scale τ . A sinusoidal light curve becomesa pair of delta functions at ± ω . A white noise spectrum issimply a constant independent of ω . The spectrum of a non-sinusoidal periodic variable is a series of delta functions atfrequencies ± nω , n = 1 , , ... , with the envelope of ampli-tudes representing the Fourier series representation of thelight curve shape.When we fit a periodic triangle wave in the periodiclimit of the QPO, we can match only the first, n = 1 peakat | ω | = ω in the power spectrum, leaving all the power inthe n > n > τ . Suppose,for example, that the ratio of the power in the n = 1 and n = 2 peaks is R , then the QPO model can have the samepower at the second peak if τ ∼ P R − / . So the more thelight curve shape deviates from truly sinusoidal, the morethe best-fitting model will prefer a finite τ even though theunderlying model is actually periodic.This strongly suggests testing a higher order QPOmodel with the second period fixed to be P/ n = 2 term). Periodic variables, such as Cepheids andRR Lyrae, with strongly non-sinusoidal light curves, canbe well fit by Fourier series with period-dependent am-plitudes (see, e.g., Pejcha & Kochanek 2012), and there-fore should be better fit by a QPO model with a n = 2term. A simple second order example of such a light curve, f ( t ) = A cos(2 πt/P ) + B cos(4 πt/P + φ ) , corresponds to theperiodic limit of a higher order ( p = 4) stochastic processwith the periods fixed in a 2:1 ratio. In § p = 2) modelto the higher-order p = 4 QPO model results for RR Lyraeand Cepheids with the primary period fixed to the primaryOGLE period and the secondary Fourier period fixed to halfthis period, fitting only an amplitude of oscillation for theprimary and secondary Fourier period, and a shared coher-ence time ( A , B , and τ ). The phase φ does not affect theautocorrelation function. Keeping these issues in mind, we now examine the resultsfrom fitting all the variables with the DRW and QPO mod-els. We first explore whether the QPO model is a betterdescription of periodic and quasi-periodic phenomena thanthe DRW model. Next we examine how the QPO model dis-
MNRAS , 1–19 (2016) xploring a model for quasi-periodic light curves τ / P a)a)a) sinesawtoothOGLE-LMC-RRLYR-00001 τ / P c)c) sinesawtoothOGLE-LMC-RRLYR-00001 f err χ / d o f b)b)b) sinesawtoothOGLE-LMC-RRLYR-00001 σ/2A d)d) sinesawtoothOGLE-LMC-RRLYR-00001 −4−2024 s i nh − (cid:0) A I C Q P O − A I C P e r i o d i c (cid:1) −4−2024 s i nh − (cid:0) A I C Q P O − A I C P e r i o d i c (cid:1) Figure 5.
The upper panels show the coherence ratio, τ/P , for the OGLE RR Lyrae OGLE-LMC-RRLYR-00001 and artificial lightcurves with sinusoidal or sawtooth waveforms. The lower panels show the QPO model preference (Equation 6; dotted lines, right scale)and the χ /dof (solid lines, left scale); positive/negative values favour periodic/QPO. The vertical dashed line corresponds to therecommended factor by which the reported OGLE-III errors for OGLE-LMC-RRLYR-00001 should be increased (Wyrzykowski et al.2009). The horizontal dashed line corresponds to a χ /dof of unity and the grey band indicates the expected 68 per cent confidenceregion assuming Gaussian statistics. The left columns show these quantities as a function of the error scaling factor f err (error bars σ → f err σ ). The right column shows these quantities as a function of the photometric error relative to the peak-to-peak amplitude, σ i / A . A point marks the position of the χ /dof for OGLE-LMC-RRLYR-00001, and a square marks its model preference. tinguishes between periodic and quasi-periodic phenomena.Finally, we examine the QPO parameter distributions for thevarious variable classes with an eye towards classification. Figure 6 shows the distribution of the variable classes in therelative probabilities of the QPO and DRW models for theAIC likelihoods. We use the sinh − x scaling of Equation 6to compress the full dynamic range into the figure – when | x | = 1 the associated model is favoured at 95 per centconfidence.Quasar variability is generally believed to be stochas-tic rather than quasi-periodic, and we find that a major-ity of the quasars (58 per cent) strongly favour the DRW model ( x > MNRAS , 1–19 (2016) Zinn et al. −15 −10 −5 0 5 10 15 sinh −1 (cid:0) AIC
QPO − AIC
DRW (cid:1) i n t e g r a l f r a c t i o n Q P O D R W Q P O D R W Q P O D R W Q P O D R W Q P O D R W Q P O D R W MirasQSOsSRVsOSARGsF Cepheids1 CepheidsRRab+RRdRRc+RRe
Figure 6.
The integral distribution of the labeled variable classes in sinh − (cid:2) (AIC QPO − AIC
DRW ) / (cid:3) , where the QPO (DRW) modelis favoured for negative (positive) values. The sinh − x mapping is linear for x ≪ x , so that the fulldynamic range of the likelihood ratio can be displayed. The mapping is normalized so that the two vertical lines at | x | = 1 correspond to95 per cent confidence likelihood ratios in favour of the QPO ( x = −
1) or DRW ( x = 1) models. The BIC distributions are very similar. model is illusory because the best-fitting QPO model is also“incoherent” in the sense that τ /P <
1. Essentially a cubictrend model plus a long period, relatively incoherent oscil-lation can fit slow quasar variability well in such cases. If weadd the requirement that the QPO model must not only bepreferred by the information criterion, but must also havea best-fitting QPO model with τ /P >
1, then the fractionof “quasi-periodic” quasars drops below 20 per cent. Evenamong these objects, however, half have a coherence timeless than the OGLE cadence. In this regime, the QPO modeldiffers little from the white-noise limit of a DRW model: theperiod is meaningless because the coherence time is so shortthat its effect is to broaden the photometric errors. Thesesolutions are found in spite of evident variability on timescales greater than the cadence. Essentially the variability islargely modeled by the cubic trend, and the QPO is actingas an additional source of white noise (see, e.g., Figure 7).As one would expect, most of the variable stars in Fig-ure 6 are better fit by the QPO model. RR Lyrae andCepheids overwhelmingly show a strong preference ( f < − .
02 per cent) and Cepheids (0 . f >
1) for DRW. Only 0 . . τ /P < § τ /P <
1. In the limit of incoherent oscillations,the QPO model increasingly resembles the DRW model, andthe fits become indistinguishable by their likelihoods and χ /dof . In this limit, the AIC model preference will alwaysfavour a DRW model, because it has fewer parameters. Al-though periods may remain detectable in a periodogram formarginally incoherent sources, the DRW model can still bestatistically preferred in the AIC formalism due to its fewer MNRAS , 1–19 (2016) xploring a model for quasi-periodic light curves HJD-2450000 [day] I [ m a g ] J044037.32-690321.7
HJD-2450000 [day] I [ m a g ] J043522.69-690352.9
HJD-2450000 [day] I [ m a g ] J044122.61-683452.2
HJD-2450000 [day] I [ m a g ] MQS J044033.48-683825.2
Figure 7.
Four examples of QSOs for which the QPO model was favoured over the DRW model. The points are the OGLE data and theerror bar indicates the median photometric error. The offset curves show the cubic trend models. The curve and shaded regions againcorrespond to the average of the best-fitting light curves and their dispersion (see § degrees of freedom. An issue to keep in mind for the OS-ARGs is that their variability amplitudes are comparable totheir photometric noise, which makes it difficult to distin-guish among models.Figure 8 reprises the similar figure in Koz lowski et al.(2010), comparing the distribution of quasars and vari-able stars in DRW parameters. If no other criteria areimposed (left panel), the selection region proposed byKoz lowski et al. (2010) contains 62 per cent of the quasars,but also contains 0.52 variable stars for every quasar.In practice, Koz lowski et al. (2010) considered addi-tional simple selection criteria (e.g., color and magnitude)that significantly increase the efficiency of the selection pro-cess, but here we only consider the time domain criteria. Ifwe impose the criterion that the light curves must also preferthe DRW model, so AIC
QPO − AIC
DRW >
0, the contami-nation is greatly reduced (the right panel of Figure 8). Theselection region now contains 55 per cent of the quasars, but only 0.20 variable stars for every quasar. The additional cri-terion loses only 7 per cent of the quasars but reduces thestellar contamination by 67 per cent. We attempted to in-crease the efficiency of quasar selection by also selecting ob-jects that had incoherent QPO models. However, this doesnot improve quasar selection because the largest contami-nant are OSARGs, which also favour incoherent QPO mod-els. The MQS quasars are, of course, selected from a regionof very high stellar density (the Magellanic Clouds!), whichmakes the density of variable stars tremendously larger thanin a typical extragalactic field ( ∼ stars in theMagellanic Clouds!). In a true extragalactic field, the purityof such a variability-selected sample would be far higher.Additional criteria such as the magnitude and color criteriaconsidered by Koz lowski et al. (2010) could also be devel-oped to avoid the issues leading to quasars favouring theQPO model. For example, we could also require more vari- MNRAS , 1–19 (2016) Zinn et al. -1 ˆσ [mag year −1/2 ] τ [ d a y ] MirasSRVs OSARGsRR Lyrae QSOs -1 ˆσ [mag year −1/2 ] MirasSRVs OSARGsRR Lyrae QSOs
Figure 8.
Distribution of several variable classes and QSOs in the DRW parameter space of ˆ σ ≡ σ DRW p /τ and τ to match theequivalent figures in Koz lowski et al. (2010). The left panel shows all objects, while the right panel shows the distribution after eliminatingsources that prefer ( x < τ , the optimal selection region for QSOs will vary based on the time baseline of a typicallight curve. For our purposes, we retain the selection region based on OGLE-III time baselines. Imposing the cut on the relative likelihoodsof the DRW and QPO models eliminates ∼ ability “power” in the stochastic process than in the cubicpolynomial or require a minimum variability power. We next test whether sources are better fit as truly peri-odic (i.e., τ → ∞ days) or as quasi-periodic (i.e., finite τ )variables by comparing the likelihoods for a QPO with fixed τ /P = 10 (as a proxy for τ → ∞ days) to those with τ freeto vary, as shown in Figure 9. Quasars, as we would expect,strongly favour the QPO model, with 97 per cent prefer-ring ( x <
0) the QPO model. Among the few exceptions,two have extreme outliers in the light curves (up to 25 σ )and one has an exceptionally short light curve of 74 epochs.As discussed earlier, when quasars are modeled as QPOsthey strongly favour either the “incoherent” limit ( τ /P < τ is smaller thanthe cadence.It is not surprising that the LPVs shown in Figure 9generally favour a QPO over a periodic model, as thesevariables generally do not have simple, coherent light curvestructures. Only 25 per cent of the Miras strongly prefer( x >
1) a periodic model, and most of the Miras that pre-fer a periodic model have genuinely regular, sinusoidal lightcurves with mildly variable ( ∼
10 per cent) amplitudes. Theless regular SRVs show a smaller periodic fraction while thelow-amplitude OSARGs show an overwhelming preferencefor the QPO model. The low amplitudes of the OSARGsalso lead to noisier light curves, where it becomes difficultto distinguish the models (see Figure 5 and the discussionin § ∼
20 per cent of Miras and less than 1 per cent of SRVsand OSARGs having τ /P > §
3, is that a large fraction of periodic RR Lyrae and Cepheidvariables apparently favour a finite coherence time. ForCepheids, more than 30 per cent and 90 per cent of firstovertone and fundamental Cepheids, respectively, stronglyprefer the QPO model even though they are very periodicoscillators. RR Lyrae show a similar trend, with 20 per centand 85 per cent of RRc and RRab, respectively, preferring aQPO. This is not due to the presence of Blazhko effect RRLyrae alone, although such Blazhko RR Lyrae do tend toprefer the QPO model (see below). For those objects thatdo prefer the QPO model, 95 per cent of RRc have coher-ence ratios τ /P > τ /P < τ /P <
10, while more than 95 per cent of firstovertone Cepheids have τ /P >
10. As discussed in §
3, thesebehaviors are a consequence of light curve shape. First over-tone Cepheids, RRc, and RRe variables all have relativelysinusoidal light curves, while fundamental mode Cepheids,RRab, and RRd variables have more triangular/sawtoothwaveforms. Qualitatively, the integral fraction curves in Fig-ure 9 favour periodic models more strongly when the lightcurve is sinusoidal and favour QPO models more stronglywhen the light curves are sawtooth-like.We introduced the higher order p = 4 QPO models withperiods locked in a 2:1 ratio and a common coherence timeto explore if this difference was driven by the assumptionof sinusoidal waveforms in the one-component QPO model.We show in Figure 10 that using the two-component QPOmodels halves the number of RRab and RRd variables thatprefer a QPO compared to the one-component model. Theincrease in periodic preference going from a one-componentto the two-component QPO models for the more sinusoidalRRc and RRe sources is marginal. An even larger shift isseen for the fundamental mode Cepheids. These results arein agreement with the behavior we found in § MNRAS , 1–19 (2016) xploring a model for quasi-periodic light curves −15 −10 −5 0 5 10 15 sinh −1 (cid:0) AIC
QPO − AIC
Periodic (cid:1) i n t e g r a l f r a c t i n Q P O P e r i d i c Q P O P e r i d i c Q P O P e r i d i c Q P O P e r i d i c Q P O P e r i d i c Q P O P e r i d i c Q P O P e r i d i c MirasQSOsSRVsOSARGsF Cepheids1 CepheidsRRab+RRdRRc+RRe
Figure 9.
Distribution of selected variables in the AIC likelihood ratio sinh − x with x = ( AIC
QPO − AIC P ) / x >
0, right) or quasi-periodic ( x <
0, left) variability. periments comparing artificial light curves with sawtooth orsinusoid waveforms. It is not surprising, therefore, that ahigher-order QPO model results in a greater preference forperiodicity. The two-component QPO is a better approxi-mation to the sawtooth-like light curves of RRab, RRd, andfundamental mode Cepheids, so the coherence time can actmore like the desired physical parameter, and less as a proxyfor light curve shape.In addition to light curve shape, we also noted in § χ /dof = 1, then the fraction of Cepheidsand RR Lyrae favouring the periodic model increase signifi-cantly, from 70 per cent and 10 per cent to 90 per cent and65 per cent in the case of first overtone and fundamentalCepheids, respectively, and increasing from 80 per cent and15 per cent among RRc/RRe and RRab/RRd to 95 per centand 90 per cent, respectively. The lesson from § τ parameter cannotsimply be interpreted as a measure of coherence because it issignificantly affected by light curve shape. This is, however,no barrier to using the QPO parameters to classify variables.Figure 11 shows the distribution of all the sources in the P − τ , P − σ , and τ − σ QPO parameter spaces. We see that the P − τ plane maps physically distinct variable types to nearlyunique spaces. Below we briefly discuss the properties of eachclass of variable. We do not examine the properties of thevariables in DRW parameter space since such an explorationwas carried out in Koz lowski et al. (2010).If we first consider the distribution of the quasars, we seethat they largely have long periods ( P > ∼
100 days) and lowcoherence ( τ /P < § § MNRAS , 1–19 (2016) Zinn et al. −15 −10 −5 0 5 10 15 sinh −1 (cid:0) AIC
QPO − AIC
Periodic (cid:1) i n t e g r a l f r a c t i o n Q P O P e r i o d i c Q P O P e r i o d i c Q P O P e r i o d i c Q P O P e r i o d i c RRab+RRdRRc+RReF Cepheids1 Cepheids
Figure 10.
Distribution of the RRab+RRd, RRc+RRe, fundamental Cepheids, and first overtone Cepheids in the AIC likelihooddifference for the QPO and periodic models when using either the standard QPO model (dotted) or the higher order model (solid). Thechanges are small for the more sinusoidal RRc+RRe and first overtone Cepheids, but shift strongly toward the periodic case for the moresawtooth-like RRab+RRd and fundamental mode Cepheids.
The various LPVs generally occupy somewhat differentregions of the QPO parameter space. That they have differ-ences in period and amplitude is, of course, already known,so we focus on the the new parameter τ . Figure 12 showsthe distribution of the Miras and SRVs in P - τ space furtherdivided by their C/O abundance class. While not perfectlyseparated by τ at fixed period, it is generally true that theMiras are more coherent (larger τ /P ) than the SRVs andthat C-rich SRVs are more coherent than O-rich SRVs. Sim-ilarly, at fixed amplitude σ , the coherence time separates theO- and C-rich SRVs well, and provides some discriminationbetween the O- and C-rich Miras. It is probably not possibleto well-separate the OSARGs from the other LPVs simplybased on their QPO parameters, possibly because there areall simply part of an evolutionary continuum as suggested bySoszy´nski et al. (2013). While our discussion in § § τ simply as a coher-ence time is risky, it is broadly true that longer period Mirasand SRVs have larger τ /P suggesting that they are more co-herent oscillators. This is also seen simply from folded lightcurves, where longer period Miras look more like the upper panels of Figure 1 and shorter period ones look more likethe lower panels.Although they are highly regular variables, theCepheids are not found to have large coherence times τ inthe QPO models for all the reasons discussed earlier. That τ appears primarily to be a measure of light curve shape ratherthan coherence of oscillation, does lead to a clear separationof the fundamental mode and overtone Cepheids at fixedperiod. Similarly, the RRc and RRe variables tend to havelarger τ than the RRab and RRd variables, but these classesare so well separated in period that the QPO parameter addslittle new information. Blazhko effect RR Lyra should ap-pear as low-coherence QPOs, although they are too rare toexplain the overall distributions. If we look at the parame-ters of the 52 confirmed Blazhko RR Lyrae from Chen et al.(2013) we see that they are located preferentially at lowercoherence ratios among RRab, but are nearly uniform in dis-tribution for RRc and RRe. With a higher order model thatwould make all the normal RR Lyrae coherent oscillators,the coherence time would likely be an efficient means of iden-tifying Blazhko variables. In particular, Figure 13 shows the MNRAS , 1–19 (2016) xploring a model for quasi-periodic light curves Figure 11.
Variable distributions in projections of the P − σ − τ QPO parameters. Note that chemical distinctions are included forSRVs and Miras. In the P − τ plane, dotted lines of constant coherence, τ/P , are shown to guide the eye, where τ/P = 1 is the boundarybetween coherent and incoherent periodicity. The first two columns show τ on the y-axis, whereas the third column is against P . Greystars indicate RR Lyrae classified as Blazhko by Chen et al. (2013). P − τ plane for the RR Lyrae and Cepheid variables usingthe p = 4 QPO model. Both RR Lyrae and Cepheids nowhave significantly higher coherence times, as expected fromthe discussion in §
3, and the Blazhko RR Lyrae are betterseparated.
We have carried out a large scale survey of the propertiesof quasars and variable stars using two different models ofstochastic variability. The models are the two lowest ordermodels within the CARMA family of solutions to stochasticdifferential equations: (1) the damped random walk (DRW)model characterized by an amplitude σ and a damping time τ ; and (2) the quasi-periodic oscillation (QPO) model char-acterized by an amplitude σ , a period P and a coherencetime τ . In general, variable stars are better described bythe QPO model while quasars are better described by theDRW model. This is not a claim that the DRW is a perfectmodel of QSO variability, although it is completely adequatefor our present purposes. There are some ambiguities in thedetails of the model classification because the QPO modelbecomes indistinguishable from the DRW model in the lim-its where the oscillations are very incoherent ( τ /P < τ , becomes shorter than the observing cadence.By combining the results from the two stochastic models, we MNRAS , 1–19 (2016) Zinn et al. P [day] τ [ d a y ] O-Rich SRVsC-Rich SRVsC-Rich MirasO-Rich Miras τ / P = τ / P = τ / P = τ / P = τ / P = Figure 12.
A subset of the top left panel of Figure 11 showing only Miras and Semi-Regular Variables. -2 -1 P [day] -2 -1 τ [ d a y ] F Cepheids1 CepheidsRRab+RRdRRc+RRe τ / P = τ / P = τ / P = τ / P = τ / P = Figure 13.
The P − τ plane for RR Lyrae and Cepheids using the higher-order p = 4 QPO model. Grey stars indicate RR Lyraeclassified as Blazhko by Chen et al. (2013). can significantly increase the efficiency in separating quasarsand variable stars over using the DRW model results alone.In principal, a quantity like the coherence time in theQPO models has the potential to provide important physi-cal information on the nature of pulsations or, if applied tospotted stars, to the evolution time scale of the spots. Wefound, however, that the estimated value of τ in the QPOmodel is also quite sensitive to the deviations of the light curve shape from a sinusoid. Qualitatively, this is most eas-ily understood in terms of the power spectrum of the vari-ability. The QPO model has a power spectrum which is aLorentzian centred on the principle frequency with a widthset by the damping time scale. For a variable with significantdeviations from a sinusoid, such as the “sawtooth–like” longperiod Cepheids or RRab variables, the best p = 2 QPOmodels have a finite coherence time in order to broaden the MNRAS , 1–19 (2016) xploring a model for quasi-periodic light curves Lorentzian to capture some of the variability power in theovertones of the fundamental period. If we use a higher or-der QPO model with terms at both the fundamental periodand the first overtone, in order to better mimic the shapeof such light curves, we find a longer coherence time. Thus,for the coherence time scale of these models to have thatphysical meaning, the model must be of high enough or-der to be able to match the shape of the underlying lightcurve. This is related to the models of a few individual vari-able sources by Kelly et al. (2014), where they continued in-creasing the order of the stochastic model until no increasein the order was statistically justified. For the RRab theyconsidered, they found that a higher-order CARMA(7 , O ( N ) operations usingthe algorithm discovered by Ambikasaran et al. (2014) evenwhen evaluating the full model likelihood rather than simplyusing the forecasting approach of Kelly et al. (2014).Finally, the combination of the DRW and QPO mod-els provides a systematic approach to evaluating period-icity in quasars. There is considerable interest in iden-tifying periodic behavior in quasar light curves as aprobe of quasar binaries and there are now a number ofclaimed detections (e.g., Sillanpaa et al. 1988; Maness et al.2004; Rodriguez et al. 2006; Eracleous et al. 2012; Liu et al.2015; Yan et al. 2015; Graham et al. 2015; Li et al. 2016;Zheng et al. 2016; Bhatta et al. 2016). The problem is thatstandard means of evaluating periodicity, such as a Lomb–Scargle periodogram (Lomb 1976; Scargle 1982), are essen-tially comparing the probability of modeling the source asa sinusoid to the probability of the same light curve beinggenerated by a known level of white noise characterized bythe photometric errors. As discussed by Press (1978) andmore recently by Vaughan et al. (2016), this is not a goodstatistical test when the underlying source class is known tobe variable and with a relatively red noise-like power spec-trum showing more variability on longer time scales. We seethis in our analyses of the MQS quasars, where effects likenoise levels or the structure of the variability can make theQPO model of a quasar more likely than the DRW model.The QPO model we consider here, which spans variabil-ity from truly periodic back to the completely non-periodicDRW model provides a good statistical basis for evaluat-ing periodicity in quasar light curves (modulo the caveats above) because it has a continuous parameter τ to evaluatethe significance of the periodicity. Our approach also has theability to systematically remove light curve means or othercalibration issues for light curves combining multiple sourcesof data. ACKNOWLEDGMENTS
CSK is supported by NSF grant AST-1515927. We wouldlike to thank Profs. M. Kubiak and G. Pietrzy´nski, for-mer members of the OGLE team, for their contributionto the collection of the OGLE photometric data over thepast years. The OGLE project has received funding fromthe Polish National Science Centre grant MAESTRO no.2014/14/A/ST9/00121 to AU. SK acknowledges the finan-cial support from the Polish National Science Center grantno. 2014/15/B/ST9/00093.
REFERENCES
Aigrain S., Parviainen H., Pope B., 2016, preprint,( arXiv:1603.09167 )Ambikasaran S., Foreman-Mackey D., Greengard L., Hogg D. W.,O’Neil M., 2014, preprint, ( arXiv:1403.6015 )Andrae R., Kim D.-W., Bailer-Jones C. A. L., 2013, A&A,554, A137Auvergne M., et al., 2009, A&A, 506, 411Bakos G. ´A., L´az´ar J., Papp I., S´ari P., Green E. M., 2002, PASP,114, 974Balaji B., Croll B., Levine A. M., Rappaport S., 2015, MNRAS,448, 429Ball N. M., Brunner R. J., 2010,International Journal of Modern Physics D, 19, 1049Bedding T. R., Zijlstra A. A., 1998, ApJ, 506, L47Bedding T. R., Kiss L. L., Kjeldsen H., Brewer B. J., Dind Z. E.,Kawaler S. D., Zijlstra A. A., 2005, MNRAS, 361, 1375Bhatta G., et al., 2016, ApJ, 832, 47Bloom J. S., et al., 2012, PASP, 124, 1175Borucki W. J., et al., 2010, Science, 327, 977Brockwell P. J., Davis R. A., 1996, Introduction to Time Seriesand Forecasting. Springer, New YorkBrockwell P. J., Davis R. A., 2002, Introduction toTime Series and Forecasting, 2nd edn. Springer,
Cano Z., Wang S.-Q., Dai Z.-G., Wu X.-F., 2016, preprint,( arXiv:1604.03549 )Chen B.-Q., Jiang B.-W., Yang M., 2013,Research in Astronomy and Astrophysics, 13, 290Christensen-Dalsgaard J., Kjeldsen H., Mattei J. A., 2001, TheAstrophysical Journal Letters, 562, L141Christy R. F., 1966, ApJ, 144, 108Debosscher J., Sarro L. M., Aerts C., Cuypers J., VandenbusscheB., Garrido R., Solano E., 2007, A&A, 475, 1159Dexter J., Agol E., 2011, ApJ, 727, L24Drake A. J., et al., 2009, ApJ, 696, 870Eracleous M., Boroson T. A., Halpern J. P., Liu J., 2012, ApJS,201, 23Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013,PASP, 125, 306Frasca A., Fr¨ohlich H.-E., Bonanno A., Catanzaro G., Biazzo K.,Molenda- ˙Zakowicz J., 2011, A&A, 532, A81Garc´ıa R. A., et al., 2014, A&A, 572, A34Graczyk D., et al., 2011, Acta Astron., 61, 103MNRAS , 1–19 (2016) Zinn et al.
Graham M. J., Djorgovski S. G., Drake A. J., MahabalA. A., Chang M., Stern D., Donalek C., Glikman E., 2014,Mon. Not. Roy. Astron. Soc., 439, 703Graham M. J., et al., 2015, MNRAS, 453, 1562Hamuy M., Phillips M. M., Suntzeff N. B., Schommer R. A., MazaJ., Smith R. C., Lira P., Aviles R., 1996, AJ, 112, 2438He S., Yuan W., Huang J. Z., Long J., Macri L. M., 2016, preprint,( arXiv:1609.06680 )Howarth J. J., 1991, Journal of the British Astronomical Associ-ation, 101, 101Ivezic Z., Tyson J. A., Abel B., Acosta E., Allsman R., AlSayyadY., et al. 2008, preprint, ( arXiv:0805.2366 )Kaiser N., et al., 2002, in Tyson J. A., Wolff S., eds, Society ofPhoto-Optical Instrumentation Engineers (SPIE) ConferenceSeries Vol. 4836, Survey and Other Telescope Technologiesand Discoveries. pp 154–164, doi:10.1117/12.457365Kasliwal V. P., Vogeley M. S., Richards G. T., 2015, MNRAS,451, 4328Kato T., Uemura M., Ishioka R., Nogami D., Kunjaya C., BabaH., Yamaoka H., 2004, PASJ, 56, S1Kelly B. C., Bechtold J., Siemiginowska A., 2009, ApJ, 698, 895Kelly B. C., Becker A. C., Sobolewska M., Siemiginowska A.,Uttley P., 2014, ApJ, 788, 33Kholopov P. N., Samus N. N., Kazarovets E. V., Perova N. B.,1985, Information Bulletin on Variable Stars, 2681Kim D.-W., Protopapas P., Bailer-Jones C. A. L., Byun Y.-I., Chang S.-W., Marquette J.-B., Shin M.-S., 2014, A&A,566, A43Klebesadel R. W., Strong I. B., Olson R. A., 1973, ApJ, 182, L85Koen C., 2005, MNRAS, 361, 887Koen C., 2012, MNRAS, 419, 1197Koz lowski S., 2016, preprint, ( arXiv:1611.08248 )Koz lowski S., et al., 2010, ApJ, 708, 927Koz lowski S., et al., 2013, ApJ, 775, 92LSST Science Collaboration et al., 2009, preprint,( arXiv:0912.0201 )Law N. M., et al., 2009, PASP, 121, 1395Li Y.-R., et al., 2016, ApJ, 822, 4Liu T., et al., 2015, ApJ, 803, L16Lomb N. R., 1976, Ap&SS, 39, 447MacLeod C. L., et al., 2010, ApJ, 721, 1014MacLeod C. L., et al., 2012, ApJ, 753, 106Mackenzie C., Pichara K., Protopapas P., 2016, ApJ, 820, 138Maness H. L., Taylor G. B., Zavala R. T., Peck A. B., PollackL. K., 2004, ApJ, 602, 123McQuillan A., Aigrain S., Mazeh T., 2013, MNRAS, 432, 1203Mushotzky R. F., Edelson R., Baumgartner W., Gandhi P., 2011,The Astrophysical Journal Letters, 743, L12Pejcha O., Kochanek C. S., 2012, ApJ, 748, 107Pepper J., Gould A., Depoy D. L., 2004, in Holt S. S., Dem-ing D., eds, American Institute of Physics Conference Se-ries Vol. 713, The Search for Other Worlds. pp 185–188( arXiv:astro-ph/0401220 ), doi:10.1063/1.1774522Pojma´nski G., 1997, Acta Astron., 47, 467Pojma´nski G., 2002, Acta Astron., 52, 397Pollacco D. L., et al., 2006, PASP, 118, 1407Press W. H., 1978, Comments on Astrophysics, 7, 103Pskovskii I. P., 1977, Soviet Ast., 21, 675Reinhold T., Reiners A., Basri G., 2013, A&A, 560, A4Richards J. W., Starr D. L., Miller A. A., Bloom J. S., ButlerN. R., Brink H., Crellin-Quick A., 2012, The AstrophysicalJournal Supplement Series, 203, 32Ricker G. R., et al., 2014, in Society of Photo-Optical In-strumentation Engineers (SPIE) Conference Series. p. 20( arXiv:1406.0151 ), doi:10.1117/12.2063489Robinson E. L., 1976, ARA&A, 14, 119Rodriguez C., Taylor G. B., Zavala R. T., Peck A. B., PollackL. K., Romani R. W., 2006, ApJ, 646, 49 Rybicki G. B., Press W. H., 1992, ApJ, 398, 169Sako M., et al., 2008, AJ, 135, 348Sako M., et al., 2011, ApJ, 738, 162Sarro L. M., Debosscher J., L´opez M., Aerts C., 2009, A&A,494, 739Scargle J. D., 1982, ApJ, 263, 835Shappee B., et al., 2014, in American Astronomical Society Meet-ing Abstracts arXiv:1310.6774 ) The mean of the best-fitting light curves, ˆ s , approximatingthe observed light curve, ~y , may be calculated using theexpression derived in Rybicki & Press (1992),ˆ s = SC − ( ~y − L~q ) , (1)with variance from the true light curve given by h ( ~s − ˆ s ) i = S − S T C ⊥ S, (2)where C ⊥ ≡ (cid:16) C − − C − L ( L T C − L ) − L T C − (cid:17) − . (3)As previously discussed in Rybicki & Press (1992) andKoz lowski et al. (2010), the CARMA formalism can be gen-eralized to take into account arbitrary trends. This assumes MNRAS , 1–19 (2016) xploring a model for quasi-periodic light curves that the data are modelled as ~y = ~s + ~n + L~q, (4)where ~s is the expected signal, described in this case by aGaussian process, ~n is a Gaussian noise term, and L~q de-scribes a trend. L is a N × ℓ matrix, with ℓ representing thenumber of temporal trends to include in the fit and N isthe number of data points. The column vector ~q is of length ℓ , and contains the polynomial coefficients for the trend.In this notation, we can include a mean and linear trendin the model, q + ( t − t ) q , by setting q = ( q , q ) and L i = 1 , L i = t i − t .As another example, to remove jumps within a lightcurve, like the calibration shifts between the OGLE-III andIV QSO light curves, one would again have q = ( q , q ), butdefine L to be L = (1 ,
0) for the M OGLE-III points of thelight curve and L = (0 ,
1) for the N − M OGLE-IV pointsof the light curve.
This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000