AGN variability in the age of VRO

Over the next ten years, the Vera C.\ Rubin Observatory (VRO) will observe $\sim$10 million active galactic nuclei (AGN) with a regular and high cadence. During this time, the intensities of most of these AGN will fluctuate stochastically. Here, we explore the prospects to quantify precisely these fluctuations with VRO measurements of AGN light curves. To do so, we suppose that each light curve is described by a damped random walk with a given fluctuation amplitude and correlation time. Theoretical arguments and some current measurements suggest that the correlation timescale and fluctuation amplitude for each AGN may be correlated with other observables. We use an expected-information analysis to calculate the precision with which these parameters will be inferred from the measured light curves. We find that the measurements will be so precise as to allow the AGN to be separated into up to $\sim 10$ different correlation-timescale bins. We then show that if the correlation time varies as some power of the luminosity, the normalization and power-law index of that relation will be determined to $\mathcal{O}(10^{-4}\%)$. These results suggest that with VRO, precisely measured variability parameters will take their place alongside spectroscopy in the detailed characterization of individual AGN and in the study of AGN population statistics. Analogous analyses will be enabled by other time-domain projects, such as CMB-S4.

The rich and ubiquitous property of AGN variability has been exploited not only to classify [24,25] but also to identify AGN [26][27][28][29][30]. However, attempts to connect variability measurements to physical mechanisms have lacked sufficient data, leaving such connections mostly tenuous [31,32]. More success has been found in modelling the stochastic nature of these processes [33]. Such modelling has seen that most AGN exhibit variability that is well described by a damped random walk (DRW), as shown by analyses involving structure functions, autocorrelation functions, and power spectra [34][35][36][37][38], although, some examples of non-DRW AGN have been found [39][40][41]. The DRW model has been used to extract the variability timescale through the use of the structure function [42]. In addition, numerical investigations have been used to model AGN variability across multiple timescales [43,44]. Recently,the variability of 67 AGN were characterized to follow a power law with an index measured to one part in ten [45].
In this paper, we explore the prospects for studying AGN variability with VRO. Previous analyses of AGN variability have been limited due to small AGN sample size or infrequent visit times, both of which will be remedied with the dawn of VRO. For a single AGN, we forecast that VRO will measure both the variability amplitude and timescale up to 10 3 σ. With such precise measurements, we then model a power law relationship between an AGN's bolometric luminosity and its variability with index β b in VRO frequency band b. We find that this index will be measured possibly to one part in a million. Therefore, these results suggest that the next decade of observations will lead to a wealth of knowledge in AGN variability. This paper is structured as follows. In Sec. II we present the formalism for measuring the variability amplitude and timescale in the context of a single AGN and arXiv:2110.13149v1 [astro-ph.GA] 25 Oct 2021 a population of AGN. We then follow up this formalism in Sec. III and present estimators in order to quantify variability. Moreover, using these estimators we make an expected information forecast to VRO's sensitivity in measuring both variability parameters for a single AGN, along with their covariance. Then, we present the analogous calculation for the theoretical best sensitivity to measuring a power law relationship between AGN bolometric luminosity and variability. We discuss these results and conclude in Secs. IV and V, respectively.

II. FORMALISM
Assume the intensities from a population of AGN have been measured over time. We pursue a description of the variability of these intensities through the use of the two-point correlation of their intensities. In this vein, we first present the autocorrelation function for a single AGN. Then, for a population of AGN, we extend the presentation of a single AGN and model a relationship between the bolometric luminosity of an AGN with both its variability parameters through a power law.

A. Single AGN
Let I j b (t) be the observed intensity of AGN j in frequency band b andĪ j b be its time average . With these two quantities, define δ j b (t) = I j b (t)/Ī j b − 1 to be the observed variability of AGN j in frequency band b. We describe the statistical properties of the observed variability of AGN j in frequency band b in terms of the observed two-point variability correlation function δ j b (t 1 )δ j b (t 2 ) . If the underlying mechanism creating the observed signals is independent of time for the duration of observation, then the two-point correlation function is homogeneous in time (i.e. stationary) and thus only a function of the time lag t = |t With this assumption, we model the observed twopoint variability correlation function for a single AGN as with ξ j bb (t) the two-point variability correlation function taking the form of a damped random walk [34,36,46], A jb the variability amplitude andt o jb the variability timescale of AGN j in the observer's frame, ∆t b the temporal resolution of the experiment, and δ(t) the Dirac delta function. Moreover, σ 2 jb is the variance of a white noise process representing photometric error in an AGN's intensity measurements. We assume this noise does not correlate with any AGN's variability. Note that the observed two-point variability correlation function contains The dimensionless variability power spectrum ωP j bb /(2π) for an AGN with variability amplitude A jb = 1 using Eq. (4). The peak of this power spectrum occurs at ω =t −1 jb , with amplitude A 2 jb /(2π). For frequencies smaller than this peak it rises as ω, and for frequencies larger it falls as ω −1 . Since we plot the angular frequency in units oft −1 jb , its value is arbitrary.
instrumental noise, while the two-point variability correlation function does not. Due to either AGN physics or instrument properties, all quantities mentioned depend on the observing frequency band b. Furthermore, due to cosmic redshifting, the variability timescalet r jb of an AGN located at redshift z as measured in its rest frame is related to the observer frame analog through the expressiont o jb = (1 + z)t r jb . With this expression of the observed correlation function, its Fourier transform is with P j bb (ω) the variability power spectrum for AGN j in frequency band b. Here and in what follows, we use Fourier convention f (t) = (2π) −1 dωe −iωtf (k) and f (ω) = dte iωt f (t). We plot an example variability power spectrum for a single AGN in Fig. 1.

B. AGN Population
In addition to modelling each AGN individually, we also model the two-point variability correlation function of a single AGN from a set of population parameters. Define ξ bb (t, z, L) to be the variability correlation function for an AGN located at redshift z with bolometric luminosity L, In writing this expression, we implicitly assume that AGN variability is only characterized by two parameters: its redshift and bolometric luminosity. In addition to this correlation function, we define the power spectrum P bb (ω, z, L) analogously. Typically, AGN with higher bolometric luminosities are more massive. Since variability on timescales smaller than the light crossing time of the emitted object is suppressed, the larger an AGN, the larger its expected variability timescale. Thus, we model the relationship between AGN variability and bolometric luminosity as with L b = L bol (m b lim , z min ) a normalization constant chosen to be the dimmest expected observed AGN. In this expression, we assumed all AGN to have the same variability amplitude for simplicity. Thus, a population of AGN is described by the three parameters A b ,t r b , and β b . Recently, 67 AGN were found to follow a similar variability timescale relation, with the mass of the AGN as the the only dependent parameter, and the index β ∼ 0.23 [45].
Let dN AGN /dzdL be the redshift and bolometric luminosity distribution of this population of AGN. Then these AGN are distributed throughout the Universe according to with dn(z, L)/dL the AGN luminosity function, dV (z)/dz = 4πf sky r(z) 2 dr(z)/dz the comoving volume observed over a fraction f sky of the sky, r(z) = z 0 |dr/dz|dz the comoving radial distance to a redshift z, dr/dz = −c/(1 + z)H(z) its redshift derivative, c the speed of light, and the Hubble parameter. We use Planck 2018 ΛCDM parameters H 0 = 2.18 × 10 −18 s −1 and Ω m = 0.315 [47], along with the Full AGN luminosity function in Table 3 from Ref. [48].
Given a cosmological distribution of AGN, only those that appear bright enough will be observed. More specifically, given a limiting apparent magnitude m lim b in a band b, the distribution of observed AGN in that band is with Θ(x) the Heaviside theta function, bolometric luminosity of the source, d L (z) = (1 + z)r(z) the luminosity distance, δν b the frequency bandwidth of band b, and F AB = 3.631 × 10 −23 WHz −1 m −2 . Note that in this expression for the bolometric intensity, we assume the observed intensity is roughly constant across the entire bandwidth, and that the redshifted frequency does not alter the intensity in each band significantly. In general, the bolometric correction within the optical range is a function of the bolometric luminosity of the source. However, across all bolometric luminosities the correction changes only up to 20%, and inversion of this expression can only be done numerically. Thus, for simplicity, we adopt that K b (m b ) = 10 for all magnitudes and bands [48]. We plot the observed AGN distribution, without the theta function, in Fig. 2.

III. FORECASTS
One may measure both the variability amplitude and timescale of an AGN from measuring only its variabiliy correlation function at different lag times. However, the relation between the observed correlation function at any lag time and the true underlying stochastic process becomes increasingly inaccurate for variability timescales smaller than the cadence and larger than the observational period. On the other hand, the observed power spectrum is accurate for all Fourier modes well within these limits.
In this section, we use the expected information from the power spectrum estimators for a single AGN and a population of AGN to forecast VRO's ability to measure various variability parameters. The analysis, discussed below, leads to the signal-to-noise results for the variabil-ity amplitude in Fig. 4. Furthermore, we show the covariance between the variability amplitude and timescale in Fig. 5 and the variability index and timescale in Fig. 6. All figures are done for a representative sample of VRO's frequency bands.
We model the individual band errors as a sum of Poissonian shot noise from source, Poissonian shot noise from the sky, Gaussian instrumental noise, and systematic error. More specifically, for VRO we use the fit given by Ref. [23] and propagate the error from apparent magnitude to variability,

A. Single AGN
For notational simplicity, we assume all AGN are observed for the same duration T j = T and sampled at the same times. However, we allow for different sampling between different frequency bands b. Under the null hypothesis, AGN undergo no variability and thus the noise for the power spectrum estimator of AGN j in band b is We also combine the information of all bands through the use of bolometric corrections. First we calculate the bolometric correlation function as given by a particular band. Since the bolometric intensity is approximately linear in the band intensity, their fractional errors are the same. Then, we inverse-variance weigh each band to Eq. (13). The magnitudes plotted ranges from the 5σ apparent magnitude limit in the corresponding band, shown in Table. I, to the theoretical value for the brightest AGN that will be observed m b cut = 15.7.
obtain an estimate of the actual bolometric correlation function. Therefore, the error σ jbol null in measuring this correlation function is Since the bolometric band is a combination of measurements done in different bands with different temporal resolutions -the bolometric band has unequal, but periodic, temporal spacing in measurements. Moreover, not every temporal spacing has an equal number of measurements. Rather than model this spacing, we take an equal-time temporal resolution ∆t bol = T /(n vis − 1), with the condition that null-hypothesis forecasts using this resolution are upper bounds. Moreover, we note that given the bolometric correlation function, we can invert the bolometric corrections in order to translate the bolometric error into the error in any particular band b. Thus, through inverse variance weighing, the bolometric band represents the optimal sensitivity for any particular band.
With the null-hypothesis power spectrum noise in hand, we use the expected information matrix to infer the covariance matrix for our AGN parameters. We plot the signal-to-noise of measurements for the variability amplitude for a single AGN under the null hypothesis in Fig. 4.
To calculate the covariance between the variability amplitude A jb timescalet jb once a signal is detected, we must include the correlations from the signal. Therefore, the noise for the power spectrum P j bb (ω) estimator is now Under the non-null hypothesis, there is a covariance induced in the Fourier amplitudes inferred between different bands. Therefore, in order to asses the ability of VRO to synthesize information from different bands, we assume that all measurements are now done with a cadence ∆t b = T /(n vis − 1) and a single intensity error. Using this resolution, we plot the covariance between the variability amplitude and timescale in Fig. 5. In practice, the AGN shown in Fig. 5 are not affected by these assumptions given that we assume that only AGN that are detected at high signal to noise are included in the analysis.

B. AGN Population
Given a set of individual AGN measurements compromising an AGN population, we also infer the precision When the variability amplitude drops below the noise threshold, as indicated by the dot-dash line, the error becomes too large and the measurement fidelity significantly drops. If the variability timescale is larger than the observation time, as indicated by the dotted line, then all intensity measurements are maximally correlated and the SNR saturates to a constant signal. On the other hand, if the variability timescale is smaller than the temporal resolution, as shown by the dashed line, then each measurement is maximally independent and thus the SNR saturates once more.
with which we can measure the variability-timescale relation in Eq (7). Thus, we again carry out an expected information analysis using the power spectrum, but now parametrized by population parameters A b ,t r b , and β b and present the results in Fig. 6. Since we assume each AGN in this population is described by the same population parameters, the expected information is now the integral over the expected information gained from each of these AGN.

IV. DISCUSSION
Five assumptions are worth clarifying. First, we assumed that the AGN variability correlation function between two temporal measurements at t 1 and t 2 is only a function of the time lag t = |t 2 − t 1 |, i.e. variability is a stationary process. While this is often the case, non-stationarity has been found to exist under certain circumstances. If non-stationarity is a property of a particular class of AGN, then statistics such as the structure function or Wigner function may be utilized instead of the correlation function. The black circles indicates 1σ (68%) confidence, and the yellow 2σ (95%). We note that these results hold for most AGN magnitudes and VRO frequency bands, as the AGN included in this analysis are all assumed to be detected at high signal to noise. For low amplitude and variability timescale, only modes in the white noise regime, P ∝ (A 2t )ω 0 , of the power spectrum are probed, and so there is negative correlation between the two parameters. As the two parameters increase, the red noise regime , P ∝ (A 2 /t)ω −2 , leads to a positive correlation.
Second, we modeled the correlation function using a damped random walk model, which as we stated previously, is not accurate for all AGN classes. However, for any two parameter model the forecasts presented should be accurate to within orders of unity. Models that include a third parameter, such as a damped random walk with an additional break in the corresponding power spectrum between the white and red noise regimes, will only reduce the fidelity of measurements of the variability amplitude and timescale and are outside the scope of this paper.
Third, we assume that the relationship between an AGN and its observed variability timescale can be described by two parameters: its redshift and bolometric luminosity. In reality, we expect other AGN parameters, such as its color, to also play an important role in determining the timescale within a class of AGN. Such a description of an AGN's variability timescale, while important and necessary for a complete description, is outside the scope of this paper.
Fourth, we assumed that the observed frequency of light in a given band is the result of emitted light in the same frequency band. In reality, it is possible that light emitted in a higher frequency band will redshift across lower bands -leading to the final signal be a sum over different frequency bands. As a result, the autocorrelation of a single observed band will be the result of a The covariance between the variability timescale index and norm in VRO's i band in terms of the fractional differences (β − β fid )/β fid and (t r −t r fid )/t r fid . We take fiducial parameters βi ∈ {0.1, 0.23, 0.37, 0.5} andt r i = {30 days, 100 days}. The black circles indicate 1σ (68%) confidence, and the yellow 2σ (95%). We note that these results hold for all VRO frequency bands, as VRO is limited not by instrumental noise. We take Ai = 1. Since increases in both the index β and the normtr increase the observed variability timescale, they are anti-correlated.
cross correlation of emitted bands. Moreover, while we focused on variability two-point functions within a given band, the cross correlation between bands of VRO, as well as between VRO and other experiments will yield even more information about the structure of the AGN. Time lag measurements between UV/optical light and Xrays have already been used to measure the size regions such as the dust torus and broad-line region. We leave all such calculations for future work.
Lastly, we assumed that the true power spectrum can be recovered through measurements of the power spectrum in a finite box with finite resolution perfectly. For an actual experiment, we expect that measurements of the true power spectrum at Fourier modes close to either limit to be degraded. This degradement can be added in our expected information analysis through the introduction of an additional source of error. However, such error only has an effect on our final result when the variability timescale of an AGN becomes close to either limit. For a population of AGN, a bulk of them will most likely have variability timescales greater than a few days and less than a few years. As a result, we expect such degrading to not have a drastic impact on our results.

V. CONCLUSION
In this paper we presented a general framework for measuring the variability amplitude and timescale of any AGN. First, we measured the variability for each AGN and from there construct estimators for the variability correlation function.
Since each timescale estimator was created using the power spectrum at two distinct modes, this introduced covariance between each timescale estimator. However, despite this covariance matrix being non-diagonal, we were able to calculate its inverse. Then, with each timescale estimator and the corresponding covariance matrix, we created a single estimator for the variability timescale using inverse covariance weighting. With an estimator for the variability amplitude and timescale, we then used linear error propagation to calculate the covariance matrix between these two parameters from the initial variability two-point functions. Using this covariance matrix, we forecasted the sensitivity of a VRO to measuring these parameters. We found that both the variability amplitude and timescale will be able to be measured up to 10σ across all bands.
Finally, we calculated the theoretical best sensitivity to a VRO-like experiment measuring a relationship between the luminosity of an AGN and its variability amplitude and timescale. Namely, we used a logarithmic power law model between the luminosity of the AGN and its variability parameters. We found its index to be measured with at least 10 4 σ fidelity.