Observational nonstationarity of AGN variability: The only way to go is down!
Neven Caplar, Theodore Pena, Sean D. Johnson, Jenny E. Greene
DDraft version January 24, 2020
Typeset using L A TEX twocolumn style in AASTeX63
Observational Nonstationarity of AGN Variability: The Only Way to Go Is Down!
Neven Caplar , Theodore Pena , Sean D. Johnson , and Jenny E. Greene Department of Astrophysical Sciences, Princeton University, 4 Ivy Ln., Princeton, NJ 08544, USA Department of Physics and Astronomy, Tufts University, Medford, MA 02155, USA (Received -; Revised -; Accepted -)
Submitted to ApJLABSTRACTTo gain insights into long-term variability of Active Galactic Nuclei (AGN), we analyze an AGNsample from the Sloan Digital Sky Survey (SDSS) and compare their photometry with observationsfrom the Hyper Suprime-Cam survey (HSC) observed (cid:104) . (cid:105) years after SDSS. On average, the AGNare fainter in HSC than SDSS. We demonstrate that the difference is not due to subtle differences inthe SDSS versus HSC filters or photometry. The decrease in mean brightness is redshift dependent,consistent with expectations for a change that is a function of the rest-frame time separation betweenobservations. At a given redshift, the mean decrease in brightness is stronger for more luminousAGN and for objects with longer time separation between measurements. We demonstrate that thedependence on redshift and luminosity of measured mean brightness decrease is consistent with simplemodels of Eddington ratio variability in AGN on long (Myr, Gyr) timescales. We show how our resultscan be used to constrain the variability and demographic properties of AGN populations. Keywords: accretion, accretion disks - black hole physics - methods: data analysis - quasars: general INTRODUCTIONChanging flux levels with time are nearly ubiquitousamong Active Galactic Nuclei (AGN). Studies of theseluminosity fluctuations, i.e., AGN variability, have en-abled measurements of central supermassive black holemasses (e.g., Bentz 2015), added insights on the struc-ture of AGN accretion disks (e.g., Fausnaugh et al.2016), and provided powerful AGN selection techniques(e.g., Schmidt et al. 2010). AGN variability has beendirectly observed in large samples on timescales rangingfrom minutes to days, years, and decades (e.g., MacLeodet al. 2010, 2012; Morganson et al. 2014; Cartier et al.2015; Caplar et al. 2017; Smith et al. 2018). Throughindirect methods and simulations, AGN variability hasalso been studied on Myr and Gyr scales (e.g., Novaket al. 2011; Bland-Hawthorn et al. 2013; Sartori et al.2018).Most of the direct observational studies mentionedabove quantify AGN variability as a weakly station-
Corresponding author: Neven [email protected] ary process, i.e., with the assumption that the meanluminosity of statistically large ensembles of AGN doesnot change with time. Empirically, stochastic variabil-ity measured in these short-term studies dominated anypossible subtle changes in the mean brightness occur-ring during the duration of the studies. From the-oretical grounds, the stochastic variability is thoughtto be reflective of the details of the physics of AGNaccretion disks and other nearby structures, while themean change of luminosity would be connected to longtimescale accretion processes thought to have minimalimpact on typical survey timescales (but see Lawrence2018).The assumption of no change of mean brightness onshort timescales differs from the long-term studies ofAGN activity, which indicate large changes in AGN ac-tivity on Myr and Gyr scales. The firmest observa-tional proof comes from the studies of individual ex-tended AGN photoionized clouds, so-called “Voorwerp”objects (e.g., Sartori et al. 2016; Johnson et al. 2018),and the He II transverse proximity effect (e.g., Schmidtet al. 2018), which clearly show that some AGN exhibit a r X i v : . [ a s t r o - ph . GA ] J a n Caplar et al. order-of-magnitude changes in their luminosity on 10 -10 yr time-scales.There has been comparatively little observational re-search on the deviations from the symmetric behavior ofAGN variability. Numerous early studies conducted ob-servationally difficult searches for potential differencesin the variability properties of AGN that were becom-ing brighter or dimmer, but found little or no evidencefor statistical differences in the properties of AGN lightcurves that were fading or getting brighter (de Vrieset al. 2003, 2005; Bauer et al. 2009; Voevodkin 2011).MacLeod et al. (2012) combined the earliest statisticallysignificant sample of AGN measurements from the Palo-mar Observatory Sky Surveys (POSS) with the SDSSdata. They noted that objects from POSS are dimmerwhen observed in SDSS. They concluded that this maybe explained by a “Malmquist-like” bias, i.e., the factthat luminosity-selected sample of variable objects willnecessarily be dimmer in the later survey, even if thereis no change in the mean brightness of the underlyingsample. A similar conclusion was reached by Rumbaughet al. (2018) who studied examples of extreme variabil-ity by comparing SDSS and Dark Energy Survey mea-surements. Morganson et al. (2014) also found the de-crease of the mean brightness on decade timescales, forthe sample of AGN from the SDSS observed in Pan-STARRS1, but attributed this effect to the filter differ-ences.Here, we use the AGN sample from SDSS and measuretheir mean brightness in SDSS and HSC. The depth,size, and time separation from SDSS, and the quality ofthe HSC survey make it especially suitable for this kindof study. In this work, we aim to show that AGN exhibitchanges in their mean brightness in a redshift and lumi-nosity dependent manner on the timescales accessiblewith past (SDSS) and current (HSC) surveys.The code and the data needed to reproduce all ofthe results mentioned in this work are available atgithub.com/nevencaplar/AGN-Going-Down. OBSERVATIONS2.1.
Data
To study AGN variability on decade timescales, weidentified AGN from the SDSS (York et al. 2000) DR7Quasar catalog (Schneider et al. 2010) that were also ob-served later by the Hyper Suprime-Cam (HSC; Miyazakiet al. 2018) Subaru Strategic imaging survey. The SDSSsurvey used a dedicated 2.5 m (Gunn et al. 2006) tele-scope at Apache Point Observatory to obtain imagesin five optical bands (ugriz) over a large patch ( ∼ ) of the northern sky. For a presentation of the pho-tometric calibration and selection function of objects, we refer the reader to the detailed discussion in Schneideret al. (2010). The HSC survey is a wide-field opticalimaging program being conducted with the 8.2 m Sub-aru telescope. The second public data release (madeavailable in May 2019; Aihara et al. 2019) covers around300 deg overlapping with the SDSS footprint. The HSCdata in five optical bands (grizy) are sensitive down to ≈ g , r and i psf-magnitudes. We excluded all of the objects with anyflags showing problems in the calibration. This conser-vative cut ensures that our conclusions are not drivenby possible problems in the brightness measurements inthe HSC pipeline. This procedure yields 5919 matchedAGN found in both surveys.2.2. Main result
To measure the mean difference in the brightness be-tween the two surveys, we split the sample in bins ofredshift, each consisting of 100 objects. This numberenables us to follow the redshift evolution of the trendsin some detail, while minimizing the statistical uncer-tainty in the mean brightness change. For each redshiftbin, we then measured the mean and the median differ-ence between the observed psf-magnitudes in the SDSSand HSC surveys. We also verified that our conclusionsare unchanged when using fixed 3 arcsec aperture mag-nitudes.The resulting mean change in flux is plotted as a func-tion of redshift in Figure 1. To avoid cluttering the plot,we only show data points at each redshift for measure-ments in the g-band. We estimated uncertainties on themean value at each redshift by bootstrapping the un-derlying 100 AGN in each bin. We choose to presentg-band variability given that contribution of the host-galaxy light, however small for these bright AGN, willbe smallest in the bluest available band. However, re-sults for all three bands are very similar. We also showlinear fits to the data in all three bands, where one canexplicitly see the similarity between all of the results.We have also verified that the redshift evolution effectis present if we use median differences instead of meandifferences of magnitudes, but the magnitude of the ef-fect is somewhat decreased. For instance, the best fit forthe median difference is − .
139 + 0 . z , while for themean difference it is − .
176 + 0 . z . The fact that theeffect is still present when using the median shows thatit cannot be fully explained by a relatively small num-ber of extremely variable quasars (e.g., MacLeod et al.2016; Rumbaugh et al. 2018). We also show a linearfit to the data as a function of rest-frame time separa- bservational Nonstationarity of AGN variability m e a n p s f - m a g . ( S D SS - H S C ) f a i n t e r i n H S C redshift fit to AGN data, g-bandtime fit to AGN data, g-bandredshift fit to AGN data, r-bandredshift fit to AGN data, i-bandAGN data, g-band t rf > [restframe years] Figure 1.
Mean difference in the measured psf-magnitudes for the sample of AGN from SDSS that have been observed in bothSDSS and HSC. Blue points show the data for the g-band, while the blue dashed line and the shaded region show the linearfit as a function of redshift to the data and 1- σ uncertainty band. The red and the black dashed lines show the linear fit as afunction of redshift in r- and i- bands, respectively. We do not show the data and uncertainty bands for r- and i-band to improvethe clarity of the figure, but these are comparable to the g-band quantities. The dotted line shows a fit to the g-band data as afunction of restframe time-separation between measurements, indicated on the upper axis. tion between two measurements, i.e., as a function of14 .
85 years / (1 + z ), where 14.85 years is the mean timeseparation between observations (see Section 2.5). Thisfit also provides a good explanation for the observeddata. We discuss the proposed model in which mea-sured changes of the mean/median flux are the conse-quence of the long-term AGN behaviour and primarilydepend on the rest-frame time separation between thetwo measurements further in Sections 2.4, 2.5 and 3.To ensure that the observed redshift dependence isnot a spurious artifact due to differences between thetwo surveys, we conduct four different checks that welist here: • consideration of filter differences; • constructing a control sample; • separating the AGN sample according to bright-ness; and • separating the AGN sample according to the timeseparation between the SDSS and HSC observa-tions.We elaborate on each of these procedures in some de-tail below. For consistency, we always show the meandifference in the g-band and conduct linear fits as a func-tion of redshift, but all of our conclusions are applicable to all three bands and fitting variables (redshift or rest-frame time separation) .2.3. Filter difference and control sample
We performed two experiments to assess the poten-tial impact of differences in the photometry between theSDSS and HSC surveys that could lead to spurious ap-parent change of measured brightness. First, to assessthe potential impact of differences in the filter systemson the measured magnitudes, we predicted g SDSS − g HSC for the mean SDSS quasar spectral energy distribution(Vanden Berk et al. 2001) as a function of redshift usingthe SDSS (Fukugita et al. 1996) and HSC (Kawanomotoet al. 2018) defined system throughput including the fil-ters, telescopes, cameras, and the survey standard at-mospheres. Second, we constructed a control sampleconsisting of nonvariable stars with colors similar to theAGN. The stars were taken from the catalog of nonva-riable objects from the equatorial Stripe 82 presentedin Ivezi´c et al. (2007) which we additionally cleanedby removing suspected AGN from Flesch (2015). Foreach AGN we find the star (repetition allowed) thatminimizes the Euclidean distance between the measuredmagnitudes in the g-, r- and i- bands from SDSS. Afterthat, we treated the resulting catalog of stars in exactly Figures showing results for all of the possible combinations ofchoices for the used observed bands, fitting variables, and usingmean/median to derive results can be created from the code andthe data available in the GitHub repository
Caplar et al. m e a n g - p s f m a g . ( S D SS - H S C ) f a i n t e r i n H S C fit to AGN data, g-bandfit to control sample data, g-bandexpected difference due to SDSS/HSC filtersAGN data, g-bandcontrol sample data, g-band Figure 2.
Mean difference in the measured psf-magnitudes for the sample of AGN from SDSS, and the control sample of starsthat match these AGN in color. The blue points show the data for AGN, while the blue line and the shaded region show thelinear fit to the data and 1- σ uncertainty. This is equivalent to the data and fit shown in Figure 1. The maroon points, line,and shaded region show equivalent quantities constructed for the sample of nonvariable stars in Stripe 82 region. The black lineshows the expected redshift dependence due to filter differences between the two surveys. the same way as we have treated the AGN sample. As,by definition, we expect no change in the brightness ofthese stars when imaged in the two surveys, any sys-tematic differences between the two surveys will be ex-pressed in this comparison. These experiments captureeffects both from filter differences and from any differ-ences in the PSF magnitude measurement techniques.We show the results of this experiment and deduced ef-fects of filter differences in Figure 2. The overall decreasein mean flux is not present in the control sample of non-variable stars. In particular, we emphasize the absenceof “redshift” trend in the control sample. This is an ex-pected result, as the different “redshifts” for the controlsample correspond to only relatively small changes inthe mean color of the objects, which does not affect thecalibration of the surveys greatly. We also note that theexpected filter differences between the two surveys pro-duce a relatively small and almost redshift-independenteffect for the AGN sample. This is due to the small dif-ferences between the SDSS (Fukugita et al. 1996) andHSC (Kawanomoto et al. 2018) g-bands, and character-istic power-law SED of an AGN that results in modest u − g and g − r colors of ≈ − . . With λ eff = 4770 ˚A and FWHM = 1379 ˚A for SDSS versus λ eff = 4754 ˚A and FWHM = 1395 ˚A for HSC Split according to brightness
We then continue to study the redshift effect aftersplitting our sample in brightness. We do this for twoseparate reasons. Observationally, we expect that sys-tematic differences between the surveys would be morestrongly manifested for objects that have lower bright-ness, as various errors and uncertainties start to domi-nate closer to the brightness limit of the SDSS survey.Also, among less luminous AGN, which occur predomi-nantly at lower redshifts, flux from the host galaxy maystart to be nonnegligible (Shen et al. 2011), which mightbias our results. Additionally, given that variability isenhanced at lower luminosities, “Eddington bias”, refer-ring to the fact that intrinsically lower-luminosity AGNmight get scattered into the selection of the first, shal-low, survey and then “return” to their mean value whenobserved later, would produce a measured mean changeof brightness that would be more noticeable at lowerluminosities.On the other hand, physically, we would expect thatthe mean brightness change would be larger for more lu-minous AGN. Under the assumption that all AGN arethe members of the same population, with the same un-derlying Eddington ratio distribution, AGN are brightenough to be detected in a shallow flux-limited surveyonly during rare parts of their life-cycle. We wouldexpect that, on average, the population of such AGNwould gravitate to their mean, low-flux state as a func-tion of time. In particular, if this assumption is correct,we would expect that brighter AGN are in the moreextreme part of their life-cycle, occupying more extremeends of their long-term Eddington ratio distribution. We bservational Nonstationarity of AGN variability redshift m e a n p s f - m a g . ( S D SS - H S C ) f a i n t e r i n H S C fit to the 0-20% AGN data, sorted by brightness (brightest)fit to the 20-40% AGN data, sorted by brightnessfit to the 40-60% AGN data, sorted by brightnessfit to the 60-80% AGN data, sorted by brightnessfit to the 80-100% AGN data, sorted by brightness (dimmest)0-20% AGN data, sorted by brightness (brightest)80-100% AGN data, sorted by brightness (dimmest) redshift -0.100.1 m e a n e rr Figure 3.
Mean difference in the measured psf-magnitudes in the g-band for the sample of AGN separated according to theirbrightness at each redshift. Different colored lines show the results of the linear fit to the subsets of the data. They have beenconstructed by separating the data at each redshift in quintiles, according to their brightness. Binning is coarser than in Figures1 and 2 as to preserve statistical power of individual data points. We do not show the data points for the three inner quintilesof the data to improve the clarity of the figure. Lower panel shows mean 1- σ uncertainty bands around linear fits. Note thatscaling on the y-axis is different than in the upper panel. would therefore expect that brighter AGN will decreasetheir brightness more during any given observed time-frame.We split the data in each redshift bin into five furtherbins, according to their observed brightness in SDSS. Wethen proceeded to fit the data in each of these brightnessbins with a linear function and show the results of thefitting procedure in Figure 3. As uncertainties on thefits are quite similar for all of the 5 bins, we show themean error on the fit in the separate panel below themain panel. We see that the effect is indeed stronger forthe brighter AGN, as we expected from our theoreticalreasoning. This makes us even more confident in thephysical nature of this effect. We also wish to point outthat the observed brightening for the dimmest objects ismostly driven by the last point at the highest redshift,and it is not obvious that it is also a physical result.2.5. Split according to the time separation
As a final check, we separated our sample in the quin-tiles according to the time separation between the obser-vations. As both surveys took data over several years,we split the samples into those taken, by random chance,at the shortest and longest time intervals and comparethe results. If the change of mean brightness is mostlydue to observational effects, we would expect no differ-ence between the short and long separation datasets, while if the difference is physical we would expect to seesome difference between these two sets.This experiment is somewhat complicated by the factthat we, at this stage, are only working with the stackedHSC data, i.e., the measured brightness of any object isa combination of measurements at different times dur-ing the duration of the survey. For HSC data, we takethe mean of all of the observation times that go intoeach stacked observation and use that “mean time” asthe time of the observation. The distribution of timedifferences between the surveys is roughly normal, withthe mean at 14.85 years. We then create a sample outof the data in each redshift bin for which the time sep-aration is within the shortest time separation quintile(short separation sample) and out of the data that arein the longest time separations quintile (long separationsample). Mean time separation for the short separationsample is 12.94 years, and for the long separation sampleit is 16.89 years.We then proceeded as before to study the redshift de-pendence of each of these samples. We show the data,the results of the linear fit to the full data and long/shortseparation datasets in Figure 4. We see that, in gen-eral, long separation data do indeed tend to show largerchanges between the two surveys. Of course, the resultsare quite noisy which is not surprising given the samplesizes and underlying stochastic variability. In Figure 4
Caplar et al. redshift m e a n p s f - m a g . ( S D SS - H S C ) f a i n t e r i n H S C fit to all AGN datafit to short separation AGN datafit to long separation AGN datapredicted fit to short sep. AGN data, derived from long sep. dataall AGN datashort separation AGN datalong separation AGN data redshift -0.100.1 m e a n e rr t rf for the whole sample > [restframe years] Figure 4.
Mean difference in the measured psf-magnitudes in the g-band for the sample of AGN split according to the timeseparation between the measurements. The blue points show the data for the whole sample of AGN in the g-band, while the blueline shows the linear fit to the data. This is equivalent to the data and fit shown in Figure 1, although the binning is different tomatch binning for the short and long separation samples (shown in pink and green, respectively). The black dotted line showsthe simplest “derivation” of the short separation fit, which has been calculated from the long separation fit by reducing it bythe ratio of the mean time-separations for these two samples (12.94, 16.89 years). The lower panel shows 1- σ uncertainty bandson these linear fits for the short and long separation data. Note that scaling on the y-axis is different than in the upper panel. we also show the expected linear fit for the short sepa-ration sample, which was derived from the long separa-tion sample by multiplying the slope with the ratio ofmean time separations of each sample, i.e., with the fac-tor 12.94/16.89. This is a simplified assumption, as themean change in brightness is not necessarily linear withtime, but we see that the modifications explain well themagnitude of the observed difference. MODELING AND DISCUSSIONIn this section we discuss how the effect can be used toconstrain parameters of AGN accretion given a reason-able set of assumptions. Recently Sartori et al. (2019)developed a code that is capable of simulating Edding-ton ratio curves with a duration of Myr to Gyr, and atime resolution of 10-100 days. The inputs to the codeare a probability density function (PDF; in this case weassume it is the Eddington ratio function) and the powerspectrum density (PSD). The assumption that PDF isgiven by a full Eddington ratio function, rather then bya lognormal distribution with a given σ , is the main dif-ference from earlier modeling work (e.g., MacLeod et al.2010). We do not attempt in this work to distinguishbetween the two models. We show here an example of how the observed depen-dence can be used to constrain the PSD parameters. Wemodel the PSD as a broken power, i.e., with P SD ( f ) = A × (cid:20)(cid:18) ff br (cid:19) α low + (cid:18) ff br (cid:19) α high (cid:21) − (1)where f br is the break frequency, and α low and α high are the slopes at lower and higher frequencies, re-spectively (longer and shorter timescales, respectively).While there is agreement in the community that α high ≈ <
10 days, e.g.,Edelson et al. 2014), the deduced values for α low and f br vary greatly depending on the survey and methodused (e.g., MacLeod et al. 2010, 2012; Graham et al.2014; Koz(cid:32)lowski 2017). Physically, the determination of f br is of great interest as it would provide us with a clueabout the physical scale on which the properties of AGNaccretion change.In the left panel of Figure 5 we show the expectedmean change of the measured brightness during 14.85years, the average time difference between two measure-ments, as a function of f br and α low . Changing theseparameters effectively changes the “burstiness” of theAGN accretion episodes and therefore influences how bservational Nonstationarity of AGN variability l o g ( f r e q . o f t h e P S D b r e a k / H z ) l o g ( o b s e r v e d E dd . r a t i o ) (g,r,i)/mag Figure 5.
Left: the mean change in measured brightness for AGN sampled at 0.5 mag above the brightness cut of a hypotheticalsurvey, and measured again 10 rest-frame years later, as a function of α low and f br . Right: typical simulated “light curves”(observed Eddington ratio) curves, from each of the areas denoted with a small cross in the left panel. Colors correspond to thecolors in the left panel, where we use orange (instead of white) to color the curves from the middle, observationally plausible,region. We show three simulated curves from the observationally plausible region to demonstrate the diversity of behaviors,reminiscent of the observed diversity of AGN variability. quickly the AGN are changing their luminosity in afixed time period. This plot has been made for thesystems selected with an Eddington ratio cut 0.2 dex(0.5 mag) above the break of the Eddington ratio distri-bution, which broadly mimics the SDSS observationalcut. We can see that the observed mean brightnesschange defines a very specific range of allowed valuesin this parameter space. The two parameters are some-what degenerate - observationally, when using a limitedamount of data points, there is little difference if theprocess decorrelates quickly at longer timescales (small α low and large f br ) or slowly at shorter timescales (large α low and small f br - see also Figure 12 in Caplar &Tacchella 2019). In the right-hand side of 5 we showrepresentative “light curves” from different regions ofthe parameter space. In actuality we generate curvesthat satisfy observed Eddington ratio distributions, andwe make an assumption that variability in these “lightcurves” is equivalent to the variability in the observedlight curves. In particular, we emphasize the wide va-riety of the behaviors for the curves that are consistentwith the observed changes in the mean brightness. Thisis reminiscent of the wide diversity of observed variabil-ity behaviors for AGN.Qualitatively, as indicated before, this model also ex-plains why the most luminous objects at the lowest red-shift are more likely to get dimmer. As they alreadyoccupy the uppermost edges of the probability densityfunction (Eddington ratio distribution) when they wereobserved in SDSS, they are far more likely to get dim-mer and move to more common regions of the parameter space. In other words, for the brightest AGN, the onlyway to go is down!In the future, we aim to improve observational con-straints and finely map time dependence by incorporat-ing information from various surveys, such as POSS,Pan-STARRS, Zwicky Transient Factory and GAIA.These surveys do not achieve such depth as HSC, butmonitor the sky with high cadence. We will measure thechange of brightness as a function of time, while mod-eling the effect of the incompleteness that arises whenstudying AGN variability and AGN dimming in shallowsurveys. We aim to use this information describing theobserved bias to distinguish between the models withdifferent PDFs (full Eddington ratio function or lognor-mal distribution) and place fine constraints on the evo-lution of properties (primarily power spectrum density)describing AGN variability.ACKNOWLEDGEMENTSWe thank the anonymous referee for providing com-ments. During the preparation of this manuscript, wehave benefited from useful discussions with LaurentEyer, Andy Goulding, ˇZeljko Ivezi´c, Robert Lupton,Lauren MacArthur, Chelsea MacLeod, Sophie Reed, LiaSartori, John Silvermann, and Krzysztof Suberlak. Weespecially thank Yusra AlSayyad, who prepared the fil-ter flags used when retrieving the HSC data. We thankˇZeljko Ivezi´c and Christopher Kochanek for pointing outadditional theoretical possibilities for explaining the ob-served effect.This research made use of NASA’s AstrophysicsData System (ADS), the arXiv.org preprint server, thePython plotting library matplotlib (Hunter 2007) and Caplar et al. astropy , a community-developed core Python packagefor Astronomy (Astropy Collaboration et al. 2013).REFERENCES, a community-developed core Python packagefor Astronomy (Astropy Collaboration et al. 2013).REFERENCES