Are the Variability Properties of the Kepler AGN Light Curves Consistent with a Damped Random Walk?
aa r X i v : . [ a s t r o - ph . GA ] J un Mon. Not. R. Astron. Soc. , 1–20 (2015) Printed 3 June 2015 (MN L A TEX style file v2.2)
Are the Variability Properties of the Kepler AGN LightCurves Consistent with a Damped Random Walk?
Vishal P. Kasliwal ⋆ , Michael S. Vogeley ⋆ , & Gordon T. Richards ⋆ Department of Physics, Drexel University, 3141 Chestnut St., Philadelphia, PA 19104, USA
Accepted 2015 June 1. Received 2015 April 25; in original form 2014 August 1
ABSTRACT
We test the consistency of active galactic nuclei (AGN) optical flux variability with the damped random walk (DRW) model. Our sample consists of 20 multi-quarter
Kepler
AGN light curves including both Type 1 and 2 Seyferts, radio-loud and -quiet AGN,quasars, and blazars.
Kepler observations of AGN light curves offer a unique insightinto the variability properties of AGN light curves because of the very rapid (11 . − . over a period of 3.5 yr. We categorize the light curves of all 20 objectsbased on visual similarities and find that the light curves fall into 5 broad categories.We measure the first order structure function of these light curves and model theobserved light curve with a general broken power-law PSD characterized by a short-timescale power-law index γ and turnover timescale τ . We find that less than half theobjects are consistent with a DRW and observe variability on short timescales ( ∼ τ ranges from ∼ −
135 d. Interesting structure functionfeatures include pronounced dips on rest-frame timescales ranging from 10 −
100 dand varying slopes on different timescales. The range of observed short-timescale PSDslopes and the presence of dip and varying slope features suggests that the DRW modelmay not be appropriate for all AGN. We conclude that AGN variability is a complexphenomenon that requires a more sophisticated statistical treatment.
Key words: galaxies: active – galaxies: Seyfert – quasars: general – BL Lacertaeobjects: general – accretion, accretion discs
The continuum of the X-ray through optical radiation spec-trum observed in active galactic nuclei (AGN) arises in theaccretion disk surrounding the central supermassive blackhole (Shakura 1973; Shakura & Sunyaev 1973). It is wellknown that this region of the spectrum exhibits strong,stochastic variability at the 200 per cent (X-ray) to 10per cent (optical) level (Mushotzky, Done, & Pounds 1993;Wagner & Witzel 1995; Krolik 1999). Variability is observedover a wide range of timescales ranging from hours to monthsto years. The physical mechanisms driving this variabilityare not well understood — models range from X-ray flaresthat drive optical variability (Goosmann et al. 2006) to lo-cal variations in the plasma viscosity of the accretion disk(Lyubarskii 1997) caused by small scale angular momentumoutflows triggered by local dynamo processes (King et al.2004).The accretion disk of a typical 10 M ⊙ supermassive ⋆ E-mail: [email protected] (VPK); [email protected] (MSV);[email protected] (GTR) black hole is too small to be imaged. Models of AGNvariability can only be tested by comparing the stochas-tic properties of observed AGN light curves to the behav-ior predicted by the variability model. The most popularAGN variability model is the damped random walk (DRW)model or 1-parameter Auto-Regressive, AR(1), process(Kelly, Bechtold, & Siemiginowska 2009; Kozlowski et al.2010; MacLeod et al. 2010; Zu et al. 2013). This model isextremely simple and does a good job of fitting ground-based variability data, yet we will demonstrate that it failsto capture the full range of variability behavior exhibited byAGN. We begin by outlining the terminology and mathe-matics used in discussing variability and the AR(1) processin the next section.
The light curve of an AGN may be regarded as a discretetime series that samples a continuous process. Measurementsof the source flux F i are made at times t i for N instants intime. Ideally, these measurements are obtained at fixed in- c (cid:13) Vishal P. Kasliwal, Michael S. Vogeley, and Gordon T. Richards crements separated by a constant sampling interval δt . Insuch cases, the time interval ∆ t between any two flux mea-surements, F i & F j , is related to the sampling interval by∆ t = nδt, (1)where n = i − j with 1 n < N . This kind of samplingpattern is usually not possible in the case of ground-basedastronomical time series because of interruptions caused byfactors such as the weather and the availability of the targetin the night sky. Even in the case of space-based measure-ments of the series, it is possible to have ‘missing values’,i.e. values of t i with i = nδt , for which there is no corre-sponding measurement of the flux F i . The presence of thesemissing values complicates the analysis process. In the caseof a time series that exhibits stochastic behavior, the goalof time series analysis is to characterize the joint probabilitydistribution of F i for the measured data points. If the jointprobability distribution of the data points is independent oftime, the time series is said to be stationary , i.e. one wouldobserve the same light curve (from a statistical standpoint)at all times.The timescale on which the bulk fueling rates of theAGN vary is in the hundreds of thousands to millions ofyears and therefore cannot be directly probed by obser-vation for any individual AGN. It is reasonable to expectthat the stochastic variability observed in the flux timeseries of individual AGN is caused by frequently reoccur-ring short-lived processes local to the accretion disk, i.e.the time series is stationary on the timescales we probed.It is well known from the theory of time series analysisthat a stationary time series can be modeled by an Auto-Regressive Moving Average, or ARMA, process (Hamilton1994; Woodward, Gray, & Elliot 2012; Prado & West 2010;Box, Jenkins, & Reinsel 2006; Brockwell & Davis 2002) ofappropriate order (p,q). The general form of an ARMA(p,q)process, i.e. a process with p autoregressive terms and qmoving average terms, is given by F i = p X m =1 φ m F i − m + q X n =1 θ i − n w i − n + w i , (2)where φ m are autoregression coefficients, θ n are moving av-erage coefficients, and w i are known as ‘innovations’ in sta-tistical literature. Autoregressive terms introduce infinite-duration correlations in the data while moving average termsintroduce finite-duration correlations. It is possible to repre-sent any stationary process as a pure autoregressive or mov-ing average process, albeit perhaps of infinite order (Scargle1981). Innovations drive the stochastic behavior of the pro-cess. Each innovation is an independent random number typ-ically drawn from a Gaussian random distribution with zeromean and variance σ ∝ δt . To make the process station-ary we must restrict | φ i − m | < m . It isclear from previous work (Kelly, Bechtold, & Siemiginowska2009; Kozlowski et al. 2010; MacLeod et al. 2010) that AGNlight curves exhibit long-term correlations. The simplestARMA process that exhibits long-term correlations is thepurely autoregressive (no moving average terms) AR(1) pro-cess or damped random walk (DRW) (Ivezi´c et al. 2014)given by F i = φ F i − + w i , (3)with 0 < φ <
1. The lone autoregression coefficient φ con- nects the current value of the series to the previous value,making the AR(1) process a memoryless or Markov process.The numerical value of φ depends upon the sampling in-terval δt . For a given sampling interval, a value of φ ∼ F i ≃ F i − and so the value of φ sets a char-acteristic timescale τ after which the contribution of theprevious value of the series will be relatively insignificant,i.e. τ is the decorrelation timescale . Numerically, φ = e − δtτ , (4)where δt is the sampling interval and τ is the aforementionedcharacteristic decorrelation timescale (Gillespie 1996). If δt ≫ τ , then φ ∼ δt ≪ τ , then φ ∼ φ <
1, i.e. δt < τ . It turns out that the vari-ance of the innovations is also linked to the decorrelationtimescale: h w i i = σ (1 − e − δtτ ) , (5)where σ allows the amplitude of the variance of the innova-tions to vary independently of δt and τ . For a fixed value of σ , if the sampling interval is much longer than τ , the vari-ance of the innovations will be very close to σ . Combinedwith the accompanying small value of φ , the resulting timeseries will look close to pure white noise. If, on the otherhand, we choose a sampling interval much smaller than thedecorrelation timescale, the variance of the innovations isforced to be very close to zero and the value of φ is veryclose to 1, resulting in very obvious correlations between thedata points.Physically, the AR(1) process could result from an ac-cretion disk with local ‘spots’ that contribute more or lessflux than the mean flux level of the disk. These spots ap-pear at random and dissipate over some characteristic phys-ical timescale (Agol & Dexter 2011). Under the assumptionthat the distribution of the contribution to the total fluxfrom each spot is Gaussian, the Central Limit Theorem as-sures us that the overall changes in flux, i.e. the innovationsin Eq. (2) and in Eq. (3), are drawn from a Gaussian dis-tribution. Long term correlations exist because the spotsdo not dissipate instantaneously. The power spectral density (PSD) of this time series follows a bent power-law given by P SD AR (1) ( f ) = 4 σ τ πfτ ) , (6)while the auto-covariance function (ACVF), related to thePSD via a Fourier Transform (Wiener-Khintchine theorem),is given by ACV F AR (1) (∆ t ) = σ e − ∆ tτ , (7)where σ is the low-frequency asymptotic amplitude and τ isthe characteristic timescale at which the PSD turns over anddecays (Gillespie 1996). Hence the auto-correlation function (ACF) is given by ACF AR (1) (∆ t ) = e − ∆ tτ , (8)Does the AR(1) process actually describe the variabil-ity of AGN? Recent ground based studies (Zu et al. 2013;Graham et al. 2014) indicate that on very short timescales, c (cid:13) , 1–20 re the Variability Properties of the Kepler AGN Light Curves Consistent with a Damped Random Walk? AGN light curves show strong autocorrelation with PSDslopes steeper than 1 /f . To model stochastic light curveswith arbitrary PSD slopes, we introduce the damped power-law (DPL) model by generalizing the AR(1) PSD in Eq. (6)to P SD
DPL ( f ) = σ Eff πfτ ) γ , (9)where we have lumped the quantities in the numerator ofEq. (6) into a single variable σ Eff and introduced γ to al-low the logarithmic slope of the PSD to be a free parameter.When γ <
2, the process exhibits weaker autocorrelation onshort timescales than the AR(1) i.e. the time series looksless smooth. When γ >
2, the process shows stronger au-tocorrelation on short timescales than the AR(1) i.e. thetime series looks smoother. This DPL model is a simplifi-cation of the popular ‘Bending Power-Law’ (BPL) model ofMcHardy et al. (2004):
P SD
BPL ( f ) = Af − α L Y i =1
11 + ( ff Bi ) α i +1 − α i , (10)where α i are the logarithmic PSD slopes above the corre-sponding ‘bend’ frequencies f Bi . Allowing the logarithmicPSD slope to be a free parameter is similar to the structurefunction parametrization chosen by Schmidt et al. (2010): SF (∆ t ) = A (cid:18) ∆ t yr (cid:19) λ , (11)where λ plays the same role as γ in the DPL model–theyallow the correlation level at fixed time-lag to vary from thatpredicted by a pure AR(1) process. To test our DPL model,we use light curve data from the NASA Kepler mission.
KEPLER
MISSION
The
Kepler mission was designed to survey the Miky Wayfor Earth-sized and smaller planets around the habitablezone by using the transit-method to detect exo-planets.
Kepler is designed as a Schmidt-camera with a primaryaperture of 0.95-m, a fixed 115.6 deg FOV, and a half-maximum bandpass from 435-nm to 845-nm with > Ke-pler focal-plane consists of 42 CCDs arranged in 21 mod-ules that are readout by 4 channels each. Each 50 × × Kepler
PSF is intentionally large–6.4pixels near the center of the FOV–in order to increase thephotometric precision (Gilliland 2004) by ensuring that nopixel contains more than 60 per cent of the flux from apoint source target. This has two consequences that com-plicate simple aperture photometry: 1. point sources, suchas stars, appear extended, and 2. faint background sources‘crowd’ the aperture used to measure the flux from a tar-get of interest. The signal at each pixel on the CCDs ismeasured over an interval known as the integration time or frame which consists of a variable exposure time and a fixed readout time of 0.520 s. The flight exposure time is 6.02 syielding a flight integration time of 6.54 s per frame. Twodatasets can be created by Kepler : the long cadence (LC)data and the short cadence (SC) data . The long cadence data is integrated over 270 frames (29.430 min), while theshort cadence data is integrated over 9 frames (0.981 min).Most of the objects observed by us have no SC data–only LCdata was obtained. Furthermore, the SC data set is harderto calibrate because of the small number of SC quiet tar-gets observed (Kinemuchi et al. 2012). For these reasons weonly use LC data in our analysis. Kinemuchi et al. (2012)provides an excellent overview of the final data productssupplied by
Kepler , including strategies for dealing with thevarious types of artifacts present in the data.The LC data-set poses several calibration chal-lenges. Certain readout channels on specific CCD mod-ules have known performance issues discussed in theKepler Instrument Handbook (KSCI-19033) including out-of-spec read noise levels and gain. More troubling are the‘moire’ and ‘rolling band’ effects seen in some channels(Kolodziejczak et al. 2010). The combination of module andchannel is usually quoted in the format module
Kepler target will fall on a fixed sequence of 4 module andchannel combinations over the course of a year, i.e. since thespacecraft rolls 90 o every quarter, at the beginning of every5 th quarter targets will fall on the same module Kepler data is available in both uncorrected (SAP flux) and cor-rected (PDCSAP flux) forms. Corrections to the calibrateddata are obtained using the PDC module (Stumpe et al.2012; Smith et al. 2012).The PDC module identifies features common to hun-dreds of quiet targets on each CCD by examining the cross-correlation between targets. The most correlated targets areused to create 16 best-fit vectors called Cotrending BasisVectors (CBVs) using Singular Value Decomposition (SVD).Corrections to the calibrated data are obtained by remov-ing the CBVs from the data using a weighted normaliza-tion. The algorithm used by PDC to compute these CBVweights changed with the introduction of Ver 8.0 of the
Ke-pler
Science Operations Center (SOC) pipeline (Data Re-lease 14). Prior to pipeline Ver. 8.0, the PDC module useda Least Squares (PDC-LS) approach to compute the CBVweights. Ver. 8.0 introduced a Bayesian
Maximum A Poste-riori (PDC-MAP) algorithm compute the CBV weights andoptimally remove spacecraft-induced trends while minimiz-ing the removal of the true underlying signal (Stumpe et al.2012). The Bayesian-MAP procedure assumes that the theobserved signal may be represented as F = Hθ + ǫ, (12)where H is a design matrix consisting of the CBV vectors(typically 4), θ is the vector of the CBV weights (which c (cid:13) , 1–20 Vishal P. Kasliwal, Michael S. Vogeley, and Gordon T. Richards is what the PDC module is trying to determine), and ǫ ∼ (0 , Σ N ) is the (uncorrelated, i.i.d.) instrument noise vector.The log-likelihood of an observed light curve, assuming thismodel for the data, is then given byln p ( F | θ ) = − N ln(2 π ) + ln( | Σ N | ) + ( F − Hθ ) T Σ − N ( F − Hθ )2 . (13)Maximization of this likelihood can be performed analyti-cally, however the resulting de-trended light curve will nec-essarily have intrinsic variability removed because LS cal-culates the CBV weights by projecting the raw light curveonto each CBV—any coincidences between the true lightcurve and the CBV will spuriously add to the CBV weight.To avoid this behavior, the new PDC-MAP approachapplies a weighted Bayesian prior, p ( θ ), to the likelihoodin (13) to bias the best-fit CBV weights toward removingonly the spacecraft-induced variability signal while leavingintrinsic variability intact. p ( θ ), is constructed by using LSto de-trend quiet targets in the phase-space vicinity (RA,Dec, and apparent magnitude) of the object of interest. Theprior weight, W pr , is calculated such that it is closer to 0 forrelatively quieter targets and closer to 1 for relatively activetargets. The PDC algorithm computes b θ = arg max θ (ln p ( F | θ ) + W pr ln p ( θ )) . (14)Therefore p ( θ ) biases the CBV weights towards values cor-responding to those for neighboring quiet targets. Since theCBV weights are biased towards removing just enough fea-tures in the light curve to render featureless light curves forneighboring quiet targets, intrinsic variability in the lightcurve of a truly variable object should be relatively unaf-fected by the de-trending procedure. Further details on thePDC de-trending algorithm may be found in Smith et al.(2012).The benefit of using the Kepler -MAST supplied PD-CSAP flux light curves is that for most targets, the PDC-MAP algorithm produces optimally de-trended light curvesfree of instrument effects. At the same time, the drawbackof doing so is that it is theoretically possible for the PDC-MAP pipeline to remove intrinsic variability, particularlywhen the target exhibits intrinsic long-term trends that can-cel out the systematic trend, or when the location of thetarget in position-luminosity phase-space is lacking an ade-quate number of close neighbors. On the other hand, usingthe raw SAP flux light curves and performing a customizedde-trending using the PyKE tool kepcotrend implies thatthe CBV weights must be computed using no informationfrom neighboring targets and hence are likely to project outintrinsic variability in the light curve.We use
Kepler -MAST supplied light-curves from DataRelease 23 that have been processed using version 9.1 of the
Kepler
SOC pipeline (PDC module uses PDC-MAP algo-rithm) and only perform a minor correction to remove thelarge quarterly offsets in the data caused by spacecraft rolls.
Only ∼ Kepler
FOV prior tothe start of the mission due to the lack of coverage from ex-isting extragalactic surveys. Since the beginning of the mis-sion in May 2009, Guest Observer (GO) efforts to find more AGN in the
Kepler
FOV by multiple groups (Carini & Ryle2012; Edelson & Malkan 2012; Wehrle et al. 2013) have ledto a total of roughly 80 AGN of various types being mon-itored by
Kepler . The May 2013 failure of a critical reac-tion wheel led to the termination of the original
Kepler mis-sion. We wish to examine variability properties across a widerange of AGN type and redshift; however a large number ofthe
Kepler
AGN were selected photometrically and do nothave redshifts. As such, we restrict our analysis to the 20that are spectroscopically confirmed AGN. This allows us toreject non-AGN contaminants, and also allows us to studyvariability behavior as a function of AGN type and performcomparisons in the rest-frame of the object.Of the 7 AGN known to lie in the
Kepler
FOV priorto 2009, the best studied is kplr006932990, also known asZw 229-15. The VCV Catalog (V´eron-Cetty & V´eron 2010)lists kplr006932990 as a radio-quiet Type 1 Seyfert at a red-shift of 0.027 with V-band absolute magnitude M V = − . β reverberation mapping, Barth et al. (2011) ob-tain a virial black hole mass of M BH = 1 . +0 . − . × M ⊙ .The other Type 1 AGN known to exist in the Kepler
FOVprior to the commencement of the mission is the X-Raysource kplr09650715 (RXS J19298+4622), which is a radio-quiet Seyfert 1 at z = 0 .
127 with V-band absolute mag-nitude M V = − . Kepler
FOV:kplr008703536 (IGR J19473+4452), which is a radio-quietSeyfert 2 (Masetti et al. 2006) at z = 0 . M V = − . z = 0 .
513 with photographic magnitude M O = − . Kepler
AGN include theFSRQ blazars kplr11606854 (MG4 J191843+4937) at z =0 .
926 and kplr010663134 (Q 1922+4748) at z = 1 . Kepler
AGN isthe quasar kplr012208602 (4C 50.47) at a redshift of 1 . α = − . c (cid:13) , 1–20 re the Variability Properties of the Kepler AGN Light Curves Consistent with a Damped Random Walk? Table 1.
Kepler
AGN Sample.Identifier Physical ParametersKeplerID Alt. ID RA Dec z M g AGN TypeHrs min sec Deg min sec6932990 Zw229-15 †
19 05 25.969 +42 27 40.07 0.0275 -23.8 Sy1 † ⋆
19 25 02.181 +50 43 13.95 0.067 -21.9 Sy111178007 W2R 1858+48 ⋆
18 58 01.111 +48 50 23.39 0.079 -20.2 Sy19650715 RXS J19298+4622 †
19 29 50.490 +46 22 23.59 0.127 -21.8 † Sy1 † ⋆
19 10 02.496 +38 00 09.47 0.130 -20.4 Sy13347632 W2R 1931+38 ⋆
19 31 15.485 +38 28 17.29 0.158 -21.0 Sy2 ? ⋆
19 15 09.127 +41 02 39.08 0.222 -22.1 Sy12694186 W2R 1904+37 ⋆
19 04 58.674 +37 55 41.09 0.089 -23.8 Sy1 ? ⋆
19 22 11.234 +45 38 06.16 0.115 -22.2 Sy1.9 ⋆ ⋆
18 45 59.577 +48 16 47.57 0.152 -21.2 Sy1 ? ⋆
19 20 47.750 +38 26 41.28 0.368 -23.0 Sy17610713 1931+43 ⋆
19 31 12.566 +43 13 27.62 0.439 -24.7 QSO6690887 W2R 1926+42 ⋆
19 26 31.089 +42 09 59.12 0.154 -21.7 BL-Lac ⋆, § ⋆
19 14 15.492 +42 04 59.88 0.502 -24.6 QSO5597763 W2R 1853+40 ⋆
18 53 19.284 +40 53 36.42 0.625 -25.1 QSO11606854 MG4 J191843+4937 ¶
19 18 45.617 +49 37 55.06 0.926 -25.0 FSRQ ‡ †
19 23 27.234 +47 54 17.00 1.520 -25.2 QSO † ; FSRQ ‡ ; Jet k †
19 47 19.308 +44 49 42.32 0.0539 -21.7 Sy2 † †
19 09 46.501 +48 34 32.26 0.513 -23.2 Sy1.5 † ; FSRQ ‡ †
19 26 06.318 +50 52 57.14 1.098 -24.7 QSO † ; Jet k ⋆ Reliable Identifications of Active Galactic Nuclei from the
WISE , 2MASS, and
ROSAT
All-Sky Surveys (Edelson & Malkan 2012) † A Catalog of Quasars and Active Nuclei: 13th Ed. (V´eron-Cetty & V´eron 2010) ‡ An All-Sky Survey of Flat-Spectrum Radio Sources (Healey et al. 2007) § Kepler Observations of Rapid Optical Variability in the BL Lac Object W2R1926+42 (Edelson et al. 2013) k A New List of Extra-Galactic Radio Jets (Liu & Zhang 2002) ¶ CGRaBS: An All-Sky Survey of Gamma-Ray Blazar Candidates (Healey et al. 2008) ? Weak spectral features–unsure sub-classificationThe AGN are grouped based on the light curve categories discussed in Sec. 6. Col. 6 lists SDSS g-band absolute magnitudes computedfrom the SDSS apparent magnitude supplied in the Kepler Input Catalog (2011) with cosmological parameters H = 67 .
77 kms − Mpc − , Ω M = 0 . Λ = 0 . k-correction using Eq. (10.29) in Peterson (2003) assuming α = − . z = 0. We present tentative sub-classifications for 11 AGN inCol. 7 using imaging from the STScI Digitized Sky Survey, SDSS and spectra in Edelson & Malkan (2012). In the case of broad-lineobjects we distinguish between Seyfert 1s and QSOs based on the presence of a galactic component in STscI DSS/SDSS imaging andalso based on M g > −
24 for Seyferts. and k-corrected absolute magnitude (traditionally M i . −
22 for QSOs from Richards et al. (2006)) computed fromthe Kepler Input Catalog (2011) SDSS g-band apparentmagnitude available on MAST.Table 1 presents positional and physical parametersfor the light curves from Fig. 3. Column 1 lists the
Ke-plerID used when querying the
Kepler
MAST database(Kepler Archive Manual KDMC-10008-005). Column 2 liststhe ‘common’ name of the object along with a pri-mary reference for the object, usually drawn from eitherV´eron-Cetty & V´eron (2010) or Edelson & Malkan (2012).The next two columns supply the J2000.0 RA and Dec lo-cation of the object as listed in the Kepler Input Catalog(2011). The next column supplies the redshift of the ob-ject as listed in the primary reference for the object. Thenext column supplies the k-corrected
SDSS g-band absolutemagnitude ( M g ) based on the apparent SDSS g-band KICmagnitude ( m g ). As described in the Kepler Input Catalog(2011), the galactic extinction corrected m g was determinedduring the Stellar Classification Project (SCP) by usingDAOPHOT (Stetson 1987) to perform PSF photometry onstar-like objects. We have converted these g-band magni-tudes to absolute magnitudes at z = 0 using k-corrections from Peterson (2003) with cosmological parameters H = 67 .
77 kms − Mpc − , Ω M = 0 . Λ = 0 . α = − . Kepler exhibit huge quarterly off-sets and require ‘stitching’ before we may analyze them.The next section discusses how we stitch the light curvestogether.
As discussed in Sec. 3 and 4, quarterly discontinuities ex-ist in
Kepler data. Once every 3 months, the spacecraftperforms a roll maneuver to transmit data to the Earth.After each roll maneuver, targets fall on different CCDs re-sulting in slightly different distributions of the target fluxacross adjacent pixels. To compensate, the target apertureis redefined after every roll maneuver. Due to the large sizeof the
Kepler
PSF relative to the aperture size, differentapertures contain different portions of the flux from the tar- c (cid:13) , 1–20 Vishal P. Kasliwal, Michael S. Vogeley, and Gordon T. Richards get. This results in the severe discontinuities observed inthe target flux. Fig. 1 shows the discontinuities present inthe PDCSAP light curve of kplr006932990. The top panelis plotted in the observed frame where the sampling inter-val is 29.42 min. The bottom panel is plotted in the restframe of kplr006932990 with sampling interval 28.64 minbecause of the cosmological redshift factor of 1 + z = 1 . Kepler samplinginterval for more distant objects resulting in shorter butmore densely sampled light curves. The discontinuities inthe top panel of Fig. 1 must be removed before the longterm variability properties of this object can be studied. Onecould remove such discontinuities by re-calibrating the dataon an object-by-object basis, but this requires considerablemanual effort.
Kepler data products include both the lightcurves of the target as well as the
Target Pixel File (TPF)corresponding to the target. The TPF consist of ‘postagestamp’ snapshots of the object taken at every cadence. Thelight curve of the object is constructed by defining an aper-ture consisting of a subset of the pixels in the postagestamp. This aperture is determined by the
Kepler pipelinebased on the Kepler Input Catalog (2011) and is optimizedfor stellar targets as outlined in Data Processing Handbook(KSCI-19081-001). By re-defining the target aperture byhand, it is possible to reduce the discontinuities in the fluxacross quarters as shown in Carini & Ryle (2012). Such re-extracted light curves must then be de-trended of spacecraftinduced features by removing the CBVs (discussed in 3) withcustom weights determined using the PyKE tool kepcotrend .A simpler but cruder approach to removing quarterlydiscontinuities is to match a suitable metric across quar-ter boundaries (Kinemuchi et al. 2012). Mushotzky et al.(2011) perform a simple end-matching to stitch light curvesacross quarters, i.e. they correct the measured flux values inthe second of every sequential pair of quarters by the differ-ence between the flux value of the last point of the leadingmember, and the flux value of the first point of the trail-ing member of the pair of quarters. This method is quitesuitable for the high signal-to-noise (S/N) objects studiedin Mushotzky et al. (2011), but does not work as well as thevariability signal level drops closer to the instrument noiselevel in the dimmer targets. Edelson et al. (2013) stitch thelight curve of the BL Lac object kplr006690887 across quar-ters by using the per-quarter average flux value to determinea multiplicative factor that they then use to scale quarters.Revalski et al. (2014) use a similar procedure but determinetheir multiplicative scaling factor from the average of thelast 20 points in the quarter preceding the discontinuity andthe first 20 points in the quarter following the discontinuity.We stitch the light curves across quarterly discontinu-ities by matching the average flux value of the last 100 pointsof the leading quarter to the first 100 points of the trailingquarter in every sequential pair of quarters. The method isrobust to outliers and is relatively unaffected by low S/N. Itmay be argued that it is more correct to use a fixed durationin time at the beginning and ends of quarters to computeaverage fluxes, rather than a fixed number of points. Usinga variable number of points determined by a constant timewindow introduces variations in the reliability of the statisticfor removing discontinuities; a one day window in the refer- obs (d) F l u x ( e − ) Quarter Start Quarter Finish
200 400 600 800 1000 1200t rest (d) F l u x ( e − ) Figure 1.
PDCSAP light curve of the Seyfert 1 kplr006932990(Zw 229-15) — Top panel: unstitched; bottom panel: stitched.The huge discontinuities in flux in the top panel are caused byredefinition of the target aperture after the spacecraft performsa roll maneuver as explained in Sec. 5 and are endemic to all the
Kepler light curves. Our stitching technique serves to remove thequarterly discontinuities present in the light curve obtained fromMAST. ence frame of a z = 1 object will have twice as many fluxmeasurements than the same window in the case of a z = 0object. Fixing the number of points has the benefit that itis equally applicable to stitching together the light curvesof objects at unknown redshifts - the stitched light curvesof such objects can stand on equal footing with the stitchedlight curves of objects with known redshift. Fig. 1 shows anapplication of our method to the high S/N light curve ofkplr006932990. The top panel presents the un-stitched lightcurve while the bottom panel presents the stitched lightcurve corrected to the rest-frame of the AGN. As can beseen in the figure, our stitching method removes the dis-continuities present in the original light curve. Similar re-sults may be expected from the end-matching technique ofMushotzky et al. (2011). In contrast to the high S/N lightcurve of kplr006932990, Fig. 2 shows the un-stitched andstitched light curves of the low S/N AGN kplr003337670.The relative thickness of this light curve (caused by scat-ter from measurement noise) makes it essential to considera suitably large number of points at the ends of quarterswhen stitching. Simple end-matching fails to stitch this lightcurve, and using a smaller number of points in the stitchingresults in obvious flux mis-matches. We use our stitchingalgorithm to remove the quarterly discontinuities in all 20AGN light curves and discuss the resulting multi-quarter Kepler light curves in the next section. c (cid:13)000
Kepler light curves. Our stitching technique serves to remove thequarterly discontinuities present in the light curve obtained fromMAST. ence frame of a z = 1 object will have twice as many fluxmeasurements than the same window in the case of a z = 0object. Fixing the number of points has the benefit that itis equally applicable to stitching together the light curvesof objects at unknown redshifts - the stitched light curvesof such objects can stand on equal footing with the stitchedlight curves of objects with known redshift. Fig. 1 shows anapplication of our method to the high S/N light curve ofkplr006932990. The top panel presents the un-stitched lightcurve while the bottom panel presents the stitched lightcurve corrected to the rest-frame of the AGN. As can beseen in the figure, our stitching method removes the dis-continuities present in the original light curve. Similar re-sults may be expected from the end-matching technique ofMushotzky et al. (2011). In contrast to the high S/N lightcurve of kplr006932990, Fig. 2 shows the un-stitched andstitched light curves of the low S/N AGN kplr003337670.The relative thickness of this light curve (caused by scat-ter from measurement noise) makes it essential to considera suitably large number of points at the ends of quarterswhen stitching. Simple end-matching fails to stitch this lightcurve, and using a smaller number of points in the stitchingresults in obvious flux mis-matches. We use our stitchingalgorithm to remove the quarterly discontinuities in all 20AGN light curves and discuss the resulting multi-quarter Kepler light curves in the next section. c (cid:13)000 , 1–20 re the Variability Properties of the Kepler AGN Light Curves Consistent with a Damped Random Walk? obs (d) F l u x ( e − ) Quarter Start Quarter Finish
200 400t rest (d) F l u x ( e − ) Figure 2.
PDCSAP light curve of the Seyfert 1 kplr003337670 —Top panel: unstitched; bottom panel: stitched. The noisiness ofthis light curve makes it essential to use a large (100) number ofpoints at the edges of each quarters when performing the stitch.End-matching is completely ineffective while using a smaller num-ber of points at the quarter edges results in obvious mis-matchesacross the quarter boundaries.
We divide the AGN in the sample into 5 categories of ob-jects based on the visual appearance of their light curve. Fig.3 shows all 20 AGN grouped by visual features present ineach light curve. Within each category, we order the AGNby z . All 20 light curves are plotted over the same rest-frametime range on the x-axis. Category 1 consists of 3 objectsthat have the most stochastic-looking light curves i.e. lightcurves bereft of any sinusoidal features or flares. Weak rip-ple features appear in the light curves of Category 2 objects(4 AGN). Category 3 consists of 5 objects that have pro-nounced ripple features in their light curves, while category4 objects have light curves with flaring behavior (5 AGN).The last category of 3 objects have primarily featureless lightcurves consistent with marginal levels of variability.The light curves show a wide variety of different types ofbehavior. There are indications that individual light curvescan show more than one type of feature. Category 1 ob-jects (kplr6932990 , kplr012158940, and kplr011178007) areshown in the 1 st row of Fig. 3. These objects have the most‘stochastic’-looking light curves and are mostly free of anyrecurring features and trends. They are also the closest AGNin this study ( z < . Table 2.
Kepler
AGN Detector Sequences.Identifier Detector PropertiesKeplerID Module ⋆ , 8.24, 12.40 ‡ , 18.6412158940 24.81 ⋆ , 10.29, 2.1, 16.5311178007 23.79 ⋆ , 15.51 § , 3.7 ¶ , 11 . § † , 12.38, 18.62 ⋆ k , 6.14 ⋆ ‡ , 10.31, 2.35781475 7.17, 17.57, 19.65, 9.252694186 22.75 ⋆ , 20.71, 4.11 ⋆ , 6.159215110 13.41, 13.42, 13.43, 13.44 ‡† ⋆ , 6.15, 22.75 ⋆ ⋆ , 10.29, 2.17610713 8.22 † , 12.38, 18.62 ⋆ , 14.466690887 12.37, 18.61, 14.45, 8.216595745 8.22 † , 12.38, 18.62 ⋆ , 14.465597763 23.80, 15.52 ⋆ , 3.8 ¶ , 11.36 ‡ ⋆ , 2.2, 16.5410663134 19.66 ⋆ , 9.26 † , 7.18, 17.58 ‡† ¶ , 11.3312208602 24.81 ⋆ , 10.29, 2.1, 16.53 ⋆ ‘medium’ Rolling Band & Moire † ‘high’ Rolling Band & Moire ‡ Out of Spec Read Noise § Out of Spec Undershoot k Out of Spec Gain ¶ CCD Failed - no dataThe AGN are grouped based on light curve features as in Ta-ble 1. Col. 2 lists the module category 1 AGN above and the category 3 AGN discussedbelow. The light curves of these AGN are shown in the2 nd row of Fig. 3 between the light curves of the category1 and 3 AGN to facilitate comparisons. Although thereare indications of oscillatory behavior in this intermediatecategory of objects, the oscillations have very poorly definedtime-periods as compared to objects in category 3. However,their light curves are not as purely stochastic looking as thelight curves of the objects in category 1. Supporting thisnotion of an intermediate state is the observation that theobject kplr012158940, which is grouped with the category1 AGN, appears to be in the process of switching betweencategories as it begins to display features more reminiscentof category 3 light curves toward the end of the data.Category 3 AGN (kplr002694186, kplr009215110,kplr010841941, kplr003337670, and kplr007610713) areshown in the 3 rd row of Fig. 3. These objects appear toexhibit pronounced rippling features in the light curve.The aforementioned ‘moire’ and ‘rolling-band’ effects thatare known to exist in specific detector modules and chan-nels in Kepler are the natural suspects. These effects areknown to occur on timescales of a few days and have pat-terns very similar to those observed in these light curves(Kolodziejczak et al. 2010). However, a close examinationof the PDCSAP data suggests that the oscillations observedin these light curves are seen even when the target lands on c (cid:13) , 1–20 Vishal P. Kasliwal, Michael S. Vogeley, and Gordon T. Richards
300 600 900 kplr006932990Sy1z=0.0275
300 600 900 kplr012158940Sy1z=0.0670
300 600 900 kplr011178007Sy1z=0.0790
300 600 900 kplr009650715Sy1z=0.1270
300 600 900 kplr002837332Sy1z=0.1300
300 600 900 kplr003347632Sy2z=0.1580
300 600 900 kplr005781475Sy1z=0.2220
300 600 900 kplr002694186Sy1z=0.0890
300 600 900 kplr009215110Sy1.9z=0.1150
300 600 900 kplr010841941Sy1z=0.1520
300 600 900 kplr003337670Sy1z=0.3680
300 600 900 kplr007610713QSOz=0.4390
300 600 900 kplr006690887BL-Lacz=0.1540
300 600 900 kplr006595745QSOz=0.5020
300 600 900 kplr005597763QSOz=0.6250
300 600 900 kplr011606854FSRQz=0.9260
300 600 900 kplr010663134FSRQz=1.5200
300 600 900 kplr008703536Sy2z=0.0539
300 600 900 kplr011021406FSRQz=0.5130
300 600 900 kplr012208602QSOz=1.0980
Figure 3. ‘Stitched’ light curves of all the AGN in our investigation. The AGN are grouped based on the light curve categories discussedin Sec. 6. Light curves within a category are sorted in order of increasing redshift. All 20 light curves are plotted over the same rest-frametime range on the x-axis. Refer to Table 1 for details on individual objects. a detector module and channel combination that is knownto be free of the moire and rolling-band effects. For exam-ple, kplr002694186 exhibits strong oscillatory features in thelight curve as seen in Fig. 4. As reported in Table 2, thisAGN falls on the module
Kepler
Instrumentation Handbook). Yet the sameoscillatory signal is observed even during the quarters whenthe target falls on the un-affected module c (cid:13) , 1–20 re the Variability Properties of the Kepler AGN Light Curves Consistent with a Damped Random Walk? obs (d) F l u x ( e − ) Quarter Start Quarter Finish
200 400 600 800t rest (d) F l u x ( e − ) Figure 4.
PDCSAP light curve of the Seyfert 1 kplr002694186— Top panel: unstitched; bottom panel: stitched. The oscillatoryfeatures in this light curve may at first be suspected to originatein the known
Kepler ‘moire’ and ‘rolling-band’ instrumentationeffects. This AGN lands on the repeating module these only 24.81 suffers from moderate amounts of rolling-band and moire while the other 3 channels are clean. Weobserve (Fig. 2) that the rapid low amplitude oscillationsare present in the target during the first observed quarterwhen the target fell on the module rest (d) F l u x ( e − ) kplr011606854FSRQz =0.9260 obs. LC Figure 5. ‘Stitched’ PDCSAP light curve of the FSRQkplr011606854. After a featureless initial period, a burst of semi-oscillatory behavior occurs at the ∼
75 d mark culminating in theflare at the ∼
200 d mark. The light curve exhibits more oscilla-tory behavior after the flare but ultimately appears to settle backinto an almost featureless phase by the ∼
300 d mark suggestingthat at least some AGN may display ‘on’ and ‘off’ phases in theirlight curve. also has an oscillatory nature but is usually of larger am-plitude and longer period, as described in the case ofkplr003337670. Occasionally these punctuating periods re-sult in overall changes in the flux level of the object. Wecaution that these oscillatory features may not be genuinelyoscillatory, i.e. there are classes of ARMA random pro-cesses that can generate superficially oscillatory behaviorgiven appropriate ARMA parameter choices–see examplesin Woodward, Gray, & Elliot (2012).Category 4 AGN (kplr006690887, kplr006595745,kplr005597763, kplr0011606854, and kplr010663134) exhibitflares in their light curves as seen in row 4 of Fig. 3. Threeout of the five members of this class are blazars, implyingthat this type of behavior may be characteristic of blazars.However, no flares can be observed in the 4 th blazar in oursample, kplr011021406 (which is also one of the objects thatexhibits no variability). A solution to this puzzle may liein the light curve of kplr011606854 in Fig. 5. We see thatafter a featureless initial period, there is a burst of semi-oscillatory behavior at the ∼
75 d mark culminating in theflare at the ∼
200 d mark. The light curve exhibits moreoscillatory behavior after the flare, but ultimately appearsto settle back into an almost featureless phase by the ∼ c (cid:13) , 1–20 Vishal P. Kasliwal, Michael S. Vogeley, and Gordon T. Richards rest (d) F l u x ( e − ) kplr008703536Sy2z =0.0539 obs. LC Figure 6. ‘Stitched’ PDCSAP light curve of the Sy2 galaxykplr008703536. We detect no variability signal in this object -the observed variations in this light curve may be ascribed to in-strument noise suggesting that no endemic calibration problemsexist in the dataset.
FSRQ/Sy 1.5 kplr011021406. The absence of features inthese PDCSAP light curves and in a large number of stellarlight curves available on the MAST suggests that the latestPDC module (Data Release 21 onwards) corrects most ofthe systematic errors in the
Kepler
SAP light curves, leav-ing behind no endemic spacecraft event related artifacts.The only remaining sources of non-physical variability suchas moire patterns, etc. must therefore be restricted to spe-cific pixels and CCDs, allowing us to assume with high con-fidence that our stitched light curves are primarily astro-physical signal. The lack of variability observed in these twolight curves confirms that not all AGN exhibit variability(Sesar et al. 2007), though our sample is too small to beable to put constraints on the fraction of AGN that exhibitvariability. The QSO kplr012208602 exhibits a very weakvariability signal, occasionally displaying weak, intermittentsemi-sinusoidal changes in flux. We posit that these changeswould be un-observable using ground-based studies due tothe weakness of the variability signal. The light curves forall 3 objects are plotted in the last row of Fig. 3.We caution that per-channel/per-pixel effects may beresponsible for some of the features observed in our lightcurves, particularly amongst the category 2 light curves.A true determination of the reality of some of these fea-tures will have to await the very systematic, customizedper-object calibration discussed at the beginning of Sec. 5.In the next section we discuss our method for analyzing themulti-quarter light curves.
We probe the variability properties of the
Kepler
AGN lightcurves for consistency with the AR(1) process of Sec. 2 by using structure functions to determine best-fit values of theparameter γ in the damped power-law (DPL) model of Eq.(9). The DPL model is used to generate mock light curveswith the same cadence and missing observation propertiesas the AGN under investigation. We estimate best fit modelparameters by comparing the structure function of the ob-served light curve to the ensemble of structure functionscomputed from the mock light curves consistent with themodel parameters.Many definitions of structure functions have ap-peared in astronomy literature (see Graham et al.(2014) for an overview). We use the definition pro-vided by Simonetti, Cordes, & Heeschen (1985) andRytov, Kravtsov, & Tatarskii (1987). The n-th OrderStructure Function of the process F ( t ) is SF ( n ) F (∆ t ) = h Ξ ( N ) F ( t, ∆ t ) i t , (15)where Ξ ( N ) F ( t, ∆ t ) is the N-th increment of the process F ( t )at time-lag ∆ t . The N-th increment is given byΞ ( N ) F ( t, ∆ t ) = ∞ X m =0 ( − m (cid:18) Nm (cid:19) F ( t + [ N − m ]∆ t ) . (16)Using Eq. (15) and Eq. (16), the 1 st increment and 1 st OrderStructure Function of the process F ( t ) (hereafter referred toas ‘the increments’ and ‘the SF’) may be estimated usingΞ( t, ∆ t ) = F ( t + ∆ t ) − F ( t ) , (17)and SF (∆ t ) = h [ F ( t + ∆ t ) − F ( t )] i . (18)Structure functions offer benefits over directly estimatingthe PSD and ACVF/ACF of the process. Unlike the PSD,structure functions may be estimated in real space as op-posed to Fourier space, making them more robust estimatorsof model parameters than PSD estimators that suffer fromwindowing and aliasing concerns (Rutman 1978). Structurefunctions offer an advantage over estimating the ACVF andACF because of the de-trending properties of structure func-tions: an n-th Order Structure Function is insensitive to (n-1)-th order trends in the dataset. Unfortunately, the use-fulness of this property is limited to the lower order struc-ture functions. The Kepler light curves are not long enoughfor computation of significantly higher order structure func-tions. Structure functions of all orders are related to theACVF by simple linear equations. For example, the first or-der structure function can be related to the ACVF of F ( t )by SF (∆ t ) = 2 ACV F (0) − ACV F (∆ t ) . (19)To understand what information the structure functionconveys, consider the histograms of flux differences, i.e. 1 st increments for various values of time-lag ∆ t shown in Fig. 7.By definition, the structure function at time-lag ∆ t is the av-erage of the squares of the values entering each histogram,i.e. it is the variance of the flux differences at the chosentime-lag ∆ t . The structure function therefore quantifies howthe variance of the flux differences changes as a function oftime-lag. In the case of objects that show stochastic varia-tions in F ( t ), the histograms look like bell-shaped curves atsmall time-lags. As the time-lag increases, the width of thebell-curves increase until a critical time-lag τ is reached. In c (cid:13) , 1–20 re the Variability Properties of the Kepler AGN Light Curves Consistent with a Damped Random Walk? l og δ t r e s t ( l og d ) −2 0 2 δ F ( e − ) −2000−1000 0 1000 C o un t Figure 7.
Histogram plots of the 1 st increment ( F ( t +∆ t ) − F ( t ))for various choices of the time-lag ∆ t drawn from the light curveof the Seyfert 1 galaxy kplr006932990. Each histogram shows thedistribution of the 1 st increment values for a given choice of timelag ∆ t . The 1 st -order structure function is the variance of thedistribution in each histogram and quantifies how the variance ofthe increments changes as a function of time lag. the case of the AR(1) process, this critical time-lag can beidentified with the the turnover timescale in Eq. (3).Astronomical time series have measurement noisethat is usually modeled as uncorrelated white-noise, i.e. ACV F (0) = σ + σ N where σ is the contribution to thevariance of the time series at zero-lag from the actual signal,while σ N characterizes the noise level of the measurementnoise. In the case of the AR(1) process, the ACVF takes theform in Eq. (7). If the variability observed in Kepler
AGNlight curves is well described by an AR(1) process, then thestructure function should be well described by SF DRW (∆ t ) = 2 σ (1 − e − | ∆ tτ | ) + 2 σ N . (20)This may be tested by estimating the structure function of areal AGN light curve to check if it is consistent with Eq. (20)or with a more general form based on the DPL PSD of Eq.(9). Unfortunately there is no closed-form expression for theauto-correlation function corresponding to the DPL PSD,making it impossible to directly fit the functional form forthe structure function observed for a given AGN. In the nextsection, we describe a Monte-Carlo estimation technique forrecovering the DPL model parameters. We estimate the DPL parameters in Eq. (9) via Monte-Carlosimulations. We compare the real structure function of anAGN to simulations computed from mock light curves gen-erated using the DPL model of Eq. (9). We generate ‘mock’light curves using the Timmer & K¨onig (1995) method. Tocreate a single mock light curve, pseudo-random numbersare generated using the Fast Mersenne Twister SFMT19937generator seeded with hardware-generated random numbers(generated using Intel RDRAND instruction) to ensure that rest (d) F l u x ( e − ) kplr006932990Sy1z =0.0275 obs. LCeg. mock LC Figure 8.
Observed light curve of the Sy 1 AGN kplr006932990(orange) along with an example mock light curve (light green)generated as described in 8 (generated with γ = 2 . τ = 27 . the random number sequences are free of artificial corre-lations (Coddington 1994) induced by poor random seedchoices. At this intermediate stage, the mock light curve isoversampled by a factor of 10 i.e. we generate points at 10 × the required cadence in order to avoid sampling issues. Toinclude low frequency modes that are not adequately char-acterized by the length of the observed light curve, the inter-mediate mock light curve is much longer than is ultimatelyrequired to make the final mock; mock light curves gener-ated in this manner are capable of exhibiting low frequencymodes longer than the length of the observed data. FFTsare most efficient for data sequences that are a power of 2;for this reason, we pick the intermediate (including the over-sampling) to be of length 2 . This results in the intermedi-ate mock light curve being between 15 × to 45 × the lengthrequired for the final mock light curve depending on the ac-tual length of the observed light curve. We pick a uniformlydistributed random segment of the intermediate overly-longlight curve that has the same length as the observed lightcurve and generate another stream of un-correlated Gaus-sian random deviates to simulate the white-noise propertiesof Kepler instrumentation noise. After adding this ‘measure-ment noise’, we set data points corresponding to the un-observed cadences in the real light curve to 0. This proce-dure creates a final mock light curve with identical samplingand noise properties to the real light curve. Fig. 8 shows thetrue light curve (orange) along with an example mock lightcurve (light green) for the Sy 1 AGN kplr006932990 illus-trating what the mock light curves look like for the best fitDPL parameters for this object.We compute the structure functions using Eq. (18) c (cid:13) , 1–20 Vishal P. Kasliwal, Michael S. Vogeley, and Gordon T. Richards − ] )1.1e−052.3e−053.4e−05 N o r m . C o un t s SW: 7.94e-01; p-val: 9.40e-34 SF (log [e − ] )0.360.731.09 N o r m . C o un t s SW: 9.95e-01; p-val: 4.88e-03
Figure 9.
Distribution of mock SF ( δt ) estimates for 1000 simu-lated light curves of length T = 10 d, with γ = 2 . τ = 7 . δt = 5 d. The upper panel shows the (clearly non-Gaussian) dis-tribution of SF ( δt ), while the lower panel shows the distributionof log SF ( δt ). The Shapiro-Wilk test W- & p-values confirmsthat the distribution of log SF ( δt ) is closer to Gaussian thanthat of the raw SF ( δt ). modified to account for missing values - SF n = P N − ni =1 w i w i + n [ F i + n − F i ] P N − ni =1 w i w i + n , (21)where n = ∆ t/δt is the lag in cadences as opposed to physi-cal time. The weights are defined to be w i = 1 for observedcadences and w i = 0 otherwise.It is known that, for any given lag ∆ t , the distribu-tion of SF (∆ t ) generated using the Timmer-K¨onig methodis not Gaussian (Emmanoulopoulos, McHardy, & Uttley2010). However, the logarithms of the structure function,log SF (∆ t ), are distributed as per a Gaussian distribution.As can be seen in Fig. 9, histograms of mock SF (∆ t ) andlog SF (∆ t ) values for a set of 10000 mock light curves con-structed with PSD slope γ = 2, and characteristic timescale τ = 7 . d at ∆ t = 5 d for mock light curves of length T = 10 dsuggest that while the estimates of SF (∆ t ) are not Gaussiandistributed, the estimates of log SF (∆ t ) are. This obser-vation may be confirmed using a test such as the Shapiro-Wilk test. Note however, that even though the Shapiro-Wilktest confirms that log SF (∆ t ) is closer to Gaussian than SF (∆ t ), it is not always Gaussian in an absolute sense; amore through investigation of the distribution properties ofthe structure function may be warranted. We use log SF n instead of SF n to compare structure functions between ob-servations and simulations.Since the logarithms of the structure function estimatesare drawn from Gaussian distributions, we can use the meanand standard deviation vectors of the logarithms of the mockstructure functions to estimate the best-fit DPL model pa-rameters. We compute the mean µ n = h log SF n i and vari- −1 0 1 2log ∆t rest (d) l og S F ( l og [ e − ] ) kplr006932990Sy1z =0.0275γ =2.6τ =27.5 (d)χ =1.43P =83.3% model SF µmodel SF σobs. SFeg. mock SF Figure 10.
Observed structure function of the Sy 1 AGNkplr006932990 (orange), best-fit model structure function withstandard deviations (purple) and the mock structure function(light green) corresponding to the mock light curve in Fig. 8.The value of γ is most sensitive to the short timescale (∆ t . ance σ n = h (log SF n − h log SF n i ) i of the ensemble ofthe logarithms of the mock structure functions and numeri-cally minimize χ = N − X n =1 ( µ n − log SF n,obs ) σ n , (22)using the Constrained Optimization BY Linear Approxi-mations (COBYLA) optimization algorithm (Powell 1994)to search the parameter space for the best-fit estimates of γ , τ , σ S , and σ N . Coarse optimization is performed using10 mocks at each step in the optimization, after which thebest-fit parameter estimates are refined using 10 mocks ateach step in the optimization process. Such large numbersof mocks are required to ensure that χ values computed atthe same point in parameter space (using different choicesof random seeds) are within the tolerances chosen in theoptimization step.The resulting values of the reduced chi-square perdegree-of-freedom, χ DoF are close to 1. We also computethe percentile, P , of the value of χ DoF corresponding to thebest-fit parameter values using a fresh set of 1000 mock lightcurves. Percentile values P close to 100% mean that the val-ues of χ mock are almost uniformly smaller than χ obs for thebest-fit model parameters, indicating that the model is apoor fit to the data.Fig. 10 shows the observed structure function of theSeyfert 1 kplr006932990 (orange) along with the best-fitmock structure function and standard deviation (purple)and the mock structure function (light green) correspond-ing to the mock light curve in Fig. 8. On short timescales(∆ t . c (cid:13) , 1–20 re the Variability Properties of the Kepler AGN Light Curves Consistent with a Damped Random Walk? standard deviation, while on longer timescales, the structurefunction is allowed to vary more significantly. The value of γ is determined almost exclusively by this short timescale(∆ t . γ > .
0, as the PSD turns over from shorter timescales(∆ t / γ >
2) to longer timescales(∆ t '
10 d with local PSD slope γ → local valueof the PSD slope passes through γ ≡
2. Most ground-basedsurveys, such as the SDSS, have large (compared to
Kepler )photometric errors and are measuring the structure functionon timescales ∆ t ' γ ∼ χ taking ∼ γ we generated simulations of short segments of mock lightcurves with known DPL model parameters and then pro-ceeded to recover the input parameters. Based on this, weestimate that the error in the value of γ obtained is on theorder of 10%. T MIN
We estimate the shortest timescales on which we observevariability in the light curves. Our estimation of the DPLmodel parameters in Sec. 8 yields an estimate of the un-correlated white-noise level present in each light curve. Thisnoise level can be seen in the structure function in the formof a flattening-out of the structure function on very shorttimescales. For example, in Fig. 11, we observe that thestructure function of the Seyfert 1 kplr006932990 begins tolevel off to a value of ∼ − ] on timescales of ∼ . σ N ∼ − ] . By removing this noise level fromthe observed structure function, we can estimate a noise-less version of the structure function. To estimate the short-est timescale on which we observe variability in the lightcurve, we find the time-lag at which this noiseless structurefunction crosses the noise floor. In Fig. 11 we plot the ob-served noisy structure function of kplr006932990 (orange),our best fit estimates of this structure function (purple),and the noise level determined from the structure functionfit (red dashed-line). After removing this noise level from theobserved structure function, we obtain the noiseless struc-ture function (red). We also show the mock structure func-tion (light green) corresponding to the mock light curve inFig. 8 along with the noiseless version (green) of this mockstructure function. We define the variability onset timescale T min to be the intersection point (indicated by blue dashed-line) of the noiseless structure function with the noise floor.We compute the value of the variability onset timescale forthe AGN in this study along with the rest-frame samplinginterval for each object and present the results in Table 3. −1 0 1 2log ∆t rest (d) l og S F ( l og [ e − ] ) kplr006932990Sy1z =0.0275γ =2.6τ =27.5 (d)χ =1.43P =83.3% model SF µmodel SF σobs. SFobs. SF - noise rem.eg. mock SFeg. mock SF - noise rem.noise flrT min Figure 11.
Observed structure function of the Sy 1 AGNkplr006932990 (orange) and the best-fit model structure functionwith standard deviation (purple). The noiseless structure function(red) is obtained by removing the noise floor (red dashed-line)after which the shortest variability timescale T min (blue) is de-termined from the intersection of the noiseless structure functionwith the noise floor. T min typically <
10 OBSERVED STRUCTURE FUNCTIONSAND FITS
The excellent short timescale sampling properties of
Kepler make it possible to study the AGN structure functions withunprecedented levels of detail. Using the Monte-Carlo fittingmethod of Sec. 8, we fit the DPL parameters of Eq. (9) forthe 17 most variable objects. Fig. 12 shows the observedstructure functions (orange), Monte-Carlo fits (purple line),and 1 σ variation (purple shaded region).Sampling information and estimates of the DPL param-eters are given in Table 3 for each Kepler
AGN. Column 1lists the
KeplerID of the object. Column 2 lists the lengthof each light curve in quarters. Each quarter spans about 3months in duration and has ∼ γ and τ from Eq.(9) where τ is measured in rest-frame days. Column 6 liststhe best-fit estimate for the minimum timescale on whichvariability is observed T min in minutes. The ‘goodness-of-fit’, i.e. the reduced chi-square, is listed in Column 7 whileColumn 8 lists the percentile of the chi-square of the ob-served structure function fit as compared to simulationsdrawn from the DPL model i.e. it is also an indicator of‘goodness-of-fit’ in that a very large value of the chi-square c (cid:13) , 1–20 Vishal P. Kasliwal, Michael S. Vogeley, and Gordon T. Richards −1 0 1 2 kplr006932990Sy1z=0.0275 γ=2.6τ=27.5 (d)χ =1.43P=83.3% −1 0 1 2 kplr012158940Sy1z=0.0670 γ=2.5τ=40.2 (d)χ =1.15P=71.3% −1 0 1 2 kplr011178007Sy1z=0.0790 γ=2.9τ=16.8 (d)χ =0.76P=38.5% −1 0 1 2 kplr009650715Sy1z=0.1270 γ=2.2τ=47.8 (d)χ =0.64P=40.2% −1 0 1 2 kplr002837332Sy1z=0.1300 γ=2.1τ=41.0 (d)χ =0.86P=55.4% −1 0 1 2 kplr003347632Sy2z=0.1580 γ=2.5τ=18.0 (d)χ =1.07P=68.1% −1 0 1 2 kplr005781475Sy1z=0.2220 γ=3.0τ=8.1 (d)χ =1.18P=75.3% −1 0 1 2 kplr002694186Sy1z=0.0890 γ=2.1τ=134.3 (d)χ =1.72P=88.6% −1 0 1 2 kplr009215110Sy1.9z=0.1150 γ=2.2τ=55.1 (d)χ =3.22P=98.6% −1 0 1 2 kplr010841941Sy1z=0.1520 γ=2.3τ=10.9 (d)χ =1.23P=78.3% −1 0 1 2 kplr003337670Sy1z=0.3680 γ=2.2τ=125.4 (d)χ =0.16P=1.8% −1 0 1 2 kplr007610713QSOz=0.4390 γ=2.5τ=29.3 (d)χ =0.63P=35.2% −1 0 1 2 kplr006690887BL-Lacz=0.1540 γ=2.2τ=3.0 (d)χ =13.37P=99.9% −1 0 1 2 kplr006595745QSOz=0.5020 γ=1.7τ=76.8 (d)χ =1.24P=76.9% −1 0 1 2 kplr005597763QSOz=0.6250 γ=1.8τ=12.5 (d)χ =3.14P=98.7% −1 0 1 2 kplr011606854FSRQz=0.9260 γ=3.0τ=1.6 (d)χ =8.85P=99.9% −1 0 1 2 kplr010663134FSRQz=1.5200 γ=2.1τ=21.7 (d)χ =5.47P=99.9% −1 0 1 2 kplr008703536Sy2z=0.0539 −1 0 1 2 kplr011021406FSRQz=0.5130 −1 0 1 2 kplr012208602QSOz=1.0980
Figure 12.
Observed structure functions (orange) and fits using mock light curves (purple) of all the AGN. Note that all the structurefunctions have been plotted at the same scale. percentile indicates that the DPL model is insufficient toexplain the data.We observe that in all but one case (the blazarkplr00669887),
Kepler samples the light curve on shortersampling intervals δt rest (Col. 4 in Tab. 3) than T min (Col.7). The median variability onset timescale is ∼
23 h in therest frame of the object, while the minimum and maximumvalues are ∼ . ∼ . Kepler in the rest-frame of this object. This suggests that the
Kepler sam-pling pattern is more than adequate to study optical AGNvariability. We caution that T min is merely an upper esti-mate for the shortest timescale on which variability can beobserved given the noise level, therefore even better photo-metric precision than is available with Kepler is required tofind the true shortest timescale on which AGN vary.We categorize the light curves of our AGN based on thepresence of oscillatory features and flares in Sec. 6. All 3 ob-jects that we grouped in category 1 (completely stochastic-looking light curves–row 1 of Fig. 3) are fit by DPL models c (cid:13) , 1–20 re the Variability Properties of the Kepler AGN Light Curves Consistent with a Damped Random Walk? Table 3.
Sampling Parameters & Results.Identifier Sampl. Param. Analysis ResultsKeplerID L δt rest γ τ T min χ DoF P ( ⋆
11 27.02 2.1 134.3 2339 1.72 88.69215110 8 26.39 2.2 55.1 1704 3.22 98.610841941 7 25.54 2.3 10.9 2151 1.23 78.33337670 7 21.51 2.2 125.4 9946 0.16 1.807610713 8 20.45 2.5 29.3 1677 0.63 35.26690887 7 25.50 2.2 3.1 - 13.37 99.96595745 7 19.59 1.7 76.8 711 1.23 76.95597763 7 18.11 1.8 12.5 1151 3.14 98.711606854 12 15.28 3.0 1.6 2005 8.85 99.910663134 12 11.68 2.1 21.7 463 5.47 99.98703536 12 27.92 - - - - -11021406 12 19.45 - - - - -12208602 12 14.02 - - - - - ⋆ kplr002694186 has a total light curve spanning 17 quarters inthe MAST database. Data is available only for quarter 0 (as anexo-planet search program target) and then subsequently fromquarters 7 through 17 (as a GO program target). We discard thequarter 0 light curve segment because of the unreliability of thephotometry during this initial phase of the mission and reportkplr2694186 as having a light curve of length 11 quarters.The AGN are grouped based on light curve features as in Table1. Column 2 (‘ γ is the best fit slope of thepower-law portion of the power spectral density, τ is the turnovertimescale beyond which individual data points are essentially un-correlated, and T min is the shortest timescale on which variabilityis observed. Columns 7 & 8 list measures of the ‘goodness-of-fit’i.e. the reduced χ and the percentile value P of the chi-squarecompared to chi-squares computed from 1000 mock light curves. with median γ = 2 . χ DoF close to 1. This value of γ is significantly different from that expected of an AR(1)process and indicates that the light curves of these 3 objectsare smoother than a damped random walk ( γ DRW = 2 . st row of Fig. 12 is the presence of a‘dip’ feature on time-lags ∆ t between ∼
10 and 100 d. Thisdip feature may be interpreted as an excess of correlation onthese timescales. The dip starts at ∼
10 d in the rest-frame of the objects, suggesting that it is probably intrinsic to thelight curves of these objects. Such dips rarely occur in the −1 0 1 2log ∆t rest (d) − l og S F ( l og [ e − ] ) kplr009215110Sy1.9z =0.1150γ =2.2τ =55.1 (d)χ =3.22P =98.6% model SF µmodel SF σobs. SFobs. SF - noise rem.eg. mock SFeg. mock SF - noise rem.noise flrT min Figure 13.
Observed structure function of the Sy 1.9 AGNkplr009215110 (orange), noise-removed observed structure func-tion (red), example mock structure function (light green), noise-less mock structure function (green), and best-fit model structurefunction with standard deviations (purple). This structure func-tion exhibits varying PSD slope γ . The structure function risessteeply when ∆ t rest . structure functions of the simulated light curves (generallythese exhibit wiggling behavior after the turnover timescaleis reached). This suggests that, even though all three objectshave χ values consistent with the DPL model ( P ranges be-tween 40% and 80%), the DPL model is probably too simpleto characterize the variability observed in these objects.Of the objects categorized as intermediate between cat-egories 1 and 3 (stochastic-looking light curves+oscillatoryfeatures–row 2 of Fig. 3), the first 2 have PSD slopes ( γ = 2 . γ = 2 . χ DoF = 0 .
64 implying that the DPL modelmay be overfitting the light curve of this object. The lat-ter two objects, kplr003347632 & kplr5781475, have muchsteeper PSD slopes ( γ = 2 . γ = 3 . γ = 2 .
35 and the individual structure functions are shownin row 2 of Fig. 12.Objects that fall into category 3 (strong oscillatory fea-tures in light curves–row 3 of Fig. 3) based on a visual in-spection of the light curve have the largest number of lightcurve slopes consistent with that expected from the AR(1)process. Four out of 5 have γ ∼ .
2, while the median esti-mate of γ = 2 . γ range between 2 . . τ for all of theseobjects, a visual inspection of the corresponding structurefunctions in the row 3 of Fig. 12 suggests that the turnover c (cid:13) , 1–20 Vishal P. Kasliwal, Michael S. Vogeley, and Gordon T. Richards timescale is suspect in all but the case of kplr010841941, asborne out by the mock structure function error contours.Dips similar to those observed in category 1 objects are alsoseen in the structure function of some of these objects. As inthe case of category 1 AGN, the dip feature begins around ∼
10 d in the rest-frame of the objects and is usually ab-sent by ∼
100 d. kplr002694186 & kplr009215110 exhibitstructure functions that appear to have different slopes ondifferent timescales. This presence of multiple slopes maybe responsible for the relatively poor quality of the fits ofthe structure functions of these objects ( χ DoF is 1 .
72 and3 .
22 respectively). Fig. 13 shows the structure function ofthe Sy 1.9 AGN kplr009215110. On short timescales . γ ranging from a fairlyflat 1 . . χ DoF values such as 13 .
4, 8 .
9, &5 . χ values than is the normfor mock light curves generated using the DPL model, sug-gesting that the DPL model (and therefore the AR(1) pro-cess) are wholly unsuitable for modeling the light curves ofobjects that show flaring behavior.The last row of Fig. 12 shows the individualstructure functions of kplr008703536, kplr011021406, andkplr012208602. A visual inspection of the correspondinglight curves in Fig. 3 suggests that these objects do notexhibit significant levels of variability. This observationis borne out by the structure functions of these objects,which are mostly flat and featureless for the duration ofobserved time-lags. On long timescales ( ∼
100 d) thestructure function begins to exhibit low amplitude ‘wig-gles’. It is well known from the theory of ACF estima-tion (Brockwell & Davis 2002) that similar wiggles occur ontimescales comparable to the length of the time-series andare not significant.
11 COMPARISON
We observe that not all AGN exhibit DPL parameter γ ∼ Kepler
AGN.
Kepler
AGN Variability Studies
We have shown that the AGN in this study exhibit PSD withshapes that are superficially similar to that expected from asimple variability model, such as the damped random walkor AR(1) process—a flat power spectrum on long timescalesthat turns into a 1 /f γ power-law section after frequencieshigher than 1 /τ . However while the AR(1) requires the valueof γ to be exactly 2, we find a wide range of γ values incon-sistent with the AR(1) model. Both Carini & Ryle (2012) and Mushotzky et al. (2011) obtained values of the PSDslope for the Seyfert I AGN Zw 229-15 or kplr006932990.Carini & Ryle (2012) applied various PSD models and foundthat the best fit estimate of the PSD slope is γ ∼ . τ = 92 ± /
21 d for the ‘knee’ model to τ = 43 ± /
10 dfor their ‘broken power-law’ model, which is in good agree-ment with our estimate of γ = 2 . τ = 27 . γ = 3 . Kepler data product. Carini & Ryle (2012) were able to useground-based photometry from the reverberation mappingcampaign of Barth et al. (2011) along with a customizedaperture to stitch together the 4 quarters of light curve dataused in the study. Mushotzky et al. (2011) calculate PSDfits for 4 individual quarters of
Kepler data but use theun-processed SAP light curve in their estimate of the PSD.We used PDCSAP fluxes (Data Release 21+) from all 14available quarters of data in our analysis. These fluxes, cali-brated using the newest version of the PDC module, shouldsuffer only marginally from instrumentation issues. Thus, itis likely that the range of PSD estimates observed is a re-sult of the subtle but important differences between the dataproducts used. Given that all 3 studies have used differentanalysis techniques but arrive at similarly steep PSD slopeestimates, it is clear that the PSD slope is robustly steeperthan 1 /f and is irreconcilable with the AR(1) process. Mostrecently, Kelly et al. (2014) attempt to use the Kalman filterto apply a general ARMA model to a single quarter of datafrom kplr006932990 and find that the observed variability isbest-fit by an ARMA process with 6 Auto-Regressive and4 Moving Average components—a much more complicatedstochastic process than a simple damped random walk.Mushotzky et al. (2011) also obtained PSD slope es-timates for 3 other AGN: kplr012158940, kplr011178007,and kplr002694186. Their analysis was performed on theSAP fluxes observed during individual quarters. The av-erages of their estimates of the PSD slopes: γ = 2 . γ = 2 .
92 (kplr011178007) agree fairlywell with our estimates in Table 3; however, they observea much steeper γ = 2 .
845 for kplr002694186 than our es-timate of γ = 2 .
1. Both kplr012158940 and kplr011178007exhibit strong variability–most of the ‘signal’ observed inthe uncalibrated SAP light curves for both objects was in-trinsic to the source and retained in the calibrated PCD-SAP light curves. On the other hand, the calibrated PDC-SAP light curve of kplr002694186 lacks the large amplitudeflux variations seen in the original SAP light curve, indicat-ing that those variations were instrumental effects that mayhave inadvertently introduced power at lower frequencies inthe PSD computed by Mushotzky et al. (2011), artificiallymaking the PSD slope steeper than that computed in ouranalysis.Edelson et al. (2013) have examined the variabilityproperties of the light curve of the blazar kplr006690887. Us-ing segments of data spanning ∼ . c (cid:13) , 1–20 re the Variability Properties of the Kepler AGN Light Curves Consistent with a Damped Random Walk? power law PSD with slopes γ = 1 . ± .
14 and 2 . ± .
34 ifthe highest frequencies are excluded from the analysis. Us-ing the full light curve, we find that our fit is of very poorquality with γ = 2 .
2. However, flaring is not an AR(1) pro-cess. An AR(1) process is completely incapable of produc-ing the flare features observed in this AGN. Ground basedstudies of blazar variability, such as Ruan et al. (2012), haveattempted to model blazar variability with an AR(1) pro-cess using 101 blazar light curves drawn from the LINEARnear-Earth asteroid survey. The LINEAR survey collectedtime series data over a period of 5.5 yr with two 1.01 metertelescopes giving a r-band survey depth of 18 mag at 5 σ .Sources in the survey within 10 o of ecliptic have 460 obser-vations on average, while sources further from the eclipticonly have 200 observations over the duration of the survey,although the cadence can occasionally be as high as onceevery 15 min. While Ruan et al. (2012) find that the blazarlight curves observed by them are well fit by the AR(1)process with decorrelation timescales ranging from < ∼ ∼
100 d, we caution that theLINEAR survey suffers from both highly irregular, sparsesampling patterns, and large photometric uncertainties.Wehrle et al. (2013) and Revalski et al. (2014)have estimated the PSD slopes of the radio-loudAGN kplr011021406, kplr011606854, kplr01228602, andkplr010663134 in Table 3. Wehrle et al. (2013) perform aPSD analysis of the un-corrected SAP flux light curve ofeach object on a per quarter basis. They also supply thesame analysis for ‘corrected’ and ‘end-matched’ versionsof the light curve for each quarter, where they correctedthe SAP data by removing some of the instrumental effectsand windowed the data by removing a linear trend fromthe data to make the first and last point of each quarterhave identical flux value. In the case of kplr011021406,they concluded that this Seyfert 1.5/FSRQ is variable andderive a mean PSD slope γ = 1 . .
2. Individualestimates of their PSD slopes for each quarter range from1 . . γ = 2 . γ = 3 . γ = 1 . γ = 1 . γ = 1 . Previous variability studies using ground-based data suggestthat the AR(1) process is adequate to characterize AGNvariability given the measurement uncertainties and sam-pling pattern typical of ground-based studies. The originalestimation performed by Kelly, Bechtold, & Siemiginowska(2009) using the standard state-space representation of theAR(1) process (Brockwell & Davis 2002) indicated that theR- and B-band optical variability observed in their sam-ple of 109 quasars and Seyferts was well modeled by AR(1)processes with a range of decorrelation timescales stronglyindicative of a connection between thermal fluctuations oc-curring on the AGN accretion disk and the observed vari-ability. The 109 AGN used in this study include 59 MACHOquasars observed with sampling rates of once every 2–10 dover 7.5 yr, 42 PG quasars observed once every 40 d over7 yr, and 8 Seyferts observed once every 1–10 d resultingin their data set having no cadence higher than ∼ − .The best-fit characteristic timescales for their objects rangefrom ∼ ∼ ∼ Kepler structure functions in Fig. 12, thisdataset and other ground-based datasets do not sample thestructure function very well below ∼
10 d making these stud-ies less sensitive to the short timescale ( . ∼ ∼ M i . −
24) thanthe AGN observed by
Kepler . MacLeod et al. (2010) usedthe PRH algorithm to model the light curves of the ∼ c (cid:13) , 1–20 Vishal P. Kasliwal, Michael S. Vogeley, and Gordon T. Richards
Kozlowski et al. (2010); most quasars have decorrelationtimescales in the range of ∼ − ∼
50 yr baseline with a handful of points per light curve)with the Stripe 82 data by MacLeod et al. (2012). These au-thors found that the AR(1) process continued to fit the ob-served light curves well and yielded decorrelation timescalesin the range of 5 − γ -ray AGN variability observedby the Fermi space telescope.Recently Zu et al. (2013), Andrae, Kim, & Bailer-Jones(2013), and Graham et al. (2014) have rigorously testedthe AR(1) process on the timescales probed by the sur-veys used in the previous studies. Using the PRH formal-ism adopted by Kozlowski et al. (2010), but using moregeneralized forms for the structure of the covariance ma-trix employed, Zu et al. (2013) concluded that while theOGLE light curves used in the study are consistent withthe AR(1) process on the longest timescales probed by theOGLE data ( ∼ /f PSD. Furthermore, they concluded that while the AR(1)process produces acceptable fits to the light curves, on thewhole the simple models used in that study are not well con-strained by the light curves. Zu et al. (2013) concede thatthe quality of the photometry may be responsible for themixed results obtained by them when fitting the various co-variance matrix models, i.e the scatter observed by them inthe slope parameter fit may be either intrinsic to the lightcurves and evidence for varied behavior between the ob-jects, or a consequence of errors in the determination of thephotometric uncertainties. A comprehensive and systematicBayesian study of SDSS Stripe 82 quasar light curves per-formed by Andrae, Kim, & Bailer-Jones (2013) concludedthat more sophisticated models, such as AR(2), GARCH(1),and ARMA(1,1), are not required to model light curves. Ofthe 6304 quasars examined in this study, an overwhelmingmajority–5047 light curves–were best modeled by an AR(1)process while only a handful required a more sophisticatedapproach.While a large number of studies have indicated thatthe AR(1) process is sufficient to model AGN variability,these studies have used comparatively sparsely sampled lightcurves (10 −
100 points over several years) that do not wellsample the light curve on short timescales (∆ t . SF (∆ t i,j ) = h r π | ∆ m i,j | − q σ i + σ j i ∆ t , (23)where ∆ m is the magnitude difference over the time interval ∆ t . While our structure functions is quadratic in flux differ-ence, the structure function of Schmidt et al. (2012) is linear in magnitude difference. Our correlation strength parameter γ appears in the light curve PSD model, while the correla-tion strength parameter λ of Schmidt et al. (2012) appearsin the structure function model making it impossible to di-rectly compare our γ with their λ . However, Schmidt et al.(2012) find that Stripe 82 quasars are fit by λ values rang-ing from 0.0 to 1.2 with a peak at λ = 0 .
25. The rangeof λ values found by Schmidt et al. (2012) suggests that asimple stochastic process model like the AR(1) process is un-likely to work well for all objects. Graham et al. (2014) usea wavelet analysis technique known as the Slepian WaveletVariance (SWV) to study 18000 quasars from the CatalinaReal-time Transient Survey (CRTS) and SDSS Stripe 82.By comparing the expected Slepian Wavelet Variance froman AR(1) process to that observed for real quasars in S82and CRTS, Graham et al. (2014) find that there is a clearindication that the AR(1) process fails to characterize AGNvariability on short timescales.While one cannot ignore the success of the AR(1) pro-cess at modeling the time variability behavior of AGN on thetimescales probed by MACHO, OGLE, SDSS, etc., we cau-tion that it may be impossible to observe deviations from theAR(1) process without probing the short timescales avail-able through Kepler .
12 CONCLUSIONS
Individual
Kepler
AGN light curves show a wide rangeof behavior that can be loosely grouped into 5 categories:stochastic-looking, somewhat stochastic-looking+weak os-cillatory features, oscillatory features dominant, flare fea-tures dominant, and not-variable. Some light curves appearto transition from one variability state to another, suggest-ing that AGN light curves may not be stationary in the sta-tistical sense (Sec. 6 and Fig. 3). We estimate PSD slopesranging from γ = 1 . γ = 3 . γ ≡ . γ = 2 . γ < .
9. Four out of 5 of the AGNthat exhibit oscillatory features in the light curve have PSDslope close to that of the AR(1) process, while the AGN thatexhibit flares in the light curve are poorly fit by the DPLmodel (Sec. 10 and Table 3).A broad dip feature can be observed in the structurefunction on timescales ranging from ∼ ∼ γ in the DPL model varies on differenttimescales (see Fig. 13). Lastly, not all AGN show variabilityat the levels probed by Kepler , confirming previous findings c (cid:13)000
9. Four out of 5 of the AGNthat exhibit oscillatory features in the light curve have PSDslope close to that of the AR(1) process, while the AGN thatexhibit flares in the light curve are poorly fit by the DPLmodel (Sec. 10 and Table 3).A broad dip feature can be observed in the structurefunction on timescales ranging from ∼ ∼ γ in the DPL model varies on differenttimescales (see Fig. 13). Lastly, not all AGN show variabilityat the levels probed by Kepler , confirming previous findings c (cid:13)000 , 1–20 re the Variability Properties of the Kepler AGN Light Curves Consistent with a Damped Random Walk? that flux variability is seen in most but not all AGN (seeFig. 12).There is no clear relationship between the PSD slopeand object type or redshift, although the sample size istoo small to draw any definitive conclusion. We concludethat while enough of our light curves are consistent withan AR(1) process for the damped random walk model tobe an appealing choice–especially for poorly sampled lightcurves–there are enough interesting features present in thelight curves to warrant a more detailed analysis.The AGN light curves seen in Kepler suggest that AGNvariability is a very complex phenomenon with individuallight curves looking very different. We will perform a recali-bration of the
Kepler
AGN light curves to determine whichlight curve visual features are persistant. We will apply moresophisticated statistical models drawn from the field of timeseries analysis such as the ARMA process model Hamilton(1994). The AR(1) process or DRW model can determine ‘if’an AGN is varying, but is not helpful in determining ‘why’the AGN is varying.
ACKNOWLEDGMENTS
We acknowledge support from NASA grant NNX14AL56G.We thank Coleman Krawczyk for his help with data visu-alization and Dr. N.P. Ross for his insightful comments onthis paper.This paper includes data collected by the Kepler mis-sion. Funding for the Kepler mission is provided by theNASA Science Mission directorate. All of the data presentedin this paper were obtained from the Mikulski Archive forSpace Telescopes (MAST). STScI is operated by the Associ-ation of Universities for Research in Astronomy, Inc., underNASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Sciencevia grant NNX13AC07G and by other grants and contracts.
REFERENCES
Agol, E. & Dexter, J. 2011 ApJL, 727, L24Andrae, R., Kim, D.-W., & Bailer-Jones, C.A.L. 2013,A&A, 554, A137Barth, A.J. et al. 2011, ApJ, 732, 121Borucki, W.J. et al. 2010, Science, 327, 5968, 977Box, G.E.P., Jenkins, G.M., & Reinsel, G.C. 2006,
TimeSeries Analysis: Forecasting and Control
John Wiley &Sons, HobokenBrockwell, P.J. & Davis, R.A. 2002,
Introduction to TimeSeries and Forecasting
ACS CCD Gains, Full Well Depths,and Linearity up to and Beyond Saturation
InstrumentScience Report ACSGoosmann, R.W., et al. 2006, A&A, 454, 741Graharm, M.J., et al. 2014, MNRAS, 439, 1, 703Hamilton, J.D. 1994,
Time Series Analysis
Princeton Uni-versity Press, PrincetonHao, L., Strauss, M.A., et al. 2005, AJ, 129,1783Healey, S.E. et al. 2007, ApJS, 171, 61Healey, S.E. et al. 2008, ApJS, 175, 97Ivezi´c, ˘Z., Connolly, A.J., VanderPlas, J.T., & Gray, A.2014,
Statistics, Data Mining, and Machine Learning inAstronomy: A Practical Python Guide for the Analysis ofSurvey Data
Princeton University Press, PrincetonKelly, B.C., Bechtold, J., & Siemiginowska, A. 2009, ApJ,701, 508Kelly, B.C., Treu, T., Malkan, M., Pancoast, A., & Woo,J.-H. 2013 ApJ, 779, 2, 187Kelly, B.C., Becker, A.C., Sobolewska, M., Siemiginowska,A., & Uttley, P. 2014, ApJ, 788, 33Fraquelli D. & Thompson, S.E. 2014,
Kepler Archive Man-ual (KDMC-10008-005)
Christiansen, J.L. et al. 2013,
Kepler Data CharacteristicsHandbook (KSCI-19040-004)
Van Cleve, J.E. & Caldwell, D.A. 2009,
Kepler InstrumentHandbook (KSCI-19033)
Fanelli, M.N. et al. 2011,
Data Processing Handbook (KSCI-19081-001)
Kinemuchi, K., Barclay, T., Fanelli, M., Pepper, J., Still,M., & Howell, S.B. 2012, PASP, 124, 963King, A.R., Pringle, J.E., West, R.G., & Livio, M. 2004,MNRAS, 348, 111Kolodziejczak, J.J. et al. 2010, Proc. SPIE, 7742, 77421GKozlowski, S., Kochanek, C.S., et al. 2010, AJ, 708, 927Krolik, J.H. 1999,
Active Galactic Nuclei: From The Cen-tral Black Hole to the Galactic Environment
PrincetonUniversity Press, PrincetonLiu, F.K., Zhang, Y.H. 2002, A&A 381 757Lyubarskii, Y.E. 1997, MNRAS, 292, 679MacLeod, C.L. et al. 2010, ApJ, 721, 1014MacLeod, C.L. et al. 2011, ApJ, 728, 26MacLeod, C.L. et al. 2012, ApJ, 753, 106Masetti, N. et al. 2006, A&A, 448, 547McHardy, I.M., Papadakis, I.E., Uttley,P, Page, M.J., Ma-son, K.O. 2004, MNRAS, 348, 783Mushotzky, R.F., Done, C., Pounds, K.A. 1993, ARA&A,31, 717Mushotzky, R.F., Edelson, R., Baumgartner, W., &Gandhi, P. 2011, ApJL, 743, L12Peterson, B.M. 2003,
An Introduction to Active GalacticNuclei
Cambridge University Press, CambridgePowell, M.J.D. 1994,
Advances in Optimization and Nu-merical Analysis
Kluwer Academic: DordrechtPrado, R. & West, M. 2010,
Time Series: Modeling, Com-putation, and Inference
CRC Press, Boca RatonRevalski, M. et al. 2014, ApJ, 785, 1, 60Richards, G.T., et al. 2006, AJ, 131, 2766 c (cid:13) , 1–20 Vishal P. Kasliwal, Michael S. Vogeley, and Gordon T. Richards
Ruan, J.J. et al. 2012, ApJ, 760, 1, 51Rutman, J. 1978, Proc. IEEE, 66, 9, 1048Rybicki, G.B. & Press, W.H. 1992, AJ, 398, 169Rytov, S.M., Kravtsov, Y.A., & Tatarskii, V.I. 1987,
Prin-ciples of Statistical Radiophysics
Springer-Verlag, BerlinScargle, J.D. 1981, ApJS, 45, 1Schmidt, K.B., Marshall, P.J., Rix, H.-W., et al. 2010, AJ,714, 1194Schmidt, K.B. et al. 2012, ApJ, 744, 147Sesar, B. et al. 2007, AJ, 134, 2236Shakura, N.I. 1973, SvA, 16, 756Shakura, N.I. & Sunyaev, R.A. 1973, A&A, 24, 337Simonetti, J.H., Cordes, J.M., & Heeschen, D.S. 1985, ApJ,296, 46Smith, J.A., et al. 2012, PASP, 124, 919, 1000Sobolewska, M.A., Siemiginowska, A., Kelly, B.C., & Nale-wajko, K. 2014 ApJ, 786, 143Stetson, P.B. 1987, PASP, 99, 191Stumpe, M.C. et al. 2012 PASP, 124, 919, 985Timmer, J. & K¨onig, M. 1995, A&A, 300, 707V´eron-Cetty, M.-P. & V´eron, P. 2010, A&A, 518, A10Wagner, S.J. & Witzel, A. 1995, ARA&A, 33, 163Walsh, D., Beckers, J.M., Carswell, R.F., & Weymann, R.J.1984, MNRAS, 211, 105Wehrle, A.E. et al. 2013, ApJ, 773, 2, 89Woodward, W.A., gray, H.L., & Elliot, A.C. 2012,
AppliedTime Series Analysis
CRC Press, Boca RatonZu, Y., Kochanek, C.S., Kozlowski, S., & Udalski, A. 2013,ApJ, 765, 106 c (cid:13)000