THE X-RAY VARIABILITY OF A LARGE, SERENDIPITOUS SAMPLE OF SPECTROSCOPIC QUASARS

In order to obtain a list of quasars observed by Chandra, we searched the Chandra archive to find all ACIS-S or ACIS-I observations of SDSS Data Release 5 (DR5) quasars (Schneider et al. 2007) that used no gratings and were public as of 2009 January 13. We use this sample to design and calibrate our complex data-analysis pipeline for an initial study. In a follow-up study, we will expand the sample to include recently observed quasars and additional AGNs selected using new metrics such as optical/UV variability. For each spectroscopic quasar, we identified candidate Chandra observations that had telescope aim points within 15 arcmin of the given quasar. Then, we investigated these candidate observations to determine whether a quasar fell on an active detector chip. In cases where a source fell within 32 pixels of a chip boundary, we discarded the candidate observation because of increased uncertainty of the instrument response in those regions. Where necessary, we corrected aspect offsets using the prescription available online.⁴ Our final sample of multiply observed sources includes 763 Chandra observations of 264 SDSS quasars.

We reduced each observation of an SDSS quasar using CIAO 4.1.2⁵ following the procedures listed online for the tool psextract.⁶ This tool does not handle cases where zero counts are present in the source extraction region. We flagged these cases for special handling. For each quasar observation, we generated instrument response files (auxiliary response files (ARFs) and redistribution matrix files (RMFs)) either directly (using the mkarf and mkrmf tools for zero-count sources) or indirectly (through psextract). These response files include a correction for the buildup of contaminant on the ACIS chips. The instrument response accounts for spatial and temporal variation as a function of detector position and time.

For point sources brighter than r = 20 mag, SDSS astrometry is accurate to ∼50 milliarcsec,⁷ or about 10% of a Chandra pixel. We can therefore use SDSS astrometry to localize sources on Chandra CCDs. We performed "forced photometry" on the X-ray images, extracting source counts from a circular region with radius equal to the 90% encircled energy fraction at 1.5 keV for a given off-axis angle. This was done even in cases where the source is not detected in an image. Extraction radii were computed using data tabulated on the Chandra X-ray Center Web site.⁸

Backgrounds were generally extracted from an annular region centered on the source position. We selected the inner and outer radii of the annulus to be 15 + r_s and 45 + r_s pixels, respectively, where r_s is the source extraction radius. In our experience, this prescription provides sufficient area to estimate backgrounds reliably at small off-axis angles without becoming large enough to be contaminated by numerous sources at large off-axis angles. In cases where background regions fell partly off-chip, we visually selected a new, circular region that was near to the source and appeared not to be contaminated by other sources. In addition to visual inspection for source contamination, we also checked each background programmatically, searching for sources detected by the Chandra tool wavdetect that fell inside or near our background regions. We generated a list of wavdetect sources using a conservative sigthresh threshold of 10⁻⁵, which roughly corresponds to 10 false detections per chip. In cases where sources were found to contaminate our background regions, we selected a nearby background region that was free of contamination.

Figure 1 shows the distribution of the number of times a source was observed by Chandra in our sample. Of 264 sources, 182 were observed two times, while the remaining 82 were observed three or more times. Four sources were observed 15 or more times. Figure 2 shows example light curves for these sources.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Some analyses depend on the time-domain and redshift coverage of the archive data. To illustrate this coverage, Figure 3 shows the maximal rest-frame time span (max Δt_sys) between observing epochs for each source in our sample. Time spans between archived observations range from hours to ≳5 yr for lower-redshift sources. For sources with z > 2, time spans of up to ≈1 yr are covered. Of course, we can also examine timescales shorter than the maximum value of Δt_sys for sources observed more than two times.

Zoom In Zoom Out Reset image size

In an earlier study, Gibson et al. (2008b) found that Chandra detects nearly all (≈100%) non-BAL SDSS spectroscopic quasars in observations having exposure times >2.5 ks or off-axis angles <10 arcmin. (The Chandra detection rate of SDSS spectroscopic quasars is lower, but still very high, for sources with shorter exposures and/or larger off-axis angles up to 15 arcmin (Gibson et al. 2008b).) Motivated by these criteria, we define a high-quality sample, which we call "Sample HQ," consisting of 167 sources. To construct Sample HQ, we culled all observations of any source that had exposures <2.5 ks and off-axis angles >10 arcmin. We also required that sample-HQ observations be performed at an off-axis angle >1 arcmin, in order to eliminate bias that could be introduced by the target-selection criteria of the Chandra observatory. A typical source has about 5 more counts per epoch (64 versus 59 counts, on average) in Sample HQ than in the full sample.

For each quasar, we calculate (or estimate) the monochromatic luminosity $L_{\nu }({\rm 2500 \,{\rm \AA}})$ in units of erg s⁻¹ Hz⁻¹, which we denote L₂₅₀₀. This value is either calculated directly from our fits to the quasar continuum at (rest-frame) 2500 Å or is extrapolated by normalizing the composite quasar spectrum of Vanden Berk et al. (2001) to match the observed spectrum in the SDSS bandpass. In the former case, continuum fits were performed using the method of Gibson et al. (2009), in which a reddened power law or a polynomial was fit to each spectrum, excluding regions with broad emission lines, BALs, or ionized iron emission. A plot of L₂₅₀₀ as a function of redshift is shown in Figure 4 for our full sample. The bulk of the sample spans a factor of about 40 in luminosity, from 10^29.8 < L₂₅₀₀ < 10^31.4 erg s⁻¹ Hz⁻¹. Because the SDSS survey is flux-limited, L₂₅₀₀ is strongly correlated with redshift. Throughout this work, we use quasar redshifts reported in the SDSS quasar catalog. These software-generated redshifts were verified by visual inspection. Repeat observations of SDSS quasars have shown an rms difference in redshift of 0.006 (Wilhite et al. 2005); this level of precision is sufficient for our study.

Zoom In Zoom Out Reset image size

We also make use of central black hole masses (M_BH) determined from broad emission line fits to H_β, Mg ii, and C iv emission lines for SDSS spectroscopic quasars (Shen et al. 2008). We adopt these values with the understanding that accurate determination of M_BH is an area of active research, with known discrepancies among existing methods. Shen et al. (2008) note that measurement errors for the quasar, continua, and emission lines are generally dominated by statistical uncertainties (≳0.3–0.4 dex) in the calibration of M_BH-estimation methods as well as unknown systematic effects in the use of, e.g., C iv emission as an estimator. Figure 5 shows estimated black hole masses and redshifts for quasars in our sample. For comparison, we have also plotted black hole masses for Seyfert AGNs and quasars determined through reverberation mapping and other methods (Bian & Zhao 2003; Botte et al. 2004; Peterson et al. 2004; O'Neill et al. 2005). (We note that in some cases, these masses may be controversial; this plot is simply intended to be illustrative.) Our SDSS quasar sample extends well beyond the Seyfert AGN regime to higher black hole masses (M_BH > 10⁹ M_☉) and, of course, higher redshifts.

Zoom In Zoom Out Reset image size

We identify subsamples of quasars that are known to be radio loud, following the method of Gibson et al. (2008b). Radio fluxes from FIRST (Becker et al. 1995) or NVSS (Condon et al. 1998) were used to determine the ratio of flux densities at 5 GHz and 2500 Å. Sources having a ratio of R* > 10 (log (R*) > 1) were flagged as radio loud. We classify 25 quasars as radio loud by this criterion. For 97 quasars, we are not able to constrain log (R*) < 1 given the approximate 1 mJy (2.5 mJy) limit of the FIRST (NVSS) surveys. Almost all (93) of these quasars are at most radio-intermediate (RIQs; 1 < log (R*) < 2), which generally have X-ray properties similar to radio-quiet quasars (e.g., Miller et al. 2011), at least in single-epoch analyses. Based on the fraction of radio-loud quasars in the subset of sources that we can classify unambiguously (i.e., radio sensitivity extends to log (R*) < 1), we roughly anticipate ≈7 RIQ contaminants in the radio-quiet quasar set.

We used the SDSS DR5 catalog of BAL quasars (Gibson et al. 2009) to identify any sources known to host absorbing BAL outflows along the line of sight. BAL quasars typically have strong X-ray absorption, so we treat them separately. We expect some contamination from unidentified BAL quasars at lower redshift (z ≲ 1.7) because the C iv BAL-absorption region for these sources does not lie in the SDSS spectral bandpass; BALs are identified in ∼15% of higher-redshift AGNs (e.g., Hewett & Foltz 2003; Trump et al. 2006; Gibson et al. 2009). In the full sample of 264 sources, we unambiguously identify 18 cases of BAL quasars.

For each ACIS chip that observed an SDSS quasar, we construct a weighted exposure map using the prescription given at the Chandra X-ray Center Web site.⁹ This exposure map can be used to estimate the source count flux from the count rate observed in a given epoch. We used a Galactic-absorbed power law for the exposure map weights. In order to roughly estimate the shape of the ACIS spectrum obtained for each quasar observation, we fit each background-subtracted spectrum with a power law. This model was fit to counts in the observed-frame 0.5–8 keV energy band using the Cash statistic (Cash 1979). Our goal was not to model features in each spectrum, but simply to describe the overall spectral shape in terms of a photon index (Γ) that would be used to construct the weighted exposure maps. We constructed exposure maps for photon indices of 1.6, 1.8, 2.0, and 2.2, and selected the exposure map corresponding to the value of Γ that most closely matched our fit value for each observation.

We extracted total counts from a circular aperture and estimated background counts from an annular (or offset circular) aperture as described in Sections 2.1 and 2.2. Both total and background counts were multiplied by a factor of 1/c_A as an "aperture correction." Because the fraction of encircled counts depends on photon energy, we used a different correction factor for our soft (0.5–2.0 keV), hard (2.0–8.0 keV), and full (0.5–8.0 keV) bands. We calculate the correction factor for each observation by determining the median energy of a count in the appropriate band. The correction factor c_A does not differ greatly from the aperture correction c_A = 0.90 corresponding to the 90% encircled fraction for 1.5 keV photons for our extraction radii. For the soft, hard, and full bands, typical values of c_A are 0.92, 0.87, and 0.91, respectively.

Source counts were estimated by subtracting the background from the total number of counts. We estimated the 1σ upper and lower limits on the number of source counts by applying Equations (7) and (11) of Gehrels (1986) to the numbers of total and background counts, then propagated the errors to determine errors in source count rates. We applied a small correction factor to account for the influence of Galactic absorption in each band; the correction factor was determined using the Galactic column with our power-law model fits to each spectrum. We then divided the number of counts by the average value of the weighted exposure map in the aperture to obtain count rates (in counts cm⁻² s⁻¹) in the soft, hard, and full bands for each observation of an SDSS quasar. Properties for individual sources are given in Table 1, while properties measured for individual epochs are given in Tables 2 and 3.¹⁰ The count rates used in this work are measured in the observed frame unless otherwise specified. Our statistics are often dimensionless, so if we apply a factor of (1 + z) to the time dimension when we calculate count rates or derived quantities such as count-rate errors, this factor would cancel out of the final statistic. See, e.g., Equations (2), (4), and (5). For many of our analyses, we do not classify sources as "detected" or "not-detected" (at some confidence level), but instead work directly with the background-subtracted count rates. For faint sources, these rates can be zero or even negative. In any case, there is a high detection rate in our sample. For example, full-band detections were obtained at >95% confidence for 690 of 763 observations in the full sample and 490 of 507 observations meeting Sample HQ requirements. Furthermore, for some analyses, we restrict the sample to sources with higher numbers of counts in order to maximize signal-to-noise ratio (S/N).

We calculate "count-rate luminosities" from 0.5–8 keV count rates as

$$\begin{eqnarray} L_i &\equiv & P_i (1+z)^x 4\pi D^2_L, \end{eqnarray} \tag{ 1 }$$

where P_i is the flux in counts cm⁻² s⁻¹ (calculated by dividing net source counts by the exposure map), D_L is the luminosity distance, and the exponent x incorporates both a bandpass correction and a K-correction to account for the fact that the Chandra bandpass is sampling different spectral regions as a function of redshift. We have assumed a power law with photon index of Γ = 1.8 for this correction, so that x = 0.8 (Schmidt & Green 1986). We use a fiducial photon index of 1.8 instead of our best-fit photon indices because we do not want to amplify any effects that are purely due to measurement error or modeling (e.g., absorption features that change dramatically between observations). We do not convert count-rate luminosities into traditional luminosities (in units of erg s⁻¹) because doing so would require assumptions that decrease the accuracy of our results. Calculating the typical photon energy used as a conversion factor in a given bandpass involves making an assumption about the spectral shape or performing a model fit; both of these can introduce additional errors, especially for atypical and/or faint sources.

Figure 6 compares the luminosity and redshift distribution of our sample to that of A00, Pao04, and Pap08. (Data for this figure were kindly provided by O. Almaini and M. Paolillo; the data for the M02 sample were not available.) For illustrative purposes, we estimate L_0.5–8, the luminosity in the 0.5–8.0 keV band, from the count-rate luminosities in our sample assuming that the energy of each photon is ≈1 keV. The Pao04 sample drawn from the CDFS and the Pap08 Lockman Hole sample cover lower luminosities (and shorter timescales), while our sample approximately covers the A00 sample and extends it to luminosities up to ∼10 times higher. Our sample includes more quasars at higher redshifts (with 63 quasars at z > 2), as well. Figure 7 shows the maximum rest-frame timescales probed by each study, estimated for sources in A00, Pao04, Pap08 by dividing the maximum observed-frame timescale by (1 + z). Using the Chandra archive, we can probe timescales up to 10 times longer than in the lower-luminosity sample of Pao04, and ∼100 times longer than in A00 and Pap08. Figure 7 shows (as a gray box) power spectrum break timescales estimated for some typical sources in our sample using the relation given in McHardy et al. (2006). The X-ray archives enable us to explore the regime beyond this break timescale.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For each quasar, we performed a one-parameter fit to the photon fluxes using our estimated errors to determine the best constant flux that would describe our data. As expected, we find that the observed distribution of χ²_ν is skewed to higher values than would be expected if the constant model were a good fit to the data. Of course, several factors affect the distribution of χ²_ν in our analysis, including non-normality of the count-rate distribution and the estimation of source count-rate errors.

For reasons such as these, we adopt a different method to determine whether sources can be robustly classified as "variable." For each epoch of a given source, we calculate the expected total (source + background) count rate in our extraction aperture from the measured background and best-fit constant flux. We flag an epoch as "variable" if the observed count rate is higher or lower than the number of counts corresponding to a deviation from the best-fit value at >99% confidence, according to a Poisson statistic. Any source with at least one variable epoch is considered a "variable source." For the full sample, 74 of 264 sources are classified as variable; 54 of 167 sources are variable in Sample HQ. We classify the remaining sources as "non-variable," with the caution that this term should be understood to mean that we did not detect variability within the limits of our data.

Variability tests are more sensitive for sources with larger numbers of counts. In order to examine this dependency in our data, we plot in Figure 8 the mean number of source counts as a function of redshift. The "mean counts" are calculated for each source using only epochs from Sample HQ that have high exposures and low off-axis angles. Of course, the mean number of counts depends on variability characteristics, so this figure should be considered to be only a rough exploration of our data set. Sources flagged as variable are plotted in red, while non-variable sources are plotted with black points. Sources known to host BALs are plotted as circles, while known radio-loud sources are triangles and known radio-loud BAL quasars are stars. The remaining sources, including those for which radio-intermediate status and BAL absorption could not be ruled out (Section 2.3), are plotted as squares. As a rough guide, we have also plotted solid curves indicating typical numbers of counts for quasars with 0.5–8 keV luminosities of 10⁴⁴, 10⁴⁵, and 10⁴⁶ erg s⁻¹. The curves were constructed assuming a hypothetical source with an unabsorbed, Γ = 2 power-law spectrum observed on-axis on an ACIS-I chip for 18 ks (which is a typical, median exposure time in our sample). Individual sources may differ from this hypothetical source in various respects.

Zoom In Zoom Out Reset image size

Variability is primarily detected in quasars having ≳50 counts per observation (on average), with about half of those sources flagged as variable. Our sensitivity to variability is diminished at higher redshifts (z ≳ 2), where the majority of sources have ≲50 counts. On the other hand, we detect variability in 9 of 14 sources with ⩾700 average counts and 7 of 10 sources with ⩾1000 average counts at redshifts z ≲ 1.2.

In Figure 9, we show the maximum variation from the best-fit constant count rate for the sources in Sample HQ as a function of the mean number of counts per epoch. The y-axis for the plotted points is the maximum of |r/r₀ − 1| for all HQ observations of a source, where r is the measured flux and r₀ is the best-fit constant flux for that source. The solid black line shows the "3σ limit" relation ($y = 3\sqrt{x}/x$ for a mean number of counts x), which is a rough approximation of our variability-selection criterion. The thick, red line indicates the relation y = f(x), where f(x) is the fraction of sources that are identified as variable in the set of sources having mean counts ⩽x. For the entire sample including sources with small numbers of counts, ≈30% of sources are classified as variable at high significance.

Zoom In Zoom Out Reset image size

In the following discussion, we calculate fractional variation using all available measurements for each source. We define the fractional variation, F, using the formula:

$$\begin{eqnarray} F &\equiv & (c_i - c_j) / (c_i + c_j), \end{eqnarray} \tag{ 2 }$$

where c_i is the count flux in the later epoch and c_j is the flux in the earlier epoch. Mathematically, F represents half of the full distance between two measurements, (c_i − c_j)/2, as a fraction of the average value of those two measurements, (c_i + c_j)/2. We have chosen this functional form to be symmetric (up to a sign) in c_i and c_j and to roughly represent the fractional deviation from some "average" flux. Each measurement of F between two epochs is associated with a rest-frame time between measurements, Δt_sys.

3.3.1. Fractional Variation Over Time

Given the complex nature of measurement errors in our data set, we choose to assume that the intrinsic distribution of F is Gaussian. We use that assumption to constrain the distribution of fractional variation given the observed values of F and errors in F. The standard deviation of the Gaussian distribution, σ(F), is calculated using the likelihood method of Maccacaro et al. (1988, hereafter M88) in bins of 100 pairs of epochs. The errors on σ(F) that are shown in the plot were calculated in the same way. The epoch pairs were constructed using all pairs of measurements for each quasar. A quasar that was observed N times would therefore contribute N(N − 1)/2 pairs.

Figure 10 shows our estimated values of σ(F) as a function of Δt_sys. Each point in the figure is placed at an x-coordinate corresponding to the median Δt_sys in the bin for which σ(F) was calculated. The dot-dashed line represents a linear fit of σ(F) at timescales >5 × 10⁵ s and is parameterized by

$$\begin{eqnarray} \sigma (F) &=& (0.000 \pm 0.012) \log (\Delta t_{{\rm sys}}) + (0.156 \pm 0.083). \end{eqnarray} \tag{ 3 }$$

At Δt_sys ≳ 1 day, the fractional variation is about 15.6%. At short timescales (Δt_sys ≲ 5 × 10⁵ s), there is no significant variation detected above our measurement errors. If much more data were available to constrain σ(F) as a function of time, we could map out the gap in the current plot where the fractional variation "jumps" from insignificant (on the shortest timescales) and breaks to a flatter shape at longer timescales.

Zoom In Zoom Out Reset image size

We note that one source, J142052.43 + 525622.4, has so many observations that it dominates one of the long timescale bins. It varies somewhat less than the typical source, and if it is removed, the fit model shows a mild, but insignificant increase with longer Δt_sys. Interestingly, J142052.43 + 525622.4 has weaker [O iii] emission, similar to the trend observed for non-variable sources (Section 3.10), although it is technically classified as "variable" in our sample.

The red dotted lines in Figure 10 indicate the 1σ and 3σ upper limits on σ(F) for 15 epochs of radio-loud, non-BAL quasars with Δt_sys > 5 × 10⁵ s. Because the 3σ upper limit for radio-loud quasars is only a little larger than the standard deviation of fractional variation for radio-quiet quasars, radio-loud quasars appear to be less variable than radio-quiet quasars in our sample. If real, this effect could be due to additional, relatively constant jet emission that dilutes the variable X-ray spectrum of the disk corona.

One shortcoming of this "ensemble approach" to measuring fractional variation is that sources with more observations have more influence on σ(F) because they contribute more epochs to the ensemble average. We can also measure quasar variation using a single epoch from each source and comparing it to the best-fit photon flux for that quasar. This approach has the advantage of placing each quasar on an equal footing, so that a small number of quasars do not dominate the results. However, we do not have enough epochs to map out variation in detail as a function of Δt_sys, and of course there is no particular "time" associated with the best-fit constant flux. If we simply compare one epoch per source to the best-fit count rate for that quasar, we find a time-independent standard deviation of fractional variation (c_i/c₀ − 1, where c₀ is the model flux) of about 16.7% ± 2%. This result is consistent with values of σ(F) found above for timescales ≳1 week, which is reasonable given that such timescales represent the large majority of our data set.

3.3.2. Symmetry of Variation

3.3.3. Extremes of Variation

While our assumption of a Gaussian distribution of fractional variation permits us to characterize the typical extent of variation, it does not describe the behavior of outliers or extreme variation events. Figure 11 shows the fractional variation, F, between pairs of epochs in Sample HQ. For this plot, we have omitted any cases where the fractional variation is <1σ from zero, in order to clearly visualize any extreme values of F. Black points show typical quasars, while red or green points signify known BAL or radio-loud quasars, respectively. As the plot shows, F is nearly always <50%, and F > 100% is apparently quite rare.

Zoom In Zoom Out Reset image size

Extreme variation has been observed in rare cases such as the narrow-line type 1 quasar PHL 1092, which decreased in flux by a factor of ∼200 over ≈3.2 yr in the rest frame (Miniutti et al. 2009). We can use our data set to place limits on the frequency of such events. Figure 12 shows upper limits on the rate at which new observations sharing our Sample HQ properties would be expected to show a given magnitude of fractional variation. The limits are constructed by assembling all values of fractional variation, F, over a given time frame. The plot considers three time frames, Δt_sys < 5 × 10⁵ s, 5 × 10⁵ ⩽ Δt_sys ⩽ 10⁷ s, and Δt_sys ⩾ 10⁷ s. The upper limits for these time frames are plotted as a function of |F| with black, red, and green curves, respectively. For a given value of |F₀|, we calculated the upper limit on the y-axis by determining the maximum intrinsic rate of occurrences of values |F| > |F₀| given the number of observed cases with fractional variation <|F₀|, according to a binomial statistic and using 95% confidence limits.¹¹ In cases where a source was not detected at >95% confidence in an epoch, we conservatively forced the variation to be 100%. This was done for 19 of 1157 epoch pairs.

Zoom In Zoom Out Reset image size

We caution that these upper limits depend on subsample sizes and include unmodeled variation due to measurement errors. They are therefore intended only to constrain the rates at which variation greater than a given rate F occurs. (We could limit the impact of measurement error by dropping sources that have low average count rates, but that could introduce new biases.) For example, fractional variation |F| approaching 100% should occur in <10^−1.4 ≈ 4% of observations. Fractional variation F ⩾ 25% should occur in <10^−0.6 ≈ 25% of observations.

We calculate the excess variance of measured count rates (e.g., Nandra et al. 1997; Turner et al. 1999), as

$$\begin{eqnarray} &&\sigma _{{\rm EV}}^2 \equiv \frac{1}{N\mu ^2} \sum ^{N}_{i=1} \big ((n_i - \mu)^2 - \sigma _i^2\big), \end{eqnarray} \tag{ 4 }$$

where n_i is the count rate in epoch i, μ is the mean of the values n_i, and σ_i is the average statistical error in the measurement of n_i. In order to avoid any bias introduced by having different numbers of epochs per source, we calculate σ_EV using only the earliest and latest observation of each quasar. N is therefore always 2 in this calculation. The excess variance statistic is not ideal for our data set, which consists of sources observed a small number of times over diverse time intervals. However, we consider it for comparison to other published analyses.

We limit our analysis to 131 radio-quiet, non-BAL sources that have ⩾50 counts per epoch, on average, in order to screen out a large number of sources with small or negative σ²_EV. We estimate the error on σ_EV using the formula given in Turner et al. (1999). Figure 13 shows σ_EV as a function of count-rate luminosity. Sources at z > 2 are plotted in red. At lower redshifts (z < 2), 49 of 113 sources have σ²_EV < 0.001; these are plotted as black open circles at y = 0.001. Similarly, 10 of 18 higher-redshift (z > 2) AGNs have σ²_EV < 0.001 and are plotted with red open circles along the bottom of the plot. No clear patterns are visible in the plots except a tendency for higher-redshift sources (plotted in red) to have higher luminosities, due to the flux-limited nature of the SDSS quasar survey.

Zoom In Zoom Out Reset image size

In order to test for additional structure in the data, we calculated the mean of the excess variance values in bins of 21 sources each. To calculate this mean, we used only lower-redshift sources at z < 2. This mean σ²_EV is plotted in Figure 13 as a thick green line. Error bars on the data points (placed at the median count-rate luminosity for each bin) represent the estimated error on the mean. The green line shows a general trend of decreasing σ²_EV with luminosity, as is commonly observed (e.g., A00, M02, Pao04, Pap08, and references therein). A filled green circle represents the mean σ²_EV for 18 sources at higher redshifts z ⩾ 2; its value is less than 0.001. Dashed red lines in Figure 13 indicate 2σ and 3σ upper limits on the excess variance for quasars at z > 2. While we cannot completely rule out the possibility of an increase at high redshifts, the upper limits indicate that any increase would be at most very small.

The extent to which a quasar is measured to vary can depend on the timescales over which it is observed. Local Seyfert AGNs show increasing variations with time up to at least a break timescale. As we discuss in detail in Section 4.2, it would be difficult to account for this effect. We do not know the shape of the power spectrum for quasars, and many of our observations are likely beyond the break timescale (Figure 7). Instead, we explore the possibility that the increase in variability could be timescale-dependent by calculating σ_EV using only observations that are separated by a certain range of timescales. For example, we constructed a subsample of observations with rest-frame timescales ≲1 yr. We tried a variety of system-frame and observed-frame timescales, and in no case did we find evidence for increased variation at z > 2.

Using similar methods, we find no significant correlation between σ²_EV and the average HR (defined in Section 3.6) of a source. This agrees with the result of Pap08, who found no correlation between variability amplitude and spectral slope Γ in their light curves. Pap08 note that this result disagrees with earlier results (e.g., Green et al. 1993) indicating that nearby AGNs with steeper spectra showed larger amplitude variations. While this may be an indication that luminous quasars vary differently than local AGNs, we caution that many factors affect these results, including the timescales probed, the baseline of spectral shapes spanned by a sample, and sensitivity to variation in fainter sources. In the following sections (Sections 3.6 and 3.7), we examine the related issue of how spectral shape changes for a single source as it varies.

The Chandra archive provides repeat observations of quasars having a wide range of luminosities, redshifts, and black hole masses. While individual sources may not have sufficient observations to permit sensitive tests of variability properties, we can characterize the physical dependence on variability more effectively in an ensemble of observations. We place each observation of a radio-quiet, non-BAL quasar in Sample HQ on an equal footing in this ensemble by considering the quantity

$$\begin{eqnarray} S_{ij} &\equiv & (f_{ij} - c_j) / \sigma _{ij}, \end{eqnarray} \tag{ 5 }$$

where f_ij is the 0.5–8 keV count rate observed in epoch i for source j, c_j is the best-fit constant count rate for source j over all epochs, and σ_ij is the error on f_ij. In the absence of variability, we would expect the distribution of S_ij values to be roughly Gaussian with an rms of 1. The mean value of S_ij would be ≈0 (by construction) in ideal conditions where variation is symmetric and not influenced by outliers. This is generally the case except at the lowest redshift, which is influenced by a few (≈5) outliers that have higher S_ij values. We use S_ij to test the hypothesis that intrinsic variability can be detected in the data at a certain redshift, luminosity level, or black hole mass. Of course, determining the amplitude or detailed pattern of any intrinsic variability would require simulations to determine what pattern of S_ij values could be expected from a given variability model.

In Figure 14, we plot S_ij as a function of redshift. We have binned the set of S_ij values into five bins and calculated the sample standard deviation in each bin. Red data points (placed at the median x-value of each bin) indicate the square root of the unbiased sample variance, with error bars estimated as $\Delta \sigma \equiv \sigma / \sqrt{2(N-1)}$ for a bin with N data points. (However, we note that scatter in S_ij can be driven by non-Gaussian effects such as outliers.) The unbiased sample variance reflects the influence of outliers, some of which have been clipped out of the plot range for visibility purposes. A better estimate of the overall sample variability properties can be obtained using the median absolute deviation (MAD; e.g., Maronna et al. 2006), which is much less sensitive to strongly variable outliers.¹² The green line shows this value, σ_MAD calculated for each bin, with errors roughly estimated using the same formula as for the red data points.

Zoom In Zoom Out Reset image size

As redshift increases, the amount of variability that we detect decreases. σ_MAD for S_ij decreases to ≈1 as redshift increases, indicating that any intrinsic variability is dominated by measurement errors in this analysis at z ≳ 2. As Figure 8 shows, the average number of counts per epoch for a given quasar also decreases at higher redshift, so that variability becomes harder to detect at z > 2. However, we find the same result—S_ij decreasing to the noise level (S_ij = 1)—even when we limit our analysis to Sample HQ sources that have 10 to 50 average counts per epoch. We do not, therefore, find any indication that variability increases at higher redshifts in this data set.

Of course, these statements simply express the extent to which we can measure variability in the data. Even bright AGNs at high redshifts are presumably variable at some level, and their variability patterns may be complex. In order to demonstrate how a simple pattern of variation would appear in the data, we have created a set of simulated data in which we fix a constant amount of variation (10%, 20%, or 30%) in each observation of a quasar. σ_MAD for the simulated S_ij values are plotted as blue lines in Figure 14. These decrease with redshift because measurement errors are generally increasing with redshift, but the simulated lines do not fall all the way to S_ij = 1, because they do possess intrinsic variability. The green line corresponding to real quasars declines more steeply with redshift than the blue lines, again indicating that we see no significant increase in variability with redshift.

We find similar results when plotting S_ij as a function of luminosity or M_BH. σ_MAD decreases to the noise level in the highest M_BH bins. σ_MAD decreases with luminosity as well, although there is a small (but insignificant) upturn at the highest luminosity values (Figure 15).

Zoom In Zoom Out Reset image size

Finally, we have attempted to disentangle the degeneracy between redshift and luminosity by binning our data on both quantities. Even with only four redshift bins (0–0.5, 0.5–1, 1–2, >2), this pushes the limits of available data, leaving generally 10–35 S_ij points per (L, z) bin. In any case, we still find a general trend for σ_MAD to decrease with both luminosity and redshift, and no evidence of an increase at z > 2.

Throughout this analysis, we have worked with a variability signal that is measured in the observed-frame 0.5–8 keV band where Chandra is most sensitive. To the extent that AGN variation differs from one part of the spectrum to another, our ability to detect variation is also redshift dependent. (See further discussion of this point in Section 4.2.) The instrumental response varies strongly over the bandpass, further complicating any attempt to separate energy-dependent variation from redshift dependence. For the current work, when we state that "we do not detect variability" in some instance, that statement should be understood in the context of the data that we have available, using the Chandra observed-frame bandpass and instrumental response.

HRs—here defined as the ratio HR ≡ (H − S)/(H + S), where H and S are the observed-frame count rates in the hard (2.0–8.0 keV) and soft (0.5–2.0 keV) bands, respectively—may be used to describe the overall spectral shape for sources that do not have sufficient counts for detailed spectral fitting. We have computed observed-frame HRs for all observations in Sample HQ, using the method of Park et al. (2006) to determine 1σ HR errors. We calculated separate weighted exposure maps for the soft and hard bands using the spectral slope obtained from the best fit to the full-band spectrum of each source. (Narrower bandpasses and smaller count rates do not permit separately fitting the soft and hard-band spectra.)

We examine all epochs in Sample HQ (excluding those from BAL and radio-loud quasars) to determine how HR varies with luminosity. Here we put all epochs of all sources on an equal footing, scaled by their HR and luminosity in the earliest available epoch. For each epoch, we calculate the change in HR and luminosity from the first observation of that source. These quantities [HR − HR(t = 0) and L/L(t = 0)] are plotted in Figure 16 for lower-redshift (z < 2, black squares) and higher-redshift (z ⩾ 2, open circles) quasars. There is a significant anti-correlation (at >99.99% confidence, according to a Spearman rank correlation test) between the change in HR and the L/L(t = 0) for quasars at z < 2. The correlation is not significant at z > 2, although the sample sensitivity is weakened by larger measurement errors for fainter sources as well as the fact that softer photons are shifted out of the Chandra bandpass at higher redshifts. Defining ΔHR ≡ HR − HR(t = 0) and F ≡ L_ph/L_ph(t = 0), we obtain the fit

$$\begin{eqnarray} \Delta {\rm HR} &=& (-0.298\pm 0.030)\times \log (F) + (0.003\pm 0.004) \end{eqnarray} \tag{ 6 }$$

for the full sample including all redshifts.

We constructed toy physical models of the trend between ΔHR and L/L(t = 0). The change in spectral hardness could, of course, be associated with a single power-law model that changes both shape and normalization, following the relation described in Equation (6). The goal of our toy models is to determine whether we can reproduce Equation (6) with only a single variable parameter.

Zoom In Zoom Out Reset image size

We modeled several scenarios, including (1) a double power law with a constant hard component (Γ = 1.3) added to a variable softer component (Γ = 2), (2) a constant power law with variable neutral absorption, and (3) a constant power law with an absorber that varies in ionization level. For each toy model, we attempted to determine a set of fit parameters that caused the HR and luminosity to vary together in a way that reproduced the trends observed in Figure 16. The absorption-dominated models (cases 2 and 3) generally produced a much steeper trend of HR with luminosity than we observed. In contrast, the double power-law model (case 1) reproduced the observed trend reasonably well for a hard, constant component having Γ = 1.3 and a normalization at 1 keV of about 20% of the variable, soft power-law component. This model is shown as a red curve in Figure 16. The impact of a two-component model on observations of high-redshift quasars is discussed further in Section 4.3.

For sources with large numbers of counts, we can fit spectra to determine how spectral shapes and features vary as a function of source brightness. We have selected 16 radio-quiet, non-BAL quasars having at least 500 counts per epoch, on average, and fit them with a spectral model consisting of a Galactic-absorbed power law absorbed by one absorption edge. The edge is constrained to have an energy threshold of 0.739 keV in the rest frame, corresponding to absorption from the O vii ion. The optical depth of the edge, τ, is allowed to vary. We have chosen the edge model to roughly characterize a power-law spectrum with ionized absorption using a single parameter. In fact, we constrain τ > 0 at the 1σ level in only 7 of 87 epochs. We obtain similar results using a neutral, rest-frame absorption model, so the choice of this model is not strongly biasing our results. While the single-edge model certainly does not reproduce the complex spectral features that an absorber may imprint on an emission spectrum, it does provide a basic model of variable absorption when spectra do not permit careful fitting of absorption models that have more degrees of freedom. The results of the fits are given in Table 4.

Figure 17 shows the distribution of (ΔΓ_ij, Γ_ij) pairs from our fit models. We ordered the Sample HQ observations of each source chronologically, then calculated ΔΓ_ij as the difference in photon indices between epoch i + 1 and epoch i for source j. Γ_ij is the photon index from the first epoch (i.e., epoch i) in the pair. Two outlier points have been clipped in this figure for visual clarity. The plot shows a strong anti-correlation between ΔΓ and Γ at the >99.99% confidence level, according to a Spearman rank correlation test. The median value of Γ in this sample is Γ_M ≈ 2.16. (The value Γ = 2.16 is not necessarily representative of the full population; for example, the sources chosen for spectral fitting were selected to have larger numbers of X-ray counts.) Spectra tend to steepen when they are flatter than this value, and tend to flatten when they are steeper. Figure 18 shows how the best-fit photon index changes with the power-law normalization. There is a strong correlation (99.97% confidence, according to a Spearman rank correlation test) between ΔΓ_ij and ΔNorm_ij, where ΔNorm_ij is the difference in power-law normalization values between epochs i + 1 and i. Overall, spectra steepen as they brighten (and/or flatten when getting fainter), consistent with what we have already observed using HRs for a larger sample of sources (see Figure 16 and Section 3.6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Some local Seyfert AGNs show a linear relationship between hard- and soft-band X-ray count rates (e.g., Taylor et al. 2003). When the linear relationship is extrapolated to a zero soft-band count rate, a significant non-zero hard-band flux remains. This suggests a two-component model of Seyfert AGN X-ray spectra, in which a constant hard spectral component is augmented by a softer, variable power law. The non-zero offset in the hard band comes from the underlying hard spectral component that remains when the variable power-law component is absent. The shape of the hard spectral component, determined from flux–flux plots in different X-ray bands, can resemble the expected emission from cold disk reflection (e.g., Taylor et al. 2003).

For each radio-quiet, non-BAL source having two or more observations in Sample HQ, we use the IDL FITEXY routine to generate a linear fit to the observed-frame hard-band (2–8 keV) count rates as a function of soft-band (0.5–2 keV) count rates. That is, we fit for a and b in the equation:

$$\begin{eqnarray} r_H &=& a r_S + b, \end{eqnarray} \tag{ 7 }$$

where r_H and r_S are the hard- and soft-band count rates, respectively. The FITEXY routine accounts for errors in both the hard- and soft-band count rates. Our current data do not permit a precise test of whether the flux–flux relation is best modeled as linear, although visual inspection indicates that a linear model works well in most cases with more than five observations. The linear model is not formally rejected by the fit statistic in 12 of 15 quasars with ⩾5 observations. Because many sources have only two or three observations, our goal is not to test the linear model, but instead to determine a typical value for y₀ in the ensemble that would correspond to a constant hard-band offset.

Under the approximation that the y-intercept values are Gaussian distributed, we used the method of M88 to determine $\overline{y_0}$, the mean value of the y-intercept in these fits. The mean y-intercepts, corresponding to a constant hard-band flux component, were greater than zero at >99% confidence both for lower-redshift Sample HQ sources (51 at z < 1) and also for higher-redshift sources (116 at z ⩾ 1). Then, to estimate the contribution of this hard, constant component to the hard X-ray spectrum, we re-ran the fits, but this time normalized the hard-band count rates by the mean hard-band count rate for that source. The normalized fractional offset has a value of 0.428 ± 0.140 for sources at z < 1 and 0.753 ± 0.165 for sources at z ⩾ 1. The increase in the fractional offset (although with a large error bar) suggests that the shape of a putative constant component would be flatter than the shape of the variable component, so that it contributes an increasing fraction of flux at higher energies. Local Seyfert AGNs show hard-spectrum fractional contributions of 20%–40% at energies ⩾2 keV, possibly due to Compton reflection from the accretion disk; these fractions increase at higher energies (e.g., Taylor et al. 2003; Vaughan & Fabian 2004). If a similar model holds in the quasar-luminosity regime, we might expect the fractional contribution of the constant component to increase as higher energies are shifted into the Chandra bandpass. (See Section 4.3 for further discussion about the effects on observations of high-redshift quasars.)

The counts in Chandra event lists are tagged with arrival times, permitting a search for variability in the light curves of each observation. We use a Kolmogorov–Smirnov (K-S) test to identify any light curve for which the count arrival times are significantly non-uniform. (Here we define "light curve" to mean the time sequence of all counts that fall inside the source extraction region.) At a 99% K-S-test confidence level, we identify 17 potentially variable light curves. We have removed five additional observations from consideration because visual inspection of their light curves indicates that instrumental factors such as background flaring are dominating the light curves. A time-dependent analysis of background variation is beyond the scope of this study, but we note that in our bright sample (with total count rates >0.01 s⁻¹, described below), the estimated background count rate is <10% of the total rate in ≈86% of the light curve exposure. Background effects should not strongly affect our overall results.

Six of the 17 light curves that the K-S test flags as variable belong to a single source, J123800.91 + 621336.0 (hereafter "J1238") at z = 0.44 with absolute i-magnitude M_i = −23.04. This source was observed many times as it falls in the Chandra Deep Field-North (e.g., Alexander et al. 2003). Given that we have searched ≈800 light curves for a result at 99% confidence, we cannot consider the variability criterion to be highly significant for the remaining 11 light curves. This result is consistent with our previous observation that the level of intrinsic variation between observations falls below our detection limit for timescales shorter than ∼10⁵ s (Section 3.3).

Figure 19 shows the 4000–5200 Å region of the SDSS spectrum of J1238 (in black) with the quasar composite of Vanden Berk et al. (2001) overplotted in red for comparison. The FWHM of the H_β and Mg ii λ2800 (not shown) lines is ≈1860 km s⁻¹, while the [O iii] λ5007 line is quite weak in comparison to the composite value. J1238 also shows an excess of ionized iron emission in the 4500–4600 Å region. Fitting a Galactic-absorbed power law to the observed-frame 0.5–8 keV region, we find that the spectrum is steeper than typical quasars, with 2.2 < Γ < 2.4 in most epochs (Figure 20). In several epochs, Γ deviates from the median fit value of Γ = 2.6 by 3.3–5.4σ, including two epochs with flatter spectra (Γ < 2) at lower flux levels. All of these traits indicate that J1238 is a quasar analog of the Narrow Line Seyfert 1 (NLS1) population (see, e.g., reviews in Pogge 2000; Komossa 2007). NLS1 AGNs are often observed to be highly variable in X-rays and are thought to have smaller black hole masses and correspondingly high accretion rates (e.g., Boller et al. 1996; Leighly 1999; Komossa 2007, and references therein).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

J1238 was not detected in the Very Large Array FIRST survey and is classified as radio quiet in our study. It was detected as an unresolved source in the deep radio survey of Richards (2000), with a 1.4 GHz flux of 190 ± 13 μJy, corresponding to log (R*) ≈ −0.004. Although J1238 does show radio emission at faint levels, it is far below the radio-loud limit for our sample.

As a second test for short-term variability, we implemented an algorithm to search for flares or dramatic, short-term absorption (or dimming) events similar to events that have been previously reported at quasar luminosities (e.g., Remillard et al. 1991), or from the Galactic center at much lower luminosities (e.g., Baganoff et al. 2001). The algorithm breaks each light curve into segments 1000 s wide. This segment size was chosen to provide a reasonable balance between the number of segments per light curve and the number of counts per segment. The first stage of our algorithm iteratively determines the baseline count rate of the light curve, excluding any candidate flare regions. In each light curve, we identify the segment that shows the most significant deviation from the mean count rate (according to a Poisson statistic), where the latter is calculated using segments that have not already been flagged as potentially variable and are not the current segment under test. If the deviation is significant at >99% confidence according to a Poisson statistic, we flag that segment as potentially variable. Then we iterate the process to find the next-most-variable segment, and so on. When no more segments can be flagged as potentially variable, we have determined a "clean" estimate of the baseline count rate of the source. Next, the second stage of our algorithm uses this new baseline count rate to determine which segments should be classified as variable at >99% confidence.

Our algorithm identifies hundreds of candidate variable segments out of about 21,652 segments, or 0.69 yr of combined light curves. To evaluate the significance of this result, we re-run our algorithm on simulated light curves that have the same mean count rates and no intrinsic variability. For each segment flagged as variable in either the real or simulated light curves, we record the baseline (expected) count rate and the count-rate "fractional deviation" F_c, defined as

$$\begin{eqnarray} F_c &\equiv & (C_{{\rm obs}} - \langle C\rangle) / \langle C\rangle, \end{eqnarray} \tag{ 8 }$$

where C_obs is the observed count rate in that segment and 〈C〉 is the mean count rate for the light curve.

The distribution of F_c depends on 〈C〉 because weaker variability can be detected at higher count rates. In Figure 21, we show the distribution of F_c for count rates higher than 0.01 counts s⁻¹. We also consider the case of lower count rates separately (not shown). For a baseline count rate of 0.01 counts s⁻¹, a significant "event" would have a flux increase of at least about 80%, or a decrease of about 70%. For higher baseline count rates, the amplitude would be smaller. Events in the high-count-rate subsample could therefore be relatively similar to the flare reported by Remillard et al. (1991) in the quasar PKS 0558–504, which increased by up to 67% over 3 minutes and lasted for 10–20 minutes. Additional rapid flux changes have been observed for this source (e.g., Wang et al. 2001; Brinkmann et al. 2004). A more extreme example would be PHL 1092, which was observed to brighten by a factor of ≈3.8 over a rest-frame time <3.6 ks (Brandt et al. 1999). Rapid, strong flares have also been observed in PDS 456 (Reeves et al. 2000); NLS1-type AGNs such as RX J1702.5+3247 (Gliozzi et al. 2001), IRAS 13224–3809 (Boller et al. 1997; Dewangan et al. 2002); and other AGNs.

Zoom In Zoom Out Reset image size

A K-S test finds no evidence in most cases of a significant difference in the distribution of F_c between the real and simulated data sets either for the distribution as a whole or for the separate cases of absorption/dimming (F_c < 0) and flaring (F_c > 0) events. The difference in distributions is significant (at >99.9% confidence) for the absorption/dimming (F_c < 0) events in the low-counts case; presumably this is because our simulations have not modeled the background and therefore do not reach a "floor" where the counts in a time segment are dominated by the background. Based on these results from the K-S test, we conclude that the distributions of F_c do not differ remarkably from distributions in the null case of no intrinsic flaring or dimming.

Of course, it may still be the case that we are detecting variable segments of light curves at a higher rate than predicted from the simulations, even if the distribution of F_c is not significantly different. Of the four cases under test (low-count-rate absorption/dimming, low-count-rate flaring, high-count-rate absorption/dimming, high-count-rate flaring), two show event rates that are marginally higher than expected from the simulations. According to a binomial statistic, the probability of seeing more absorption/dimming events in the high-counts case is 1.4%, while the probability of observing more flaring events in the low-counts case is 3.5%. In each of the four cases, the flare rate is consistent with zero at 99% confidence. Upper limits on the rates of 1 ks events (after subtracting off the simulated rates, in numbers per observed-frame year) are <6.7 yr⁻¹ (low-count-rate absorption/dimming), <29.3 yr⁻¹ (low-count-rate flaring), <40.5 yr⁻¹ (high-count-rate absorption/dimming), and <9.0 yr⁻¹ (high-count-rate flaring). The limits quoted are for 99% (single-sided) confidence limits. The upper limits are higher in cases that showed marginal evidence for significant event rates above the simulated values because a higher (but not highly significant) number of events was detected in these cases. For the full sample (combining both low- and high-count rates), we have the following rates: <37.3 yr⁻¹ (absorption/dimming) and <51.5 yr⁻¹ (flaring). At most a small fraction of time (≲ 0.2%) is spent by quasars in the short (ks-timescale) flaring or absorption/dimming states that our algorithm is designed to detect.

The algorithm presented here is only a first empirical attempt to characterize any possible flaring activity in the large sample of archived light curves. One possible enhancement would be to search for flares that cover longer timescales. As a preliminary test, we applied our algorithm to search for flares in light curve segments 5 ks long (rather than 1 ks). (At 0.01 counts s⁻¹, a significant event would require an increase by 36% or a decrease by 32%.) Our results are qualitatively similar, but the upper limits we can place are much weaker due to the smaller numbers of light curve segments.

The SDSS spectra of our sources have a wide range of spectral shapes and signal-to-noise values. In order to conduct a basic comparison of the optical spectral properties of variable and non-variable sources in the spectral region that contains the Hβ and [O iii]λ5007 emission lines, we construct composite spectra of all sources with redshifts 0 < z < 0.8 and ⩾50 source counts per epoch, on average. We have selected this region of spectrum because the Hβ and [O iii] emission lines are commonly used to estimate black hole masses and accretion rates, and to identify subclasses of AGNs such as NLS1s. The minimum count-rate requirement ensures that we can detect variability at the ≳50% level at ≈3σ confidence (Section 3.2).

The composite spectra are constructed by calculating the geometric mean of the spectra in a given (rest-frame) wavelength bin. The input spectra are normalized to 1 at 5100 Å in each case; note that the output spectrum depends somewhat on the choice of normalization wavelength. We chose to normalize at 5100 Å because it is a commonly used spectral region that is less affected by complex emission or absorption features. We omitted from the calculation any radio-loud quasars or known hosts of BAL outflows. The variable composite was derived from SDSS spectra of 22 sources, while the non-variable composite represents 15 sources.

In Figure 22, we plot the composite spectra in black for variable (top panel) and non-variable (bottom panel) sources. As before, "variable" sources are those for which the X-ray count rate is inconsistent with a constant value for at least one epoch. For comparison, we have plotted a renormalized SDSS quasar composite spectrum from Vanden Berk et al. (2001) in red in each panel.

Zoom In Zoom Out Reset image size

We find at most mild differences in the composite spectra for the Hβ λ4861 line or the broad pattern of ionized Fe emission in the plotted region. However, [O iii] emission for non-variable sources is noticeably stronger than typical in three lines (4364, 4960, and 5007 Å) where [O iii] is clearly present in emission. Visual inspection of the individual spectra supports this observation, although there is much diversity.

To investigate this trend further, we estimated the monochromatic luminosity νL_ν at 5100 Å and also the equivalent width (EW) of the [O iii] λ5007 emission line for each spectrum. We used a simple linear continuum fit across the Hβ and [O iii] emission line region to estimate the EW. No attempt was made to deblend ionized Fe emission; although we note that the composite spectra do not indicate a strong difference in Fe emission between the two subsamples, in general. Figure 23 shows the results of these calculations. While both subsamples (variable and non-variable) cover a similar range of luminosities and EWs, the non-variable sources are more likely to occupy the high-EW portion of the plot, with log (EW) > 1.3. A K-S test finds no evidence of a difference in the luminosity distributions, but a suggestive probability (97% confidence) that the EW distributions differ.

Zoom In Zoom Out Reset image size

Because we have limited this study to sources with at least 50 counts per epoch on average, our classifications of "variable" and "non-variable" are very similar to placing a cut on the amplitude of variability. Still, the individual observations do have different sensitivities, so the two concepts are not exactly similar. In order to investigate [O iii] strength as a function of variability amplitude, we define

$$\begin{eqnarray} S^{\prime }_j &\equiv & \max _i(|S_{ij}|), \end{eqnarray} \tag{ 9 }$$

where S_ij is defined in Equation (5) and the maximum is taken over all epochs i for source j. S'_j represents the maximum amplitude of variation (in units of σ) for each of the 37 sources we are investigating. In fact, we find that sources with larger values of S'_j are generally the sources we flagged as variable; the two criteria are nearly equivalent. If we cut the sample at the median value S'_j = 2.8, then a K-S test indicates that the [O iii] distributions of the two subsamples are incompatible at 96% confidence. If we place the cut closer to the division between the "variable" and "non-variable" sources, at S'_j = 1.5 to 1.9, then the K-S test indicates a sample difference at ≳99.5% confidence. (We find similar results, though less significant by a few percent, using a definition in which variability amplitude is defined as ≡ max _i(|(f_ij − c_j)/c_j|), with the best-fit count rate in the denominator instead of the measurement error.) As before, we find that there are some indications of a real effect in the data, but larger samples are needed to test this result more sensitively.

It is also possible that some other effects in the data are influencing our "variable" and "non-variable" subsamples. If we had enough sources, it would be ideal to control for factors such as the number of observations of a given source and the timescales over which the sources were observed. With our small sample, this is difficult to do. Although our variability criterion is not associated with a specific timescale (but is simply compared to the best-fit constant rate), we can investigate what happens if we remove 10 cases of "variable" sources that were observed more than two times, under the extreme assumption that they would not have been flagged as variable from their first two observations alone. In this case, the split in [O iii] EWs becomes less significant (with a 19% probability of greater difference in EW distributions), but this is expected due to the fact that the remaining sample is small. Effectively, this cut has removed many "variable" sources without changing the cluster of non-variable sources having stronger-than-average [O iii] emission.

$[{\rm O\,\mathsc {iii}}]$ emission and X-ray variation are two effects that are associated with the well-known "Eigenvector 1" description of AGNs. Our current study has been partially motivated by a desire to determine if and how such factors are actually related, and what the physical significance of such a relation would be. We discuss this issue further in Section 4.4.

Chandra's ability to detect significant X-ray variation depends strongly on the number of counts obtained for a source. Figure 9 demonstrates the sensitivity of Chandra to variation in archived observations of SDSS spectroscopic quasars. A minimum of ∼10 counts per epoch are needed before variability can be detected at the 100% level (corresponding to a source doubling or halving its count rate compared to the best-fit constant rate). For sources with fewer than 20–30 counts per epoch, we see minimal evidence of variability.

Overall, variability is detected in ≈30% of the sample, including sources with low numbers of counts. For sources with larger numbers of counts, variability can be measured with greater precision. Eleven of 19 sources having >500 mean counts are detected as variable at the ≳13% level, and 7 of 10 sources having >1000 mean counts are detected as variable at the ≳10% level. Based on this small sample, we expect that X-ray variability at the >10% level would be observed in most, if not all, quasars with sufficient numbers of counts. Paolillo et al. (2004) likewise detected variability in over half of the AGNs in their analysis of the CDFS and estimated that intrinsic variability could be present in a fraction approaching 90% of their sources.

The level of fractional variation (σ(F) ≈ 16%) we observe certainly contributes to the scatter in observed X-ray-to-optical flux ratios for non-BAL SDSS quasars. This amount of variation is smaller than that required to explain the scatter in X-ray-to-optical spectral slopes as due to variability alone (Strateva et al. 2005; Steffen et al. 2006; Gibson et al. 2008a), at least over timescales of ≲2 yr. Using simultaneous UV/X-ray observations from XMM-Newton, Vagnetti et al. (2010) have also found that the observed scatter in X-ray-to-optical flux ratios cannot be explained by variability alone. Some physical differences in quasars must also contribute to an intrinsic scatter in X-ray-to-optical flux ratios.

On timescales shorter than Δt_sys ∼ 2 × 10⁵ s, we do not detect quasar variability. At longer timescales, the amount of intrinsic variability is roughly constant or perhaps slowly increasing. If we are observing a "break" in the ensemble variability properties at Δt_sys ∼ 2 × 10⁵ s, it is not clear how such a break is related to timing breaks identified in the power spectra of local Seyfert AGNs (e.g., McHardy 2010). Additional work is needed to determine how such a "break timescale" may be affected by the inclusion of quasars spanning a range of redshifts and luminosities in our sample.

Although earlier studies have suggested that variability increases at higher redshifts (Almaini et al. 2000; Manners et al. 2002; Paolillo et al. 2004), we have not been able to confirm that result in our data (Sections 3.4 and 3.5). In fact, any increase in excess variance at z > 2 would be above the ≈2σ upper limit on the excess variance we actually measure (Section 3.4). We also do not detect any increase in variation at high black hole masses. Increased variability at higher redshifts was not detected even when we cut the sample to remove sources with lower numbers of counts, thus maximizing the sensitivity to any variability. Looking forward, dedicated X-ray observations targeting higher-redshift quasars that are matched to lower-redshift quasars in other sample properties would be ideal for a definitive test of redshift dependence.

We stress that our approach has been to look for an empirical signature of variation using the data and observational bandpass that are available to us. Some earlier studies (A00, M02) have generated correction factors for their measured variance that model the effects of sampling sources on different timescales. These corrections rely on an assumed shape of a quasar power spectrum to boost the observed variation of sources observed over shorter times, as quasars tend to vary more over longer timescales. However, several factors argue against attempting such corrections to our study. It is not clear what the shape of the power spectrum should be over the wide range of timescales in our data, especially since many of our data points sample timescales beyond the power spectrum break (see Figure 7). We know even less about this region of the power spectrum for quasars than we do about the slope of the quasar power spectrum on shorter timescales. Furthermore, some statistics in our study (such as S_ij) do not have a well-defined timescale associated with them. For the case of the σ_EV statistic, which we carefully constructed to depend on only one timescale per source, we did not find evidence for increased variability at higher redshifts even when we limited the range of system-frame or observed-frame timescales (Section 3.4).

Applying a time dilation correction factor of (1 + z)^1/4 (A00) would have only a small (40%) effect even at z = 3 and would still not result in a detection of "increased variability" at higher redshifts in our data. Similarly, the study of P04 notes that "no correction for time dilation has been applied to these data since it would require us to make a somewhat arbitrary assumption about the power density spectrum (PDS) of our sources...." That assumption becomes even more arbitrary for the longer timescales in our data. There is, of course, some danger that if we mischaracterize the correction for our sample, we will be creating a redshift-dependent trend in the data. For these reasons, we prefer to leave our results in terms of observational, rather than modeled, constraints. Using the data we provide, readers can straightforwardly test the effects of modeling different correction factors.

Another topic that has received somewhat less attention is the possibility that the variability signal could decrease with redshift because different spectral regions are being shifted across the Chandra bandpass. We have performed a basic study of how spectral variation may change as a function of photon energy. In our case (Sections 3.8 and 4.3), we found some evidence that hard spectral regions may show increasing influence of a constant component that shifts at higher redshifts into an energy range where Chandra is most sensitive. At our current stage of understanding, attempting to correct for this effect would introduce systematic uncertainty, and we prefer to stay with our empirical approach.

Increased variability at higher redshifts would have intriguing physical implications, although we note that it is not obviously expected in the currently favored "cosmic downsizing" picture for SMBH growth (e.g., Cowie et al. 2003; Marconi et al. 2004; Brandt & Hasinger 2005, and references therein). In this picture, we observe multiple largely independent generations of growing SMBHs as we look back in cosmic time, with the more massive SMBHs generally growing earlier. At a fixed high luminosity, we then do not expect to see SMBH masses dropping (and thus perhaps variability strength rising) toward higher redshifts.

Quasar variation appears similar to that observed for lower-luminosity Seyfert AGNs in several respects. The X-ray spectra of quasars vary in shape as well as brightness. In general, the spectrum of a given source is flatter when that source is in a fainter state, and the spectrum steepens as the source brightens. The change in HR is anti-correlated with the change in luminosity (Figure 16), demonstrating the steepening-brightening trend. Direct spectral fits for multiple epochs of high-S/N sources (Figure 17) demonstrate a strong anti-correlation between ΔΓ and Γ that tends to move spectra toward Γ ∼ 2.

The power-law spectral index (derived from our simple fit model) may vary significantly for a given source. Based on our fits to sources with large numbers of counts, intrinsic variation by ΔΓ ∼ 0.1–0.2 is common. Some caution is therefore warranted when using the spectral slope as an indicator of quasar physical properties such as black hole mass or accretion rate, as the spectral slope (measured according to our simple model) can vary with X-ray luminosity over time. Of course, the variation mechanism may contribute to a relationship between Γ and accretion rate (L/L_Edd, where L_Edd is the Eddington rate) that has been found in single-epoch observations of quasar ensembles (e.g., Shemmer et al. 2008; Risaliti et al. 2009). For the epochs in Figure 18, we measure a best-fit relationship:

$$\begin{eqnarray} \Delta \Gamma &=& (0.534 \pm 0.104) \log (N_1 / N_0) - (0.039 \pm 0.012),\qquad \end{eqnarray} \tag{ 10 }$$

where N₁/N₀ is the ratio of power-law normalizations between the two epochs, and ΔΓ ≡ Γ₁ − Γ₀ is the change in Γ. (We note that this fit may not adequately describe the properties of individual sources, as we are simply assuming that one model describes the variation in all epochs. We do not have sufficient data to fit each source individually.) The slope in Equation (10) is a bit steeper than the slope (0.31 ± 0.01) in Equation (1) of Shemmer et al. (2008) that related Γ to L/L_Edd for individual AGNs. The overall effect may be similar to the "intrinsic" and "global" Baldwin effects (e.g., Baldwin 1977; Kinney et al. 1990; Pogge & Peterson 1992), in that the variation in an individual object may be more extreme than the overall trend observed across single-epoch measurements of a large number of sources. At the least, intrinsic variation contributes to scatter observed in the sample of sources.

Simple absorption models did not adequately describe the observed anti-correlation between spectral hardness and luminosity, although we cannot rule out the possibility of more complex, multi-parameter absorption models. However, a two-component emission model did reproduce the trend reasonably well without fine-tuning. The emission model consists of a variable power law with Γ = 2 and a hard (Γ = 1.3), constant component normalized to ∼20% of the variable power-law flux at 1 keV. Previously, Pao04 suggested that the increasing contribution of a flatter, reflected component may result in lower variability levels for sources with "absorbed," or harder, spectra. In our data, the relationship between hard- and soft-band X-ray fluxes, extrapolated to zero soft-band flux, also suggests that a constant, hard spectral component is present in our ensemble of sources. Similar models have been used to describe flux variation in Seyfert AGNs such as MCG –6–30–15 (e.g., Taylor et al. 2003; Vaughan & Fabian 2004).

If a typical quasar spectrum includes a non-negligible constant component that is flatter than the dominant component with power-law index Γ = 2, we might expect that this component changes the observed spectral shape for higher-redshift quasars. There has been disagreement over whether spectral shapes evolve with redshift, but so far strong evidence for evolution has not been discovered. Shemmer et al. (2005) find that the spectral slope for a joint fit to 10 of the most luminous quasars at 4 < z < 6.3 is Γ ≈ 2, similar to that of lower-redshift quasars. However, this result applies to the rest-frame 2–10 keV range, and not necessarily to a harder component that becomes more evident at higher energies. Vignali et al. (2005) report a slightly harder (Γ = 1.9^+0.10_{− 0.09}) observed-frame photon index for a joint fit of 48 quasars detected in X-rays at 4.0 < z < 6.3. Green et al. (2009) find no evidence for spectral evolution in a large sample of optically selected quasars out to z < 5.4. However, Bechtold et al. (2003) previously observed spectral flattening in a sample of 17 optically selected quasars at 3.7 < z < 6.3, and Kelly et al. (2009) observed marginally significant spectral flattening with redshift out to z ∼ 4.7.

To determine how a two-component model would affect the observed spectrum at high redshift, we randomly generated spectra for nine Chandra sources at z = 4.2. We simulated each observation from a model with two power laws having Γ = 2 and Γ = 1.3. The flatter power law was normalized to 20% of the steeper component at rest frame 1 keV; it starts to rise above the steeper component at rest-frame 10/(1 + z) keV. Then we fit each spectrum with a power-law model to determine typical measured photon indices, finding an average value of Γ ≈ 1.72 ± 0.17.

For comparison, we fit Galactic-absorbed power-law models to higher-energy regions of spectra for quasars at z > 4. The quasars were described in Shemmer et al. (2005), with spectra kindly provided to us by O. Shemmer. We fit the observed-frame energy range E₀–E₁ keV, where E₀ ≡ 10/(1 + z) keV and E₁ varied, in order to characterize the region where a hard component may dominate (in our simple model). We fit the spectra of all three EPIC cameras simultaneously, requiring them to have the same photon index, but allowing slightly different normalizations to account for cross-calibration differences between instruments. In all cases, the photon indices for the observed-frame spectra were flatter than for the rest-frame 2–10 keV region. For a fit range extending up to E₁ = 5 keV in the observed frame, the spectrum flattened by ΔΓ = 0.1 to 0.9. For higher values of E₁, the spectrum appeared to flatten more, although background contamination becomes higher at these energies as well. This result suggests that a hard spectral component is becoming more evident at high energies. Although inspection indicates that backgrounds are generally well below count rates in our fit region, it would be useful to obtain larger, high-quality samples to examine the shape of the high-energy quasar spectrum more reliably.

In addition to the hypothesis of a constant harder component in quasar spectra, the assumption (based on studies of local Seyfert AGNs) that the relationship between soft and hard X-ray count rates can be extrapolated linearly should also be tested with larger, high-quality data sets. The two-power-law model we use here is the simplest attempt to explain the putative constant component, but more sophisticated models may be required that cut off at higher energies if the flatter component is not evident in high-redshift spectra. Our results also suggest that studies of high-redshift populations should test multi-component models against individual sources in their data to constrain spectral complexity beyond a single power-law shape. The spectral shape measured from a stack or joint fit can be biased toward the brighter sources in a sample, so care should ideally be taken to account for selection effects and brightening-steepening behavior as sources vary.

In our sample, sources classified as X-ray non-variable have [O iii] λ5007 Å emission lines that are about 12 Å stronger (on average) than those of variable sources (Section 3.10). This was the most prominent effect we found when comparing composite spectra of the variable and non-variable AGNs. In the optical bandpass, variability has been observed to be related to both [O iii] and H_β emission. Quasars with stronger [O iii] emission have larger optical variability amplitudes (Giveon et al. 1999; Mao et al. 2009), although with significant scatter. The cause of this relation is not known; it has been suggested that the effect may be due to an enhanced ionization rate from a variable ionizing continuum (Giveon et al. 1999), or perhaps associated with accretion instabilities or star formation (Mao et al. 2009). However, in the X-rays, we see the opposite effect: stronger [O iii] emission is associated with lower levels of X-ray variability.

The effect may be anomalous. The set of sources for which we can sensitively measure [O iii] emission and X-ray variability is relatively small, and the effect is not highly significant. However, we do not see a significant difference (according to a K-S test) between the distributions of monochromatic luminosity (νL_ν at 5100 Å), X-ray count-rate luminosity, or mean numbers of X-ray counts for the variable and non-variable sources. There may be other factors differentiating the two subsamples. For example, some absorption mechanisms may affect both [O iii] strength and also the X-ray spectrum (e.g., Caccianiga & Severgnini 2011), weakening the variability signal. Finally, the effect may be a random deviation in our small sample.

In our search for short-term variability on kilosecond timescales (Section 3.9), we identified an unusual quasar, SDSS J1238, that was flagged as variable in multiple epochs. This quasar is NLS1-like, with atypically weak [O iii] λ5007 Å emission, a H_βλ4862 FWHM ≈ 1860 km s⁻¹, and stronger ionized Fe emission compared to the average SDSS quasar. A tendency for rapid, high-amplitude X-ray variation is a well-known property of NLS1 AGNs (e.g., Boller et al. 1996; Pogge 2000; Komossa 2007).

With EWs from a few to ∼100 Å, the [O iii] λ5007 Å lines in our study span the range of EWs originally used by Boroson & Green (1992) to define "Eigenvector 1." This empirically based collection of properties includes an anti-correlation between Fe ii and [O iii] emission strengths; it also includes trends with radio loudness and some H_β line features such as FWHM and asymmetry. NLS1 AGNs fall at one end of this eigenvector, generally having narrow H_β lines, weak [O iii] emission, and strong ionized Fe emission. The quasar J1238 also shows similar properties when compared to typical SDSS quasars. In our sample, X-ray "non-variable" quasars may exemplify the (radio-quiet) population at the opposite (strong-[O iii]) end of Eigenvector 1, while quasars showing moderate X-ray variability have [O iii] emission levels more typical of ordinary SDSS quasars.¹³ The sample of strong-[O iii] emitters that have been observed multiple times in X-rays should be expanded in order to test the potential relation between X-ray variation and Eigenvector 1 more sensitively. Some new sources can be added from the archives as X-ray and optical surveys progress, but targeted observations would be most effective to obtain a significant ensemble of the strongest-[O iii] emitters.

A link between X-ray temporal properties and [O iii] emission would be interesting because it would represent another physical connection between small-scale (disk corona) and large-scale (NLR) AGN physics (e.g., Brandt & Boller 1998). (Of course, with larger data sets we could also test more sensitively whether additional parameters, such as Balmer line profiles, are related to X-ray variation.) In one category of models, a third parameter modulates both X-ray and NLR emission, creating an indirect relation between the two. For example, [O iii] emission strength is believed to be dominated by geometric factors (e.g., Baskin & Laor 2005). The Eddington ratio (L/L_Edd) of bolometric to Eddington luminosity could control both the obscuration of X-ray/UV radiation essential for NLR ionization (e.g., Abramowicz et al. 1980; Boroson & Green 1992; Chang et al. 2007) and also X-ray emission properties. Alternatively, some physical effects could create a more direct connection between the corona and NLR. Very large NLRs that produce the strongest [O iii] emission lines may require replenishment (Netzer et al. 2004), and small-scale jets have been observed to interact with the NLRs of radio-quiet Seyfert AGNs (e.g., Falcke et al. 1998; Ho & Peng 2001; Leipski et al. 2006; Wang et al. 2011). Jet activity has been associated with other Eigenvector 1 properties, including radio loudness and X-ray spectral slope, and radio loudness (indicative of jet activity) is associated with decreased variability in our subsample of radio-loud quasars (Section 3.3.1).