Structure Function and Variability Mechanism of Quasars from SDSS Stripe 82
aa r X i v : . [ a s t r o - ph . C O ] A ug Draft version
Preprint typeset using L A TEX style emulateapj v. 6/22/04
STRUCTURE FUNCTION AND VARIABILITY MECHANISM OF QUASARS FROM SDSS STRIPE 82
Alexey Voevodkin
Draft version
ABSTRACTTheoretical predictions for the ensemble quasar structure function are tested using multi-epochobservations of Stripe 82 collected by the Sloan Digital Sky Survey. We reanalyze the entire availablevolume of the g -band imaging data using difference image photometry and build high quality lightcurves for 7562 spectroscopically confirmed quasars. Our structure function includes ∼ . × pairsof measurements and covers a wide range of time lags between 3 days and 6.9 years in the quasarrest frame. A broken power-law fit to this data shows the presence of two slopes α = 0 .
33 and α = 0 .
79 with the break at ∼
42 days. The structure function compiled using only flux increases isslightly lower than that for variations of the opposite sign, revealing a slight asymmetry between theleading and trailing edge of a typical flare. The reality of these features is confirmed with monte-carlosimulations. We give simple interpretation of the results in the frames of existing theoretical models.
Subject headings: quasars
1. INTRODUCTION
Several decades after the discovery of quasars, thephysical origin of their variability is still not fully un-derstood. An unambiguous explanation of the intrinsicvariations in quasar emission will likely provide the keyto understanding the energy supply of quasars. The mostfrequently invoked model connects variability of quasarsto instabilities arising in the accretion disk around asuper-massive black hole (e.g. Kawaguchi et al. 1998;Pereyra et al. 2006). In the starburst model quasarlight curves are a superposition of very frequent super-nova explosions in the host galaxy (e.g. Terlevich et al.1992; Cid Fernandes et al. 1996; Aretxaga et al. 1997).Gravitational microlensing by compact bodies withinthe host—such as normal stars or certain candidatesfor dark matter—also produces observable flux changes(e.g. Hawkins 2002, 2007), but their contribution to theobserved level of variability remains uncertain. The mod-els can be distinguished based on their predictions for thedependence of the slope and the amplitude of the struc-ture function on the time-scale of variability in the quasarrest frame. Recently, individual quasar light curves havebeen successfully modeled as a damped random walk(DRW) process (Kelly et al. 2009, Kozlowski et al. 2010,McLeod et al. 2010). While the method allows one tostudy the structure function of individual objects as afunction of black hole mass, quasar luminosity, wave-length of observation (McLeod et al. 2010), it also fixesthe slope of the power spectrum in the high frequencyregime at − Space Research Institute (IKI), Profsoyuznaya 84/32, Moscow,Russia, 117997 Los Alamos National Laboratory, MS-D466, Los Alamos, NM,87545 narios. The most recent papers favor the disk in-stability model (Bauer et al. 2009, de Vries et al. 2005,Vanden Berk et al. 2004) with the exception of Hawkins(2002) who find evidence in support of the microlensingmodel. Here, we re-analyze the repeated g -band imag-ing of Stripe 82 collected by the SDSS using differenceimage photometry and build an ensemble structure func-tion over a wide range of time-scales in order to test thetheoretical predictions. Our results show that in theirpresent form none of the considered models explains thebehavior of the structure function on all time-scales. In-stead, a hybrid model including both a starburst and diskinstabilities naturally follows from the data.The paper is organized as follows. In Section 2 wedescribe the data and a method of extracting light curves.Basic properties of the structure function are given inSection 3. In Section 4 we build the structure functionfor the ensemble of quasars. We discuss the results inSection 5 and a summary is given in Section 6.
2. DATA ANALYSIS
Multi-epoch SDSS Imaging in Stripe 82
For the purpose of this study we reanalyze all available g -band drift-scan images of the SDSS Stripe 82 in the fi-nal SDSS data release, hereafter DR7 (Abazajian et al.2009). Quasar light curves were obtained using differenceimage photometry (Alard & Lupton 1998; Alard 2000)and a slightly modified version of the photometric soft-ware employed by the OGLE project (Wo´zniak 2000).The original motivation for difference imaging approachwas to facilitate a sensitive search for strongly lensedquasars (Kochanek et al. 2006; Lacki et al. 2009). Butthe application of the method also returns high qualityrelative light curves of all sources in a given data set,which we use to study the variability of known quasars.The full database of Stripe 82 contains 303 runs cov-ering nearly 300 deg of sky along the celestial equator( − ◦ ≤ RA ≤ ◦ and − . ◦ ≤ Decl ≤ . ◦ ). Thearea was scanned repeatedly during 10 years with someregions observed more than 80 times (Abazajian et al.2009). The SDSS imaging system consists of 30 CCD A. Voevodkin MJD (cid:0) (cid:1) (cid:2) (cid:3) R A , d e g Fig. 1.—
Spatial and temporal coverage of the SDSS Stripe 82imaging database. North runs are shown in blue and South runsin red. Note a much better time sampling in the later half of thesurvey. detectors organized in 6 raws (camcols), where a sin-gle CCD in each row is taking data in one of the 5filters ugriz . Since there are gaps between camcols,the width of the major stripe is split into North andSouth strips covered by independent scans. The SDSSdatabase of Stripe 82 contains 162 North and 141 Southruns. The spatial and temporal coverage of all runs isshown in Figure 1, where the blue and red lines corre-spond to the North and South runs. Each run is dividedinto 2048 × ×
128 pix margin on both ends (Stoughton et al.2002). The locations of the beginning and the end ofeach scan vary from one observation to another. Thisarrangement is not compatible with subtracting multipleimages of the same area, so the runs must be “glued”back together and then subdivided using consistent fieldboundaries. The rebinning starts from RA min , the firstlocation covered by at least 20 runs (determined sepa-rately for North and South runs). The rearranged framesof 2048 × max , beyond whichfewer than 20 runs are available.Some of the data in Stripe 82 were obtained in non-photometric conditions and in relatively poor seeing(DR7). The performance of the difference image anal-ysis method is sensitive to the quality of the images, andtherefore our next step is to reject images with poor see-ing, high background, and significant background gradi-ents. A small number of images taken toward the ends ofthe run have strong background gradients due to the Sunor Moon rising or setting. Occasionally, a problem withthe bias level subtraction caused a jump in the countlevel between the left and right parts of the image corre-sponding to individual amplifiers on the CCD chip. Suchimages were removed from the analysis. A good qualityreference image was prepared for each field using 10 im-ages with the best overall seeing, background, and trans-parency. This is done with the help of the SExtractorprogram (Bertin & Arnouts 1996) by requiring a median F W HM ≤ .
5” for point sources (CLASS STAR ≥ . Difference Image Photometry of SDSS Quasars
Interpolation of images to the same pixel grid was ac-complished using approximately 50 stars identified on allimages of the field. The time base line of the Stripe82 data set approaches a decade. We found that a sig-nificant fraction of stars shows a detectable proper mo-tion producing time-variable residuals in the shape of adipole in the difference images (cf. Lacki et al. 2009).In order to limit the impact of such problems on thequality of the final photometry, we selected the referencestars with the proper motion µ ≤
100 mas yr − from theBramich et al. (2008) catalog in the overlap area (about50%). The spatial variability of the PSF matching ker-nel was computed using a third-order polynomial fit tomodel coefficients and a second-order polynomial for thedifferential background. Again, slow stars from the cat-alog of Bramich et al. (2008) were preferentially used tocompute the solution. The DIA light curve is computedby adding difference fluxes measured with the DIA pack-age to the flux in the reference image. Object magni-tudes were converted to standard SDSS g magnitudes byadding a median offset between stars identified on thereference image and corresponding objects in the SDSScatalog.The photometric quality is demonstrated in Figure 2showing the magnitude scatter versus baseline magnitudefor a typical field. Here, the baseline is the median mag-nitude and the variability of the source is calculated asthe 68-th percentile of absolute deviations from the me-dian. Black points correspond to the sources identifiedin that field and blue line shows median of the variabilityof the sources from all fields.In most fields the frame-to-frame accuracy reaches 7mmag for bright unsaturated sources. For a significantfraction of stars the photometric scatter is worse thanthat and limited by the residuals due to proper motions.The photometric uncertainty is 1% or better for sourceswith m g ≤
18 and around 10% for sources with m g ≈ m g >
22. Those were excludedfrom the structure function analysis Section 4) due tolarge uncertainties in their light curves. In the next step,we match spectroscopically confirmed quasars from DR7Quasar Catalog (Schneider et al. 2010) with our photo-metric database. There are 9519 DR7 catalog quasarswithin the area covered by Stripe 82, and we can buildDIA light curves for 7614 of them. There are several rea-sons why we miss a significant number of quasars. Afterexclusion of poor quality images many fields are left withfewer than 20 epochs and are removed from further anal-ysis. An average 10% of the image area is lost due tomasking of very bright stars. Some fields are rejected forlack of stars that can be used to calibrate light curves.Finally, our pipeline occasionally fails on a small fractionof images with problems in fpC frame headers.An example quasar light curve resulting from this pro-cedure is shown in the Figure 3. There are 43 observa-tions in the light curve and time counts from the firstepoch in Stripe 82 (run 94, fpC frame 12). The numberof epochs per quasar in our sample varies between 20 andSO Structure Function and Variability Mechanism 3
16 18 20 22 24 m g -2 -1 (cid:4) Fig. 2.—
Variability diagram in the g band for a typical field inSDSS Stripe 82. The magnitude scatter of all measurements for agiven source is plotted against the baseline magnitude (median).The photometry is based on our reanalysis of the SDSS scans us-ing difference image photometry. Black points correspond to thesources identified by SExtaractor. Sources above the magenta lineare considered to be variable. Blue line shows median variabilityof the sources from all fields. T, years m g Fig. 3.—
Example difference imaging light curve for QSOSDSS 223154 .
15 + 004747 . z = 0 .
56 with the mean of 36.
3. STRUCTURE FUNCTION
The time sampling of our quasar light curves is rela-tively sparse (see Figure 3) and insufficient for the use ofspectral methods such as the frequency-power spectrumanalysis. Following the common practice in AGN vari-ability studies, we focus on the information contained inthe structure function. A number of structure functiondefinitions are in use differing in relatively minor detailsthat are not important for this study. We adopt a struc-ture function (hereafter SF) in the form: S ( τ ) = N ( τ ) X i 4. THE ENSEMBLE STRUCTURE FUNCTION OFQUASARS The final sample contains 7562 light curves. All pos-sible pairs of measurements within each light curve arecollected for a total of 4,796,468 magnitude differences.After reducing the time lag to the quasar rest frame us-ing spectroscopic redshifts from the SDSS Quasar Cata-log (Schneider et al. 2010), the data are binned in time.The bin size was adjusted to reach 20,000 measurementsper bin. The resulting ensemble SF consists of 239 binsand is shown in Figure 4 (blue points). Error barswere estimated by dividing each timescale bin into 20smaller bins, evaluating the mean SF in each, and cal-culating the r.m.s. scatter of the resulting values (seealso de Vries et al. (2005)).In order to verify that we detect quasar variability andguard against the possibility of contamination by unre-lated factors, we repeated this calculation for constantstars. For that purpose we select stars below the ma-genta line in Figure 2 with m g ≤ 22. The result is shown Structure functions of individual sources can be effectivelyused for identification of quasars (see e.g. Schmidt et al. 2010;Palanque-Delabrouille et al. 2011). A. Voevodkin (cid:5) lg( (cid:6)(cid:7) , days) (cid:8) (cid:9) (cid:10) (cid:11) (cid:12) l g ( S F ) lg( (cid:13)(cid:14) , days) (cid:15) (cid:16) (cid:17) (cid:18) l g ( S F ) Fig. 4.— Top panel: Raw ensemble structure function of 7562quasars from SDSS Stripe 82 (blue) compared to the structurefunction for ∼ Bottom panel: Noise-corrected structure functionfor an ensemble of quasars (blue) and the best fit broken power-lawmodel (red). The oscillations for time-scales above ∼ 100 days arean artifact of the time sampling. The model offers only an approx-imate description of the structure function over the full extent ofthe data. in the top panel of Figure 4 (green points). As expected,the SF of non-variable stars is flat over the entire rangeof time-scales. We conclude that the observed variabilityof quasars is intrinsic and not caused by unrecognizedmeasurement errors. The stellar structure function ex-ceeds the noise level for the quasar sample because thestars selected for this comparison are fainter on averagethan the quasars.The quasar SF in the top panel of Figure 4 displaysseveral slopes and breaks. It starts flat at τ ≤ ≤ τ ≤ 10 days gradually transitions to an ap-proximate power-law increase over 10 ≤ τ ≤ 40 days,and finally breaks to a shallower slope above τ ≃ m g lg( (cid:19)(cid:20) , days) (cid:21) (cid:22) (cid:23) (cid:24) l g ( S F ) Fig. 5.— Noise-corrected SF A (blue) and SF B (green). Blackpoints show the g -band structure function from Wilhite et al.(2008) who used the SF B estimator. Note a close agreement inboth the shape and amplitude. reaching the level of noise at different τ (see Equation 2).The slopes seen at 10 ≤ τ ≤ 40 days and τ ≥ 100 have aphysical origin and are discussed further below. The timeinterval between 40 and 100 days is poorly sampled anddifficult to interpret. This is a result of the SDSS scan-ning strategy for Stripe 82 that produces relatively fewobservations separated by 90–270 days in the observerframe (cf. Figure 3) combined with the mean redshift ofthe sample ∼ . S ( τ ) = (cid:26) β ( τ /τ ) α if τ ≥ τ β ( τ /τ ) α if τ ≤ τ . (3)We find α = 0 . α = 0 . τ = 42 . β = 0 . 13 (solid red line). The reduced χ of the fitis far from unity, i.e. Equation 3 does not provide acomplete description of the data. However, the maingoal here is to detect and estimate the two slopes in theoverall shape of the SF for an ensemble of quasars, andcompare those coefficients with theoretical predictions.We also estimate errors on the derived parameters using’jackknife’ resampling. Namely, we randomly select 1/3part of the sample, drop it, and fit ensemble SF builtfrom the rest. We repeat the procedure 100 times andestimate the error on the given parameter as r.m.s. ofproper distribution. Such estimation gives very smallerrors (0.01, 0.02, 3.9 days, and 0.01 for α , α , τ , and β correspondingly) saying that the shape of ensemble SFis quite stable to small variations of the sample.For comparison we plot the g -band SF fromWilhite et al. (2008) (black points in the bottom panel ofFigure 4). The overall agreement in shape between theSFs is very good and the discrepancy in normalization isexplained by the different SF type used in Wilhite et al.SO Structure Function and Variability Mechanism 5 m a g n i t u d e m a g n i t u d e Fig. 6.— Light curve simulations. A random simulated lightcurve (red) is compared with a real light curve of one of the quasarsin the sample (blue). The simulated light curves were generatedfrom a power spectrum in the form of a single power-law with theslope − . 66 (top) and − . 58 (bottom) corresponding to structurefunctions in Figure 7. (2008). There are two types of SF commonly used in lit-erature. Following Bauer et al. (2009) we refer to themas SF A and SF B : S A ( τ ) = p h (∆ m ) i (4) S B ( τ ) = r π h| ∆ m |i . (5)The structure function in the form SF B was introducedby di Clemente et al. 1996. Our noise-corrected SFs andthe g -band SF B from Wilhite et al. (2008) are shown inFigure 5. There is a nearly perfect agreement in bothshape and normalization between the B-type SFs derivedfrom the two data sets.The interpretation of results obtained usingSF as a statistical tool is not straightforward.Emmanoulopoulos et al. (2010) show that spuriousbreaks can appear in the SF due to sparse time samplingand windowing effects. Moreover, the slopes derivedfrom SF modeling are typically much more uncertain (cid:25) lg( (cid:26)(cid:27) , days) (cid:28) (cid:29) (cid:30) (cid:31) l g ( S F ) Fig. 7.— Confirming the presence of two slopes in the quasarstructure function. The simulated light curves were generated froma power spectrum in the form of a single power-law with the slope − . 66 (blue) and − . 58 (green) corresponding to the upper andlower portions of the best fit model (red line) from Figure 4. than suggested by naive error estimates. Does it meanthat the presence of two slopes in Figure 4 is an artifact?To answer this question and validate our results weperform a set of monte-carlo simulations. The goal ofthis calculation is to verify that for a set of light curvesgenerated from a single power-law spectrum with a givenvalue of α and having the time sampling of real quasars,the algorithm of computing SF does not introduce anysignificant breaks or biases in the estimated slope.The simulation is based on the algorithm ofTimmer & Koenig (1995) and reproduces both the cor-rect power spectrum as well as the phase mixture. Thepower spectrum in the form P ( f ) ∝ f − (2 α +1) corre-sponds to SF S ( τ ) ∝ τ α . We simulate light curvesaccording to a prescription in Emmanoulopoulos et al.(2010). The simulated time series are then normalizedto have the same r.m.s. and mean as the real light curves.Examples are shown in the Figure 6.We simulate two sets of light curves with the power-lawspectral density ( − . 66 and − . α = 0 . 33 and α = 0 . 79. Therefore,the large observed difference in the SF slope between theshort and long timescales is not due to the time sampling,but rather of a physical origin. Possible interpretationsare discussed in the next section.The ensemble SF considered here is built from quasarswith different intrinsic characteristics such as black holemass or luminosity. Therefore, it is interesting to checkwhether two slopes in SF are characteristic for the holerange of, say, black hole masses or have dependence on A. Voevodkin TABLE 1Best fit parameters Subsample α α τ , days β low mass 0 . ± . 02 0 . ± . 04 55 . ± . . ± . . ± . 01 0 . ± . 03 57 . ± . . ± . . ± . 02 0 . ± . 09 36 . ± . . ± . mass. Using black hole mass estimations from Shen et al.(2008), we select three subsamples of quasars with blackhole masses in ranges ≤ × M ⊙ (1953 sources, lowmass range), 5 × ÷ . × M ⊙ (2108 sources, meanmass range), ≥ . × M ⊙ (1276 sources, high massrange) and build SF for each range of masses. All threeSFs show existence of two slopes (see Figure 8, toppanel). We fit SFs by broken power low model (Equa-tion 3). The best fit parameters and formal estimationof errors by ’jackknife’ resampling are given in the Ta-ble 1. In order to have better understanding of the uncer-tainties in the parameters we also plot regions showingminimal and maximal values of the fits for all three SFs(see Figure 8, bottom panel). As it is seen from the Fig-ure 8, the SFs slopes tend to increase with the growth ofthe mean subsample mass but still have common values.More subtle study is needed in order to show that thetrend in slopes as a function of black hole mass is real.Here we only conclude that the presence of two slopes inthe ensemble SF is independent on the considered blackhole mass range. 5. DISCUSSION We consider three processes frequently invoked to ex-plain the variability of quasars: instabilities in the ac-cretion flow around a supermassive black hole, super-nova explosions related to the starburst phenomenon,and gravitational microlensing by compact bodies in thehost galaxy (Kawaguchi et al. 1998; Hawkins 2002). SFsof light curves predicted by each model are character-ized by different slopes. The expected SF slopes are0 . ± . 03, 0 . ± . 03, 0 . ± . 08 respectively for vari-ability generated by miscrolensing, disk instability, andstarburst mechanisms (Hawkins 2002). The SF slopesfound in Section 4 are 0.79 and 0.33 correspondingly fortime lags below and above 42 days. Taken at face value,the first slope is consistent with the starburst model, andthe second slope falls half way between the predictionsof the disk instability model and the microlensing model.Note that SF slopes in the range 0.30–0.35 have been ob-tained in other works, but the authors still attributed thevariability to disk instabilities. Bauer et al. (2009) andMeusinger et al. (2011) present an overview of recent re-sults.Asymmetries in the SF provide an additional test ofmodel predictions. The SF that only includes thosepairs of measurements for which the flux increases withtime, S + , may be different from S − , the SF character-izing a decrease in flux. Kawaguchi et al. (1998) showthat in starburst models S + > S − . Disk instabilitiesproduce S − > S + and microlensing is symmetric, i.e. S − = S + when averaged over sufficiently long time in-tervals (Hawkins 2002). The difference between S − and S + depends on the parameters of the model and there- lg( !" , days) $ % & l g ( S F ) ’ lg( () , days) * + , - . / l g ( S F ) high massmean masslow mass Fig. 8.— Top panel: Ensemble SF for quasars with black holemasses lying in low mass range (red points), mean mass range(green points ), and high mass range (blue points) (see text fordefinition of ranges). Thin black lines show best fits for each SF. Bottom panel: Regions showing the scatter of the best fits obtainedduring ’jackknife’ resampling. Red lines show the scatter of the lowmass range SFs best fits, green region and blue lines show the samefor mean mass and high mass ranges correspondingly. fore S − ≃ S + is still possible in both the starburst andthe disk instability model (Kawaguchi et al. 1998).Noise-corrected S + and S − functions are shown in Fig-ure 9 with blue and green lines respectively. For longtime lags in the range 300–1600 days, where the SF slopeis ∼ . 33, we have S − > S + , as predicted by the diskinstability model. The significance of the difference isnot very high, but allows us to exclude the microlensingmodel. The difference disappears for time-scales below100 days, where the slope is consistent with the starburstmodel. This may suggest that a typical supernova ratein a starbursts in the host galaxies of quasars is as highas ∼ 100 yr − and effectively washes out any detectableasymmetry (see Fig. 6 in Kawaguchi et al. (1998)). How-ever, the rate should be treated with caution since it isnot derived from direct measurements. Moreover, thepredictions for S + and S − were derived from simulationsof a single variability mechanism. The addition of an-SO Structure Function and Variability Mechanism 7 , days S F Fig. 9.— Testing the asymmetry in the quasar structure functionbetween the leading and trailing branch of a typical flare. Fluxincreases contribute to S + (blue) and flux drops contribute to S − (green). other variability process can wash out the difference be-tween S + and S − . This is an interesting topic for futureresearch. The difference in S + and S − was also investi-gated by de Vries et al. (2005) and Bauer et al. (2009).The former study found tentative evidence that S + > S − (light curves have fast rise and slow decay as in starburstmodel) for 400 days ≤ τ ≤ S + ≈ S − over the full range of measurements( τ ≤ 6. SUMMARY We applied the image differencing technique to the en-tire g -band imaging data set of SDSS Stripe 82 and con-structed high quality light curves for 7562 spectroscopi-cally identified quasars with no less then 20 epochs. Thevariability analysis was performed using the ensemblestructure function. While the shape and normalization ofthe quasar structure function derived here are in perfectagreement with the results of Wilhite et al. (2008) basedon an earlier SDSS data release, our structure function is estimated from a much larger data set and covers asubstantially wider range of time lags from 8 hours to6.9 years in the quasar rest frame.We find that the ensemble SF reveals two distinctpower-law slopes with the break around 42 days andconfirm the presence of these features with monte-carlosimulations. The presence of two slopes in the ensem-ble SF is independent on the black hole mass range ofquasars used to build it. The slopes estimated froma broken pawer-law fit to the data are α ≃ . 79 for τ ≤ 100 days and α ≃ . 33 for τ ≥ 300 days. Usingpredictions from the theoretical models, the variabilityof quasars on times scales τ ≤ 100 days can be explainedby a starburst model. The slope of the structure func-tion on time-scales longer than about 3 months is ratherclose to models where variability is explained by gravi-tational microlensing. However, on time-scales τ ≥