The Extreme Ultraviolet Variability of Quasars
We study the extreme ultraviolet (EUV) variability (rest frame wavelengths 500 - 920 ) of high luminosity quasars using HST (low to intermediate redshift sample) and SDSS (high redshift sample) archives. The combined HST and SDSS data indicates a much more pronounced variability when the sampling time between observations in the quasar rest frame is sec compared to sec. Based on an excess variance analysis, for time intervals sec in the quasar rest frame, of the quasars (4/40) show evidence of EUV variability. Similarly, for time intervals sec in the quasar rest frame, of the quasars (21/38) show evidence of EUV variability. The propensity for variability does not show any statistically significant change between sec and sec (1 yr). The temporal behavior is one of a threshold time interval for significant variability as opposed to a gradual increase on these time scales. A threshold time scale can indicate a characteristic spatial dimension of the EUV region. We explore this concept in the context of the slim disk models of accretion. We find that for rapidly spinning black holes, the radial infall time to the plunge region of the optically thin surface layer of the slim disk that is responsible for the preponderance of the EUV flux emission (primarily within 0 - 7 black hole radii from the inner edge of the disk) is consistent with the empirically determined variability time scale.
The spectral energy distribution (SED) of quasars typically peaks around 1100Å then decays in the extreme ultraviolet (EUV) shortward of this peak (Telfer et al., 2002; Stevans et al., 2014). The broadband optical/UV emission is widely believed to be optically thick emission from gas accreting onto a central supermassive black hole (Sun and Malkan, 1989). The emission naturally arises in an inflow from the differential shearing of the infalling gas (Lynden-Bell and Rees, 1971). For orbital velocities on the order of the local Keplerian angular velocity, the angular velocity increases towards the black hole with an ever increasing gradient. Thus, the maximum dissipation and the highest frequency thermal emission arises from the innermost regions of the accretion flow (Shakura and Sunyaev, 1973). The EUV emission just shortward of the SED peak (rest frame wavelengths 500 - 1000 ) should be predominantly from the innermost optically thick accretion flow (Sun and Malkan, 1989; Szuszkiewicz et al., 1996; Zhu et al., 2012). In this paper, we study the time variability behavior of the EUV continuum emission. We also speculate on the size of the EUV emitting region based on slim accretion disk models and our temporal analysis.
The EUV continuum is a uniquely important region of the quasar spectrum for two reasons. Firstly, an estimate of the size of the EUV region can be compared to the dimension expected from accretion disk theory in order to critically analyze the qualitative applicability of accretion disk theory to quasar accretion flows near the central supermassive black hole. In particular, as noted above, being the highest energy optically thick emission should place this region near the inner edge of the accretion disk. If accretion disk models are appropriate then an estimate of the location of the EUV emitting region should indicate that the EUV emitting region lies close to the inner edge of the disk. The EUV region is a modest fraction of the total optically thick accretion flow luminosity, (Punsly, 2014). This provides a second constraint on accretion disk models. The small total EUV luminosity seems to favor a dimension of the EUV emitting region comparable to those of the inner edge of the disk as opposed to a dimension that is an order of magnitude or more larger than the inner edge of the disk (Sun and Malkan, 1989; Zhu et al., 2012). We will critically evaluate our size estimates of the EUV region derived in this article in the context of slim accretion disk theory. The second unique property of the EUV continuum is a consequence of its connection to the launching of powerful relativistic jets in quasars. There exists a strong correlation between jet power and EUV spectral index that has been recently demonstrated in radio loud quasars (Punsly, 2014, 2015b; Punsly et al., 2016). In particular, the long term time average jet power, , is correlated with the spectral index in the EUV, ; defined in terms of the flux density by computed between 700Å and 1100Å . Larger tends to decrease the EUV emission. In general, lobe dominated radio loud quasars tend to have larger values of than radio quiet quasars. The straightforward implication is that the EUV emitting region is related to the jet launching region in quasars. Therefore, an estimate of the size and location of the optically thick EUV emitting region yields a constraint on the size and location of the jet launching region. These two aspects indicate the potential power of estimating the size of the EUV emitting region.
The dimensions of the emitting for the EUV continuum region is many orders of magnitude smaller than the highest spatial resolution of optical telescopes. Consequently, observations cannot directly constrain its size. However, time variability can be used as a crude probe of the size of the region. We use multi-epoch archival UV observations with the Hubble Space Telescope (HST) of quasars with a redshift and Sloan Digital sky Survey (SDSS) of high redshift quasars in order to find the EUV variability time scales in the quasar rest frame. Our sample is restricted to time intervals of less than one year in the quasar rest frame. The combined data-set of HST and SDSS observations presented here shows that a dramatic change in the behavior of variability occurs on times scales less than one year and the sole goal of this effort is to achieve the maximum temporal resolution of this region. As such, we do not investigate duty cycles (i.e., multiple variations over a few years that can approximately cancel out the total variation) or power spectra of the variability. It is not that such studies are not interesting in their own right, but our analysis yields a simple clean result independent of the interpretive filter of any analysis tool. Our sample of data is small for the purposes of a time series analysis and does not possess a uniform sampling of time intervals. Power spectra are often difficult to interpret when based on small number statistics and even more so with sparse and irregular temporal sampling (Kittel, 1961; Neira and Constantinides, 1996). Consequently there is not much motivation to study time intervals longer than one year in the quasar rest frame for this study.
Previous studies of quasar variability have focused on the UV and optical bands in the quasar rest frame. In Wilhite et al. (2005) composite spectra from - in the low and high phases were compared for a sample of 315 highly variable SDSS quasars. The study showed that the blue end of the spectra were more variable. However, the only variability data as a function of time for the entire SDSS sample of objects was based on a broad band flux integration. Because of the broad line contamination, this type of analysis does not extract the time dependent variability of the continuum - our primary goal in the EUV. SDSS photometry and Palomar Sky Survey photometry were compared in order to study the time dependent variability in quasars in the UV and optical bands (Macleod et al., 2012). Broad band photometry was also implemented in Welsh et al. (2011) using GALEX data in order to explore far UV and near UV variability in quasars. Broad band photometry suffers from broad line contamination so one does not study the continuum variability explicitly. There was an International Ultraviolet Explorer variability study of 21 quasars based on broadband photometry (diClemente 1995). This was typically far UV and near UV data, but objects had spectra that covered the EUV. Once these noisy spectra were integrated over wavelength there was large contamination from broad emission lines so this small study is not of much value for our purposes. Thus, this is the first study of the variability of the continuum in quasars on the high frequency side of the peak of the spectral energy distribution.
In Section 2, we describe the HST observations with a particular emphasis on the criteria for an accurate flux calibration. Section 3 is a similar analysis of the SDSS archives. The next section combines the two samples and analyzes the data scatter. In Section 5, we compare our results to slim disk accretion flow models.
2 HST Observations
Hubble Space Telescope (HST) observations can reliably detect the quasar EUV continuum below the Lyman continuum if . In principle, HST observations can be used to monitor the EUV continuum variability of quasars over the last quarter century. There are two issues that limit the use of the HST for EUV monitoring of quasars. Firstly, the long exposures necessary with HST for individual quasar UV spectra make the time allocations required for long term UV monitoring prohibitively large. Thus, monitoring of the quasar UV spectra consumes too much of the the available resource and is not a viable proposal topic. Therefore, variability can only be found serendipitously. Secondly, in order to detect variability, reliable flux calibration is required. This is a major concern for narrow apertures in general and wider apertures if the source is not accurately aligned with the slit center.
There were many quasar observations with the Faint Object Spectrograph (FOS). However, the alignment issue was particularly severe in FOS pre-COSTAR observations (see FOS Instrument Science Reports CAL/FOS-107 and CAL/FOS-122). The two main target acquisition methods used for FOS are ACQ/BINARY and ACQ/PEAK. Often, in order to save observing time, the target acquisition sequence only included ACQ/BINARY which tries to locate the object based on the input coordinates. However, the pre-COSTAR alignment was not adequate for such an expedience. Thus, many of the observations are poorly centered and seemed to indicate variability of the quasar, but in fact one is only seeing the results of instrumental error. The most reliable target acquisition sequence involves ACQ/BINARY followed by ACQ/PEAK. The routine ACQ/PEAK is an iterative process designed to position the target at the center position of the brightest grid point of the detector diode array. The nominal centering error for ACQ/BINARY alone (not followed by ACQ/PEAK) is and the nominal centering uncertainty for a ACQ/BINARY followed by ACQ/PEAK acquisition is 0.05” (Evans and Koratkar, 2004). If ACQ/BINARY was followed by ACQ/PEAK, this acquisition sequence has another smaller centering issue. ACQ/PEAK was generally not performed with the same aperture that the measurement was made in and for pre-COSTAR observations these misalignments were not well known. Due to the alignment uncertainties, Pre-COSTAR, the smallest aperture that is likely to yield a reliable flux measurement is 1.0” if a ACQ/BINARY followed by ACQ/PEAK acquisition occurs. If the the acquisition sequence is simply ACQ/BINARY then the minimum aperture is 4.3” for a reliable flux measurement. Pre-COSTAR observations in the 4.3” aperture flux measurements will typically be reliable to level from the point of view of aperture correction errors due to target mis-centering. For smaller apertures the centering uncertainty combined with the uncertainty in aperture alignments produces an absolute flux uncertainty that is too large for their consideration in this treatment if there is no ACQ/PEAK step in the acquisition sequence. If the acquisition and the observation were both in the same aperture some of this uncertainty is removed for smaller apertures (such as 1.0”). However, this instance never occurred in the observations that we considered. We looked at the acquisition files for flags that the flux measurement was suspect. Sometimes there are two spectra taken with two different gratings during the same observation that overlap in frequency. If the flux mismatch between the gratings in this overlap is less than 10% this is suggestive of a reliable flux calibration even for an observation with no ACQ/PEAK. We excluded observations with the 0.25” x 2.0” slit in all instances because of the impact of poor target centering. The efficiency of the detectors is very low at the extreme blue end of the G130H grating coverage, and there can be significant shot noise that distorts the true continuum at the shortest wavelengths (Ian Evans private communication 2015). Thus, the blue end of the G130H grating was not considered suitable for accurate absolute flux measurements.
FOS post-COSTAR aperture positions and alignments were known much more precisely. An aperture larger than 0.3” with a ACQ/PEAK during the target acquisition should yield an accurate flux measurement. If no ACQ/PEAK was performed, the 0.12” centering uncertainty equates to a 2% flux uncertainty in a 1.0” circular aperture (FOS Instrument Science Reports CAL/FOS-140) and this is the minimum aperture that is chosen for reliable flux measurements. Since the alignments were well known, post-COSTAR, images (if they exist) can also be used to determine the centering of the source in the aperture. At the end of the acquisition sequence some observers took an ACQ image. However, we were unable to verify accurate centering for any of the potentially relevant observations with this method when no ACQ/PEAK step was performed. We were only able to find a few quasar FOS observations with an adequate aperture size in both the first and second epoch separated by less than one year in the quasar rest frame. We also found one second epoch observation with the HRS spectrograph and a 2.0” aperture (Post-COSTAR) for PG 1008+133. According to GHRS Instrument Science Report No. 79, the BRIGHT=RETURN (RTB) feature that was chosen for the spiral search will automatically return the spacecraft pointing to the dwell position with the most counts in the aperture. This combined with the 2.0” aperture should minimize slit losses. Based on this centering, the uncertainty in the absolute flux should be at most 10% - 15% (Heap et al., 1995; Soderblom and Sherbert, 1997).
After the initial experimental stage of testing the COS (Cosmic Origins Spectrograph), that ended in late September 2009, the flux calibrations are reliable. Calibration routines are updated based on the change in sensitivity of the instrument (eg. see Instrument Science Report COS 2011-02), thus systematic errors are minimized. The 2.0” slit is suitable for proper alignment as well. Thus, these observations represent the bulk of our HST sample. The Space Telescope Imaging Spectrograph (STIS) generally implemented narrow apertures in order to resolve absorption features and these types of observations were not intended for accurate flux calibrations. Slight offsets in the centering of the object would result in loss of light that would mimic variability. Thus, in spite of numerous potential observations, STIS was deemed unsuitable for this work for the potential observations of quasars that we found.
|Source||z||Date||TaaIn quasar rest frame||Instrumental Setup|
|Epoch 1||Epoch 2||Epoch 1||Epoch 2|
|PG 0117+213||1.490||07/09/1993||12/15/1994||1.77||FOS G160Lbb4.3” aperture Pre-COSTAR with peak up andPost-COSTAR 0.3” aperture with peak up.||4.3||FOS G270Hbb4.3” aperture Pre-COSTAR with peak up andPost-COSTAR 0.3” aperture with peak up.||0.3|
|UM 675||2.140||10/28/1990||02/03/1992||1.25||FOS G160L||4.3||FOS G270H||4.3|
|FIRST J020930.7-043826||1.130||12/23/2010||1/6/2011||0.06||COS G130M||PSA||COS G230L||PSA|
|SDSSJ 022614.46+001529.7||0.620||9/22/2010||7/29/2011||1.65||COS G130M||PSA||COS G130M||PSA|
|PKSB 0232-042||1.440||1/9/2010||2/16/2010||0.13||COS G130M||PSA||COS G160M||PSA|
|PKSB 0232-042||1.440||7/31/2015||8/8//2015||0.03||COS E225M||PSA||COS E225M||PSA|
|LBQS 0302-0019||3.290||11/4/1995||12/14/1995||0.06||HRS G140L||2.0||HRS G140L||2.0|
|PG 1008+133||1.300||11/3/1993||6/3/1995||2.14||FOS G160Lee4.3” aperture Pre-COSTAR with peak up and Post-COSTAR 8/1994 2.0” aperture with spiral search for peak flux||4.3||HRS G140Lee4.3” aperture Pre-COSTAR with peak up and Post-COSTAR 8/1994 2.0” aperture with spiral search for peak flux||2.0|
|PG 1115+080A||1.740||3/23/2012||2/12/2014||2.13||COS G130M||PSA||COS G140L||PSA|
|3C 263||0.650||1/2/2010||2/22/2010||0.28||COS G130M||PSA||COS G160M||PSA|
|PG 1148+549||0.980||12/26/2009||12/30/2009||0.02||COS G130M||PSA||COS G130M||PSA|
|PG 1206+459||1.160||12/29/2009||1/5/2010||0.04||COS G160M||PSA||COS G160M||PSA|
|PG 1206+459||1.160||1/5/2010||6/14/2010||0.63||COS G160M||PSA||COS G160M||PSA|
|PG 1329+412||1.940||08/17/1995||06/02/1996||0.84||FOS G190HccNo peak up, 1.0” aperture Post-COSTAR, both epochs. The 1995 observation showed 2.6% flux variation between sub-exposures in the sampled range 2230 Å –2280 Å.||1.0||FOS G270Hbb4.3” aperture Pre-COSTAR with peak up andPost-COSTAR 0.3” aperture with peak up.||1.0|
|1334-005||2.820||02/06/1993||06/02/1995||1.89||FOS G160LddNo peak up first epoch, but a 4.3” aperture. Peakup second epoch, 1.0” aperture Post-COSTAR.||4.3||FOS G190HddNo peak up first epoch, but a 4.3” aperture. Peakup second epoch, 1.0” aperture Post-COSTAR.||1.0|
|PG 1338+416||1.220||5/24/2010||5/30/2010||0.02||COS G130M||PSA||COS G160M||PSA|
|PG 1407+265||0.940||2/12/2010||2/21/2010||0.04||COS G130M||PSA||COS G130M||PSA|
|LBQS 1435-0134||1.310||8/8/2010||8/22/2010||0.06||COS G130M||PSA||COS G160M||PSA|
|SDSSJ 161916.54+334238.4||0.470||8/19/2010||6/17/2010||1.76||COS G130M||PSA||COS G140L||PSA|
|PG 1630+377||1.480||12/13/2009||8/01/2010||0.78||COS G130M||PSA||COS G160M||PSA|
|PG 1630+377||1.480||8/1/2010||11/26/2010||0.41||COS G160M||PSA||COS G160M||PSA|
|HS 1700+6416||2.740||07/26/1994||5/29/1996||1.52||HRS G140L||2.0||HRS G140L||2.0|
|HS 1700+6416||2.740||10/10/2011||4/30/2014||2.11||COS G140L||PSA||COS G140L||PSA|
|HS 1700+6416||2.740||4/29/2015||5/22/2015||0.05||COS G130M||PSA||COS G130M||PSA|
|Source||T||Infall RadiusaaThe infall time is equated to a radius based on the slim disk models of Sadowski (2011) that are discussed in Section 5. Assumes a black hole spin of||VariabilitybbAs computed per Equation(4).|
|PG 0117+213||9.77ccMass estimated from H profile in Shields et al. (2003), using the formula of Assef et al. (2011).||0.47||388||3||0.20|
|UM 675||6.95ddMass estimated from H profile in Dietrich et al. (2009), using the formula of Assef et al. (2011). was computed using Equation (3) and the spectrum in Steidel and Sargent (1991).||0.58||390||3.35||0.25|
|FIRST J020930.7-043826||1.41eeMass estimated from MgII in Finn et al. (2014).||0.52||92||2.3||0.05|
|SDSSJ 022614.46+001529.7||0.30ffMass estimated from MgII profile in SDSS spectrum, using the formula of Trakhtenbrot & Netzer (2012).||0.68||11670||14.1||0.10|
|PKSB 0232-042||1.48ggMass estimated from MgII profile in 6df survey spectrum, using the formula of Trakhtenbrot & Netzer (2012).||1.72||192||6.5||0.05|
|PKSB 0232-042||1.48ggMass estimated from MgII profile in 6df survey spectrum, using the formula of Trakhtenbrot & Netzer (2012).||1.72||46||3.35||0.00|
|LBQS 0302-0019||2.98hhMass estimated from H profile in Syphers and Shull (2014), using the formula of Assef et al. (2011). was computed using Equation (3) and the SDSS spectrum. see Table 2.||2.02||58||3.85||0.05|
|PG 1008+133||2.40ffMass estimated from MgII profile in SDSS spectrum, using the formula of Trakhtenbrot & Netzer (2012).||1.02||1913||9.9||0.50|
|PG 1115+080A||4.35iiSee the detailed discussion in the text||0.64||1095||5.5||1.15|
|3C 263||1.00jjMass estimated from H profile in Marziani et al. (2003a), using the formula of Assef et al. (2011).||0.83||591||5.75||0.18|
|PG 1148+549||1.86ffMass estimated from MgII profile in SDSS spectrum, using the formula of Trakhtenbrot & Netzer (2012).||0.61||24||2.35||0.00|
|PG 1206+459||4.79ffMass estimated from MgII profile in SDSS spectrum, using the formula of Trakhtenbrot & Netzer (2012).||0.55||16||2.15||0.00|
|PG 1206+459||4.79ffMass estimated from MgII profile in SDSS spectrum, using the formula of Trakhtenbrot & Netzer (2012).||0.55||282||2.85||0.33|
|PG 1329+412||2.19ffMass estimated from MgII profile in SDSS spectrum, using the formula of Trakhtenbrot & Netzer (2012).||1.36||820||8.6||0.10|
|1334-005||3.06kkMass estimated from the SDSS CIV profile using the formula of Park et al. (2013).||0.89||1320||7.1||0.15|
|PG 1338+416||1.66ffMass estimated from MgII profile in SDSS spectrum, using the formula of Trakhtenbrot & Netzer (2012).||0.72||30||2.2||0.00|
|PG 1407+265||1.41eeMass estimated from MgII in Finn et al. (2014).||1.39||61||4.65||0.00|
|LBQS 1435-0134||6.03ffMass estimated from MgII profile in SDSS spectrum, using the formula of Trakhtenbrot & Netzer (2012).||0.61||20||2.2||0.00|
|SDSSJ 161916.54+334238.4||0.83ffMass estimated from MgII profile in SDSS spectrum, using the formula of Trakhtenbrot & Netzer (2012).||0.22||4540||4.45||0.20|
|PG 1630+377||3.63ffMass estimated from MgII profile in SDSS spectrum, using the formula of Trakhtenbrot & Netzer (2012).||1.15||463||5.8||0.00|
|PG 1630+377||3.63ffMass estimated from MgII profile in SDSS spectrum, using the formula of Trakhtenbrot & Netzer (2012).||1.15||241||4.6||0.00|
|HS 1700+6416||7.94kkMass estimated from the SDSS CIV profile using the formula of Park et al. (2013).||1.26||411||5.8||0.05|
|HS 1700+6416||7.94kkMass estimated from the SDSS CIV profile using the formula of Park et al. (2013).||1.26||570||6.7||0.00|
|HS 1700+6416||7.94kkMass estimated from the SDSS CIV profile using the formula of Park et al. (2013).||1.26||17||2.4||0.00|
|Source||z||Date||T bbTime interval in the quasar rest frame||Central||Flux 25 Å||Variability||Infall|
|SDSSJ||Epoch 1||Epoch 2||Epoch 1||Epoch 2||Radius|
|[MJD]||[MJD]||[ s]||[Å]||[10 erg s cm]||[10 erg s cm]||||[M]|
|030449.84-000813.4 aaLBQS 0302-0019. the mass estimate is from Table 2.||3.288||51817||51873||0.11||905.0||61929||54288||2.98||0.14||2.02||4.02|
|030449.84-000813.4aaLBQS 0302-0019. the mass estimate is from Table 2.||3.288||51914||53378||2.95||905.0||58113.6||40073||2.98||0.45||2.02||4.26|
|030449.84-000813.4aaLBQS 0302-0019. the mass estimate is from Table 2.||3.288||51873||51914||0.08||905.0||54288||58113.6||2.98||0.07||2.02||15.10|
In Table 1 we list all of the pairs of epochs of quasar observations (separated by less than one year in the quasar rest frame) that we believe are likely to have good calibrations. We also arbitrarily pick a minimum time interval of 2 days in the quasar rest frame in order to clearly distinguish separate observations (some observations span many orbits). We also note that the sample in this section and the SDSS sample in the next section might have three observations for an individual quasar. We always pair the observational epochs consecutively in time and exclude pairing up the first and third epochs. This procedure ensures the stochastic independence of the data for subsequent statistical analysis. The first and third pairing is not independent of the other two consecutive pairings. Also, a first and third pairing of epochs inordinately weights the intrinsic variability properties of a single quasar within our sample. In practice, adding the first and third pairings to our sample does not affect the analysis and scatter plots to follow primarily because the one year cutoff prevents long term duty cycle behavior from canceling out the variability. Excluding the first and third pairings is formally the best procedure from a statistical point of view.
The data is split between Tables 1 and 2. Table 1 is organized as follows. Columns (1) and (2) identify the quasar and the redshift. Columns (3) and (4) are the dates of the two observations. The next column is this time difference, , in the rest frame of the quasar in units of s. Column (6) describes the spectrograph and grating used in the first epoch. The aperture used in the observation is noted in the next column. Columns (8) and (9) contain the corresponding information for the second epoch. Table 2 describes the variability time scale in terms of the estimated central black hole mass. Column (1) is the same as Table 1. The next column relates this time interval to the size of the central black hole. We used virial estimates based on the full width at half maximum (FWHM) in order estimate the mass of the central black hole, in column (2). The preferred lines are the low ionization lines MgII or H. These lines were available for all but two high redshift quasars. For these objects we used a virial mass estimate based on the CIV broad emission line. Each individual measurement has about 0.3 dex uncertainty (per the references in Table 2). Thus, each estimate individually is not accurate enough for our purposes. It is generally believed that by achieving a sufficient sample size (as in this section and the next in combination) that the large scatter will be imprinted on a detectable backdrop of underlying physical trends associated with black hole mass. By considering the light travel time across the black hole mass, , in geometrized units, we were able to convert the time interval in column (5) of Table 1 into a time in geometrized units of in column (4) of Table 2. In column (3) of Table 2, we compute the Eddington ratio, ,
where is the Eddington luminosity associated with the central black hole mass and is the bolometric luminosity of the accretion flow. From Punsly (2014), the luminosity near the peak of the spectral energy distribution at Å (quasar rest frame wavelength), provides a robust estimator of the bolometric luminosity associated with the thermal accretion flow, ,
where is the flux density in the quasar rest frame evaluated at . Note that his estimator does not include reprocessed radiation in the infrared from distant molecular clouds. This would be double counting the thermal accretion emission that is reprocessed at mid-latitudes (Davis and Laor, 2011). If the molecular clouds were not present, this radiation would be directed away from the line of sight to Earth. However, it is reradiated back into the line of sight towards Earth and this combines with the radiation that has a direct line of sight to Earth from the thermal accretion flow. As such, it would skew our estimate of the bolometric luminosity of the thermal accretion flow and needs to be subtracted from the broadband spectral energy distribution. Because some of these quasars are at high redshift, there might be Lyman absorption at Å so we compute a similar bolometric correction at Å that follows from the composite SED in Telfer et al. (2002)
If possible, we use the average of the two numbers. However, if the value derived from the spectrum near Å is considerably smaller, we simply used the latter estimate. This estimate of is used in Section 5 in order to model the variability in terms of slim disk models. In particular, Column (5) gives the radial coordinate that is consistent with the infall of gas to the black hole (in the appropriate slim disk model as discussed in detail in Section 5 and Figure 8) in a time given by column (4). The last column is the relative EUV continuum variability between epochs 1 and 2,
where and are the EUV fluxes in the first and second epochs, respectively. For the sources observed with HST, the variability is computed from the flux density of the EUV continuum estimated from the entire available spectral overlap between the two epochs. We do restrict this to below in order to avoid the confusion imposed by the very strong broad lines (as discussed in the next section on SDSS spectra where it is more of an issue). The process for the HST spectra is to estimate the continuum level in each HST spectrum individually. These two continuum levels are then compared in the wavelength region of spectral overlap (below ). We also do not consider noisy edges of the spectra. The wavelength range used for the the computation of F1 and F2 varies from epoch to epoch, it depends on the redshift, the gratings that are used and the noise level at the edge of the spectra. In practice, below , the variability of the continuum can be extricated from the emission line variability in the HST spectra (see for example Figure 1). In the next section, we consider high redshift SDSS spectra for which there is strong Lyman absorption and there is no way to fit the EUV continuum from spectral data. In this case, we use photometry in a wavelength band that is relatively free of broad emission lines, - . The actual wavelength used is listed in Table 3, since small shifts in wavelength can greatly affect the integrated flux due to very deep absorption.
Figure 1 shows the COS spectra of two of the most variable quasars in the HST sample, PG 1115+080 and PG 1206+459. PG 1115+080 is the most variable quasar in our combined sample and requires a detailed discussion. PG 1115+080 is a gravitationally lensed quasar. There are four images, the weaker components B and C and a bright component A, that resolves into components A1 and A2. A1 is the located 0.48” from A2 and constitutes 56% of the total emission at (Tsvetkova et al., 2010). The PSA aperture used with these COS observations has a 2” resolution and cannot resolve component A. According to Instrument Science Report COS 2010-14, the target acquisition ACQ/IMAGE should be able to place the peak brightness associated with the component A in the center of the slit to within 0.1” with very high reliability and there is unlikely significant light loss. This implies that there is neither confusion over which lens component is detected nor whether component A was properly centered in both epochs. It is concluded that the large variability in Table 2 is real.
The second concern with the lensing is possible magnification that can occur and its effects on estimates and . To minimize this, we use the broad emission line region to make these estimates. Being much larger than the accretion disk, we expect that it is less prone to amplification. From the quasar composite spectrum in Telfer et al. (2002), one can use the composite equivalent widths of the broad emission lines in combination with Equations (2) and (3) and the composite continuum spectral index in order to achieve a bolometric correction based on line strengths. A straightforward calculation provides two approximate bolometric corrections, one based on the MgII line strength and one based on the CIV line strength, and , respectively
We do not know the relative fraction of the broad emission line luminosity that originates in component A1 of the lensed quasar, PG 1115+080. However, as stated above, A1 produces 56% of the total emission detected at Earth at (Tsvetkova et al., 2010). We assume, in the absence of any other information that of the MgII and CIV line strength is also from A1. Thus, from the estimator in Equation (5), and for MgII and CIV, respectively. Equation (5) yields reasonable agreement between the two estimators based on the two different line strengths in the case of PG 1115+080. Taking the average, we estimate that component A1 has and this represents of the quasar. Similarly, we estimate using the line strength and not the continuum with the expectation that were are minimizeing the affects of magnification from the lens. We use the more reliable low ionization line MgII and the formulas from Shen and Liu (2012) to find
Alternatively, the formula of Trakhtenbrot & Netzer (2012) yields a different estimate
Taking the average of Equations (6) and (7), we arrive at a mass estimate, . The associated Eddington ratio of 0.64 in Table 2 is typical for broad absorption line quasars such as PG 1115+080 (Ganguly et al., 2007). Any magnification of the broad line region from lensing would imply that the intrinsic Eddington ratio is actually smaller. Thus, we expect that the magnification is at most a factor of 2 or 3 (yielding an Eddington ratio of 0.35 - 0.4) in order for PG 1115+080 to be consistent with typical Eddington ratios found in broad absorption line quasars. We acknowledge that there is additional uncertainty associated with the and for PG 1115+080 due to the gravitational lensing. Fortunately, the results of this paper do not depend on these estimates for PG 1115+080. The trends persist even without this quasar. We choose to keep it in our sample since the discussion above indicates that the large variability is not an artifact of the gravitational lensing and a centering uncertainty in the apertures.
The sample is composed predominantly of radio quiet quasars. There are four radio loud quasars. Two of them are extended quasars with large steep spectrum radio lobes, 3C 263 and PKSB 0232-042. The EUV spectra were analyzed by us previously (Punsly, 2015b). The other two are the unresolved flat spectrum radio sources, FIRST J020930.7-043826 and UM 675. Considering the strong EUV emission lines present in the HST spectrum, we do not expect significant blazar contamination of the continuum in FIRST J020930.7-043826 (Finn et al., 2014). The spectral energy density of UM 675 is increasing slightly across the near UV with strong broad lines with a spectral energy density two orders of magnitude higher than the microwave luminosity and is therefore likely dominated by the accretion disk in the EUV (Steidel and Sargent, 1991). In all four cases, the variability in the EUV spectra should be indicative of the thermal accretion flow.
3 SDSS Extreme Ultraviolet Variability Data
The SDSS archive offers a large database for quasar spectra and therefore great potential for a study of EUV variability. The main obstacle results from the large redshifts that are required to shift the quasar EUV continuum into the optical window. Even a moderate sensitivity study (say a threshold of 20% change as a detection for the sake of argument) requires decent signal to noise spectra. High redshift quasars have diminished fluxes due to the large distances, but also the intervening Ly absorbing clouds greatly attenuates the signal in the EUV. The other issue that exacerbates the situation is the desire to monitor the EUV continuum and not the more distant emission line region. Thus, we must separate out the broad emission lines from the continuum. The long wavelength end of the EUV continuum, with the least Ly absorption and highest fluxes in general, is contaminated by strong, extremely broad emission lines, OVI, Ly and the blend of Ly, SVI, CIII and NIII. In order to sample the EUV shortward of these strong emission lines requires sampling the rest frame spectrum below . In the range of to there is a region with no strong emission lines. Due to the Ly absorption troughs, integration over a finite window is necessary in order to reliably compare fluxes from epoch to epoch. In order to accommodate a finite window of integration, the shortest useable wavelength for early SDSS spectra is (where is the observed wavelength). In order to place, the quasar rest frame wavelength at , requires a redshift, . The challenge is the following. As the wavelength decreases below , the Ly absorption from intervening gas, more often than not, increases as the wavelength decreases. This decreases the signal to noise of the spectrum at the shortest wavelengths. In most cases, the signal to noise ratio is too low for accurate variability studies. Thus, we attempt to find the shortest wavelength that has adequate signal to noise with the constraint that the observed wavelength corresponds to . Alternatively stated, the quasars must be bright and the Ly absorption not particularly deleterious. This eliminates most high redshift quasars.
Thusly motivated, our search criteria was all quasars in the redshift range with two epochs of observation that are separated by less than one year in the quasar rest frame. We institute an “integration window” of , which corresponds to . The window minimizes any affects caused by wavelength calibration and sky subtraction differences between observations and the many narrow absorption lines, allowing for a consistent method of determining the observed flux density. We started with 835 potential pairs of observations (with a time separation in the quasar rest frame between 2 days and 1 year) of which the 54 pairs of observation that had sufficient flux in both epochs to qualify for the sample are listed in Table 3.
A minimum signal to noise ratio needs to be found in order to eliminate the possibility of random noise creating numerous false positive variations that mask any structure in the temporal dependence of variability. A basic test is that for the very smallest intervals, on the time scales of days in the quasar rest frame, the variations should be distributed near zero with a “small” standard deviation. In the window, we find that an observed average flux density of is a good lower flux density cutoff. In the blue end of the SDSS spectra, the rms noise is more of a concern than sky subtraction errors for weak sources. Sky subtraction errors are evident from the fact that Ly absorption troughs can extend below zero. We estimate sky subtraction errors in the blue to be typically . Judging from the intensity below 0 in some of the deep Ly troughs, the effect can reach in a few cases. The worst case rms noise near is . Thus, an observed average flux density of has a signal to noise ratio, and more typically one has . The SNR over the range of (equivalent to a window of ) will be . Thus, random noise should exceed of the average flux density with a probability of less 0.01. Thus, random noise should not be an issue in our variability analysis, leaving systematic affects (such as large sky subtraction errors) as the main source of false positive variability.
There are two main issues that affect the veracity of SDSS for detecting modest variability. The first was establishing an acceptable SNR for the data. This was addressed above. The second is the calibration of quasar spectra that were observed as part of the Baryon Oscillation Spectroscopic Survey (BOSS) (Dawson et al., 2013). For BOSS quasar data, there is a calibration issue in the blue part of the spectrum. Specifically, in order to increase spectral sensitivity for quasars, the hole drilled in the plate over the fiber is centered on the displacement associated with 4000Å. However, the calibration of the star has a hole centered about the displacement associated with 5400Å. Due to atmospheric differential refraction (ADR), this inconsistency artificially raises the flux level in the blue for the quasar and underestimates the flux density of the standard star at 4000Å . The quasar affect is small, the observed wavelength is , so the ADR displacement of the light from the center of the fiber is minimal. However, the calibration standard star was observed with the fiber at a displacement centered at 5400Å and there is significant displacement from the center of the fiber at and these fluxes are significantly underestimated. These effects can enhance the detected quasar blue flux by more than 15%. Pictures of the plate and a discussion of this circumstance can be found in Dawson et al. (2013). We have corrected for the BOSS ADR affect in Table 3 with the following methodology. We retrieved the percentile (median) seeing values and the airmass from the header of the spectra. We then computed the displacement due to ADR between and for atmospheric air pressure at Apache Point (2800 m), and corrected for the relative enhancement of the flux with respect to flux assuming a 2D Gaussian shape for the seeing image of the quasar (seen within the 2” aperture of the BOSS plate holes), and integrating the light loss following the technique of Filippenko (1982) for a circular aperture.
The results of our search for EUV variability can be found in Table 3. The table is organized as follows. The first two columns are the source and the redshift. The next two columns are the date of the paired observations. Column (5) is the time interval in the rest frame of the quasar. Column (6) is the wavelength at the center of the 25Å window in the quasar rest frame. We try to make it as short of a wavelength as possible with the restriction that the Ly absorption does not suppress the average flux below our limit . We also try to avoid placing deep absorption troughs near the window edges as this can increase the potential for systematic errors in integration from epoch to epoch. We have not optimized this central wavelength. Columns (7) and (8) are the integrated fluxes in these windows evaluated in the quasar rest frame. Note that this can be converted to the average in the observer’s frame by dividing by . This number can fall below our limit of because these are BOSS corrected fluxes. The limit applies to uncorrected fluxes. This is appropriate since the BOSS observations were designed to increase signal to noise at the blue end of the spectra at the expense of an accurate absolute flux calibration. Column (9) is estimated from the CIV broad emission line in the SDSS spectrum (Park et al., 2013). Each individual measurement has about 0.3 dex uncertainty (Park et al., 2013). Thus, each estimate individually is not accurate enough for our purposes. However, recall the strategy that we declared in the previous section concerning the HST mass estimates. It is generally believed that by achieving a sufficient sample size (as in this section and the previous in combination) that the large scatter will be imprinted on a detectable backdrop of underlying physical trends associated with black hole mass. The next column is the variability computed with Equation (4). Column (11) is . The last column is the infall radius in the slim disk model associated with and for which the infall time (computed in Section 5) equals the time interval between observations in column (5).
Our criteria of an observed average flux density of , eliminates random noise as likely contributor to false positive identifications of variability at level. This leaves systematic affects as the most likely cause of potential false positive identification of variability. We explore the propensity of such behavior in the SDSS database by utilizing an empirical calibration method based on field galaxies (Yip et al., 2009). Using the assumption that the background field galaxies should not be variable, we derive a wavelength dependent correction. The correction varies across the plate, so it is only valid in a small region near each quasar. In general two plates cannot be compared unless they are identical plates at two different dates. Otherwise, we find that there are not enough overlapping galaxies with spectra on both plates in small neighborhoods of each quasar. We applied this method of calibration to our quasars in our EUV centered windows in order to get an estimate the accuracy of the flux integrations in out sample. We performed the following steps
A neighborhood of the plate was excised around each quasar, with a right ascension from the quasars and a declination arc-minutes from the quasar. The tight constraint on the declination follows from the fact that declination influences most of the light loss associated with the ADR since these sources were in general observed at small hour angle. We found that using the full plate instead of the contiguous region to the quasar can create errors as large as 15%.
The individual spectra of galaxies with were primary and non-primary observations of the same objects and were downloaded from the SDSS DR12 sky server.
By forming a wavelength dependent ratio of the flux densities of the field galaxies in the first epoch to the those in the second epoch in each quasar neighborhood, we determined a re-scaling wavelength dependent function that we applied to the measured quasar fluxes in the second epoch.
Unlike Yip et al. (2009), it was found that dividing the individual field galaxy spectra and then summing to retrieve the re-scaling function does not give stable results. Results are influenced by noise. Stable results are obtained by adding up the field galaxy spectra of each plate and then dividing the two plate cumulative spectra in order to obtain the re-scaling function.
We considered the plates with at least 20 field galaxies in common in the local quasar neighborhood defined in item (1) above. We then computed the re-scaling function and applied it to the EUV window for the quasar. The re-scaled flux was times the “uncorrected” flux, i.e. the fluxes that appear in Table 3. If we insist on at least 30 galaxies in common between the two plates the re-scaling factor in the EUV window is . The implication is that the SDSS data taken directly from the pipeline with our SNR minimum in the EUV window is typically accurate to within 10%. We verify the magnitude of the systematic uncertainty by considering the short time durations in Table 3. For short time durations, we expect virtually no variability in the EUV. For time intervals seconds (less than 3 weeks in the quasar rest frame), we find that the ratio of the EUV fluxes between epoch 1 and epoch 2 are which is consistent with most of the variation being from the systematic errors found from our field galaxy re-calibration of a few selected quasars.
A similar method has been used based on background stars instead of field galaxies (Wilhite et al., 2005). We tried to implement this method as well. However, we were able to find background stars in common in a local neighborhood of the quasar only in 3 cases. We preferred to stay with galaxies that allow the re-scaling of a larger number of plate pairs.
4 The Combined Data Sample
In this section, we combine the data in the HST and SDSS samples to improve the statistical analysis. This is viable because our selection criteria produced a common bias that selected very luminous quasars. This is a necessary consequence for an SDSS sample based on very high redshift quasars that are sufficiently luminous to shine brightly through the Lyman forest. Similarly this is a consequence of the practicality of allocating HST time for multiple observations of quasars in the EUV. It would be difficult to obtain telescope time for high resolution spectroscopy in the EUV for any source that was not very bright in the EUV. Such observations require one or more full orbits to accomplish even for the EUV bright quasars in our sample. Thus, we acquired quasars in both samples with very large and as evidenced by the calculations presented in Tables 2 and 3. In particular, and for the SDSS and HST samples, respectively. No statistical difference can be found in the parent population with a Kolmogorov-Smirnov test or a Wilcoxon rank sum test. The median are and for the SDSS and HST samples. However, the variance of the HST sample is much larger, as expressed by the mean values of and . No statistical difference is found between the two samples in a Wilcoxon rank sum test since the medians are similar, but the difference in variances make for a 0.012 probability that the samples are drawn from the same parent sample by random chance according to a Kolmogorov-Smirnov test. For our purposes, this is not very important since the distribution is statistically indistinguishable between the two samples, the larger variance of and similar medians just allows for a wider range of to be explored with a similar as a result of combining the samples. Thus, it is very appropriate from a statistical analysis point of view to combine these two samples in the following.
The left hand frame of Figure 2 combines the data in Tables 1 - 3. Approximately 30% of the objects are HST quasars in the redshift range, and 70% are SDSS quasars with . It is a scatter plot of the EUV variability computed with Equation (4) versus the time between individual epochs of observation. The errors are derived from the analysis and observation selection criteria of the last two sections. Recall that each observation has an uncertainty in the absolute flux measurement. This results in uncertainty ( is the relative variability that was defined in Equation (4))) in the variability calculation. The most striking aspect is the abrupt change in the distribution of the degree of variability when the time scale exceeds a threshold of . Significant variations are much more common above this threshold. More precisely stated, the variability is clearly larger for time intervals s compared to time intervals s. We try to explore the relevance of this visual appearance in terms of sampling density and measurement uncertainty in the remainder of this section. This is a rather startling find because there is a range of and (although rather narrowly distributed as discussed in the introductory paragraph of this section) evident in Tables 2 and 3. In order to “normalize” these times, we convert the times into geometrized units. The right hand frame shows the scatter plot, with the time converted to units of light travel time across the central black hole mass, in geometrized units. The transition to variability looks abrupt in terms of geometrized units as well. There appears to be an abrupt change in the distribution of the degree of variation that occurs at times in the right hand frame of Figure 2.
4.1 The Affects of the Minimum Flux Cutoff
In this section, we discuss the ramifications of altering the SDSS minimum flux cutoff. If chosen properly, this should not affect the results significantly, but should exclude some outliers with a variability enhanced by the random noise level. To start with, the left hand side of Figure 3 is a histogram of the distribution of variation as a function of time using the data in the left hand frame of Figure 2. Each bar of the histogram represents the average variation of the individual quasars in that bin of time (each individual quasar variation is computed per Equation (4)). There are two items that can be explored with this depiction of the data. Firstly, there is clearly much more variability for time separations s compared to time separations s. The variability is larger at level of statistical significance based on a Kolmogorov-Smirnov test. The other issue that is explored in this frame as well as the right hand frame of Figure 3, is the robustness of our methods and the threshold of detectable variability. Recall our efforts in the last section to quantify the accuracy of the SDSS flux calibrations using the technique of using adjacent field galaxies as an accurate empirical secondary calibration. Since most of the sample are SDSS quasars, the secondary flux calibration analysis indicates that bins in the histogram with a variability are consistent with no detected variability (just systematic flux measurement errors). The implication of this is that (possibly with the exception of a few mildly variable outliers) there is no measurable EUV variation for separation times s. More precisely since the range between s and s is sparsely sampled (see the left hand frame of Figure 2), a more rigorous statement supported by the data is, there is no measurable EUV variation (possibly with the exception of a few mildly variable outliers) for separation times s. Yet, significant variation is common for separation times .
The right hand frame of Figure 3 is used to explore the issue of the average flux density lower cutoff for the SDSS sample, . One issue with the cutoff is that it would automatically reject large variation for the fainter objects. Thus, highly variable sources are underestimated is general in all bins in the left hand frame of Figure 3. For the bins on the right hand side of the histogram with significant average variation, adding a couple more sources with a modest or large variation will not modify the average noticeably. However, adding one highly variable source to the bins on the left can modify the average variation proportionately much more than in the bins on the right. In order to explore this we looked at changing the average flux density lower cutoff for the SDSS sample to . This increased the number of epochs from 54 to 64 in the SDSS sample. The histogram for the lower cutoff is shown in the right hand frame of Figure 3. The general trend is the same as Figure 2. However, one highly variable source in the second bin raised this average considerably. Even so, there is clearly much more variability for time separations s compared to time separations s. The variability is still larger at level of statistical significance based on a Kolmogorov-Smirnov test.
4.2 Excess Variance Analysis
In order to accurately assess variability one must incorporate the uncertainty of the measurements, , into the measure of variability. In this section, we use the excess variance in order to segregate the role of the uncertainty (as determined above to be mainly systematic in nature due to our sample selection criteria) in our measurements of the intrinsic variability. As discussed above, in regards to the background field galaxy analysis, the systematic uncertainty in the absolute flux measurements of the SDSS sample in Table 3 is , or and . Our detailed discussion of the selection criteria for the HST data found a similar systematic uncertainty. Thus, on a statistical basis, one can only detect relative variability larger than . This concept is typically expressed as the excess variance, , (Nandra, George and Mushotzky et al., 1997)
where is the average of the flux measurements, . This is an odd application of excess variance considering that we only have two values of ”i”. However, it precisely describes the circumstance of interest: what constitutes variation larger than that expected by random chance? We note that we can express in terms of the relative variability, , in Tables 2 and 3 using Equation (4),
Figure 4 are scatter plots of versus the time between observations. On the left hand side, the time is measured in seconds and on the right hand side, the time is measured in estimated geometrized units of the black hole mass. The plots are similar to those of Figure 2. However, the uncertainty associated with each value of relative variability, , in Figure 2 is extricated by the defining relationship for in Equation (8). The minimum criteria for a variability detection consistent with the uncertainty in each individual flux measurement is . With this criteria for variability, 4/40 of the quasar observation pairs with a time separation sec are variable and 21/38 of of the quasar observation pairs with a time separation are variable. Figure 5 is a histogram highlighting the almost switch-on like behavior of the variability once the statistical uncertainty of each measurement is accounted for in the definition of variability. We note that averaged over a bin is similar to the structure function used in the UV photometry studies of quasars (Welsh et al., 2011; Macleod et al., 2012). However, there is not much evidence in the excess variance analysis of the EUV continuum variability of a gradual increase of variability for small time intervals to larger time intervals that was seen for the UV photometric structure functions. We eliminate the most variable source from the study in the right hand frame of Figure 5 in order to see if this is skewing the analysis. It does not. There is still more of a switch-on behavior than a gradual turn on of variability. We explore this in more detail in the remainder of the subsection.
|Bin||2.0 - 3.1 s||2.0 - 2.5 s||1.0 - 2.0 s|
|2.5 - 3.1 s||N/A||0.328||0.252|
|1.0 - 2.0 s||0.101||0.083||N/A|
|0.5 - 2.0 s||0.046||0.057||N/A|
|0.0 - 2.0 s||0.003||N/A|
|0.0 - 1.5 s||0.003||N/A|
|0.0 - 1.0 s||0.002||0.051|
|0.0 - 0.5 s||0.002|
We can characterize some of the behavior of denoted in Figures 4 and 5 by computing a probability matrix based on the Kolmogorov-Smirnov test. The results are given in Table 4. The matrix is a collection of Kolmogorov-Smirnov tests. Each entry compares two bins with at least 10 pairs of observations, The value that is recorded is the probability that the excess variance in the two bins are drawn from the same population. There are four statistical inferences that can be drawn from the probability matrix,
is larger when the time between observations in the quasar rest frame is s then it is when the time between observations in the quasar rest frame is s at statistical significance.
There is no statistical significance difference in whether the time between observations in the quasar rest frame is s or if the time between observations in the quasar rest frame is s.
There is no statistical significance difference in whether the time between observations in the quasar rest frame is s or if the time between observations in the quasar rest frame is s.
There is a marginal statistical difference in when the time between observations in the quasar rest frame is s as compared to the time between observations in the quasar rest frame is .
Consider these points in the context of Figure 5. Based on points 3 and 4 and Figure 5, there is an increased level of above s, but for time intervals s indicating that this increase is fairly mild and the variability is not statistically significant compared to the level of the measurement uncertainty. From point 1 and the Figure 5 for time intervals s, has increased to a level that far exceeds the threshold for statistically significant variability for many quasars. Yet, by point 2 and Figure 5, the propensity and magnitude of variability does not increase for time intervals s. This last point does not support a monotonic increase in as the time interval increases unless the rate of increase diminishes drastically for time intervals s. The overall behavior is well described by a broad threshold for measurable variability that occurs typically in the range of time intervals between observations of s above which the variability is roughly constant for time intervals less than 1 year. This behavior is similar to what has been seen in UV and optical quasar structure functions that are based on broad band photometry (Welsh et al., 2011; Macleod et al., 2012). The photometric structure functions show a steady increase in variability as the time between observations increases to a few hundred days. Beyond this time interval, the frequency of variability is roughly constant. The one difference is that the EUV continuum variability shows signs of increasing more abruptly after 100 days rather than an extrapolation of a steady increase from the small time intervals as in the photometric UV structure functions. This might be a result of the small number statistics in this study in the range s, or a fundamental difference associated with the fact that the EUV region is at the inner edge of the accretion disk. We note that a structure function cannot be used on the data presented here since we do not eliminate low variability sources as in photometric structure functions (Welsh et al., 2011). Consequently, the variability can be less than the estimated error (see the negative values of excess variance in Figures 4 and 5) rendering a putative structure function at time intervals, s, to values that are imaginary numbers (i.e., basically the square root of the negative excess variance).
In summary, the excess variance analysis supports a rapid change in the likelihood of significant variability for time intervals between observations that exceed some threshold above s. The threshold is loosely constrained due to small number statistics in the crucial range of s. The likelihood of variability levels off above above s after crossing this poorly resolved broad threshold. The threshold is a broad region because it clearly depends on many properties such as black hole mass and accretion rate. For this extreme luminosity quasar sample, virtually all the black hole masses are large and all the accretion rates very high making the threshold region somewhat concentrated in time and this enhances its appearance in Figure 5.
5 Slim Disk Models of EUV Variation
The previous section dealt with an empirical description of the data. This section is an attempt to make a speculative assessment on the implied size of the EUV emitting region based on models of slim disks (Sadowski et al., 2011). Being located just shortward of the peak of the SED, the EUV continuum that is sampled in Tables 1 - 3 likely represents the highest frequency optically thick thermal emission from the accretion flow. This emission therefore comes from the innermost regions of the accretion disk. We do not hypothesize on the nontrivial radiative transfer of this emission through an atmosphere. This includes, amongst other possibilities, electron scattering, Eardley et al. (1978); Czerny and Elvis (1987), and absorption that drives an outgoing wind (Laor and Davis, 2014). We just assume that a thermal EUV component was emitted from the disk, but acknowledge that its amplitude might be altered by processes in an enveloping atmosphere.
From Tables 2 and 3, we expect that the Eddington rate would be too large for the standard thin disk model to be accurate (Shakura and Sunyaev, 1973; Novikov and Thorne, 1973). For , several assumptions of thin disk theory and the approximations associated with the vertical integration of slim disks start to break down (Szuszkiewicz et al., 1996). Thus, we consider the most recent slim disk model solutions (Sadowski et al., 2011). Empirically, there is a change in several spectral properties at that provides observational evidence to an accretion mode change at . This has been characterized by the designation of Population A () and Population B quasars (Marziani et al., 2003b). Our sample is entirely Population A type spectra. The classic argument of Bardeen (1970) indicates that black holes should be spun up near their maximum allowed value if subject to prolonged accretion. Prolonged accretion seems reasonable for luminous quasars considering the He II proximity effect, high redshift quasars can accrete at a high rate for Myr (Syphers and Shull, 2014). The Bardeen (1970) conclusion has been corroborated empirically by considering the X-ray background, which requires high efficiency accretion associated with near maximal spin (Elvis et al., 2002). The assumption of the Bardeen argument is that the rotation of the galaxy determines a preferred sign of angular momentum for the accreted matter. We acknowledge that there might be cases in which a complicated accretion history might not lead to a large net spin up (i.e., possible flips in the sign of angular momentum associated with merger scenarios), but we do not consider this typical. Consequently, based on the Bardeen calculation and the He II proximity effect, we assume that most of our quasar sample have a high spin rate. We still consider low spin in the following, but high spin is the nominal configuration. The axisymmetric, time stationary spacetime (Kerr-Newman) metric is uniquely determined by three quantities, , and , the mass, angular momentum per unit mass, and the electric charge of the black hole respectively. For the sake of this paper, we consider negligibly small. In Boyer-Lindquist coordinates , the event horizon, , is
Based on the Bardeen calculation and the He II proximity effect, we assume that most of our quasars are near maximal rotation rate, . In particular, we will only discuss accretion disks with . As mentioned above, there can be individual quasars in which the accretion disk angular momentum vector changed directions over time, but these are considered outliers that will not affect the overall trends of the sample. Our fiducial value of spin is conservatively modest in our opinion, .
5.1 The Location and Inflow Velocity of the EUV Emitting Gas
The first thing to consider is the location of the optically thin surface layer that is the source of the EUV continuum emanating from the appropriate slim disk models for the Eddington rates and masses that are estimated in Tables 2 and 3. Figure 6 plots the expected location of the EUV emitting plasma for various slim disk models. The range of masses are typical of the luminous quasars that comprise the samples in Tables 2 and 3 and so are the ranges of Eddington luminosity. We choose a typical line of sight for a quasar that is from the accretion disk normal (Antonucci, 1993; Barthel, 1989). In the following, changing the line of sight to makes small changes that do not affect any of the conclusions. The plot is a cumulative distribution of the flux density at and that is emitted from the disk. For example, a value of 0.5 means that half of the observable flux density is emitted inside this radius. The flux density is computed in the cosmological rest frame of the host galaxy. Thus, it includes the Doppler shifting (including both the transverse Doppler shift and the gravitational Doppler shift) of the disk emission. The disks are characterized by a vertically averaged viscosity parameter of . The emitted flux distribution from the slim disk is very weakly dependent on the value of (A. Sadowski private communication 2016). For all the models in Figure 6, 50% of the emission comes from inside a radius that is between 5M and 6M. For higher black holes spins these curves shift to the left and the 50% point shifts that way as well. As an extreme example, if and , the 50% point is at 3.8 M for the flux density. If an extreme change happens to the emissivity of the optically thin surface layer that is the source of 50% - 60% of the the flux density, this could clearly create the 25% - 30% variability that typifies the behavior of most of the variable quasars in Tables 2 and 3 and Figures 2 - 5. Based on Figure 6, we look for a variability time scale associated with the optically thin surface layer at in the slim diks models that is on the order of s or 1500 - 2500 M in geometrized units based on Figures 2 - 5.
The primary parameter to consider is the infall velocity of the optically thin surface layer. The infall time from a particular radius in the EUV portion of the disk can be computed from the radial velocity in Boyer-Lindquist coordinates, . Figure 7 plots the radial velocity for various relevant sets of parameters. The figure allows for an exploration of the dependence of on the viscosity parameter, the Eddington ratio and the black hole spin. One thing to notice is the local maximum in at small radius and high spin. This is a consequence of gravitational redshifting and the centrifugal barrier. The Boyer-Lindquist coordinates are those of observers that are stationary with respect to asymptotic infinity. These observers never see any trajectory cross the event horizon. The trajectories approach the event horizon asymptotically slowly, the freezing of the flow (Punsly, 2008),
When the spin increases, the inner edge of the disk approaches the event horizon. Thus, for high spin, the gravitational redshift effects will be large near the inner edge of the disk and in the plunge region (the region where gravitational attraction towards the black hole overwhelms centrifugal force). Since the infalling gas never actually reaches the event horizon in Boyer-Lindquist coordinates, the time for infall is not well defined. Furthermore, the solutions of Sadowski et al. (2011) do not extend all the way to the event horizon. These issues are circumvented by the condition that infall is measured to the r coordinate that is halfway between and the inner edge of the disk. We note two points:
In general, the time for the gas to traverse this gap from the inner edge of the disk to this midpoint is at most of few percent of the infall times that we compute.
Furthermore, due to gravitational redshifting the flux detected by observers at asymptotic infinity is negligible for emission between this midpoint and the event horizon.
For these two reasons, this simple definition of the infall end point leads to very robust results that are not sensitive to precise location of the end point. The primary result of Figure 7 is that the magnitude of within the accretion disk is orders of magnitude less than the speed of light. Thus the infall times in geometrized units are orders of magnitude larger than the radius from which the infall begins! We also note the scalings displayed in Figure 7, indicate that increases (decreases) with increases (decreases) of and and decreases (increases) of .
Notice that most of the curves in Figure 7 are based on an extremely large value of the viscosity parameter, . The motivation for this is the following. The slim disk models are very optically thick except for a surface layer of optically thin gas. This is the gas that is responsible for the observed EUV flux. The slim disk has vertical structure and the velocity that is calculated in Sadowski et al. (2011) is vertically averaged. We are actually only interested in the radial velocity of the optically thin surface layer. This is beyond slim disk models and requires full 3-D calculations. Although the vertical structure of slim disks is not well known, recent 3-D numerical simulations indicate that near the inner edge of the disk, surface layers move inward approximately twice as fast as the vertical average (Sadowski, 2016). This can be represented by a higher effective in the surface layers that is twice as large as the vertical average (A. Sadowski, private communication 2016). Hence, the focus on in the Figure 7.
The first thing to consider is the implication of Figure 7 for the infall times for the surface layers of the slim disk that are responsible for the EUV continuum in Figure 6. The cumulative flux density in Figure 6 and the infall times in bottom frame of Figure 8 are both functions of the radial coordinate, . We exploit this in the top frame of Figure 8 by plotting these implicit functions of , without any direct reference to . The top frame of Figure 8 is a plot of the fraction of the cumulative flux density inside of a radius, R, versus the infall time from a radius, R, based on the radial velocities that are calculated as in Figure 7. For typical values of black hole mass and Eddington ratio s, is 50% to 65% of the total emitted from the disk. Large changes to the surface layer in this region will affect a sufficient fraction of the EUV emitting gas in order to account for the variability in Tables 2 and 3 as well as Figures 2 - 5. Thus, we have established that the infall time associated with the EUV emitting surface layer of the appropriate slim disk model is consistent with the variability time scale indicated in Figure 5. The bottom frame of Figure 8 uses the values from Figure 7 to integrate trajectories and determine the infall times for different starting points. The infall time in geometrized units is plotted versus the radial coordinate of the initiation point of the infall for various designated slim disk models. The plot indicates that the infall time decreases (increases) with increases (decreases) of and and decreases (increases) of . The horizontal dotted lines represent the range of infall times in geometrized units for the significantly variable sources in our combined sample of quasars in the right hand frame of Figure 1.
Figure 9 is more ambitious. It is a plot of the infall radius, , for the quasars in Tables 2 and 3 based on the estimated mass, and spin of the black hole, versus the variability. These values were computed from the slim disk models and were tabulated in column (3) of Table 2 and column (12) of Table 3. On the top of the figure is a plot of the relative variability versus the maximum infall radius, . The maximum infall radius was obtained by associating the infall time from (as depicted in bottom frame of Figure 8) with the interval between observations. The calculation assumes the fiducial model of . The bottom plot in Figure 9 is a plot of the excess variance as a function of . The advantage of these plots is that they remove the dependence of black hole mass and Eddington rate. The disadvantage is that the data is randomized by the significant systematic uncertainty in the black hole mass estimates. However, as an ensemble of data, one expects that the general trends should be valid. The scatter plots indicate that in general, significant variability is seen for infall radii, . This location is consistent with the distribution of EUV flux depicted in Figure 6.
There is a curious extreme point in the scatter plots in Figure 9. Consider the value of 5.5 M for the maximum infall radius of PG 1115+080. This has a natural explanation within the slim disk models. Recall from Table 3 the extremely large black hole mass estimate for this object . This is a cooler disk than our nominal disk with . This shifts the 50% cumulative distribution point down to 5M in Figure 6. We also see from Figure 6, that the 50% point is located at smaller values of radius than the the 50% cumulative distribution point. This is probably more relevant considering the spectrum of PG 1115+080 in Figure 1.
5.2 A Crude Model of EUV Variability
The discussions associated with Figures 6 - 9 indicate that there is a possible link between plasma infall times from the EUV emitting surface layer of the slim disk models and the EUV continuum variability detected in Tables 2 and 3. With this motivation, we consider the most straightforward explanation of EUV variability in the context of slim disk models of quasars with the aid of the schematic diagram in Figure 10. We acknowledge that this is a very speculative and crude model. In the model, a significant change in EUV luminosity is associated with an evacuation of the inner disk plasma and a replenishment with new plasma with a different emissivity. More specifically, since the equatorial regions of the disk are optically thick, a significant change in EUV luminosity is associated with an evacuation of the optically thin surface layer of inner disk plasma and a replenishment with a new surface layer of plasma with a different emissivity.
Figure 10 shows the basic configuration. The black hole is depicted as a black disk. The optically thin surface layer of the accretion disk is an annulus with an emissivity that is depicted in grey scale, where black is the highest emissivity and white is zero emissivity. In the schematic, there are two concentric circles that demarcate the boundary of the EUV region. The boundaries can move slightly as the accretion state changes, but this is a higher order affect that is beyond this schematic representation. We describe the particular circumstance depicted in each frame from left to right. The frame on the left indicates an enhanced emissivity in dark shading that envelopes the EUV region and extends beyond it. This would be the configuration of a significant elevation of EUV flux in magnitude and time. Within standard accretion disk theory, this increase is associated with an increase of the effective temperature of the surface layer (Sadowski et al., 2011). The second frame shows an epoch with decreased emissivity. If we were to compare the two epochs, it would correspond to some of the larger variations in our combined sample of the last section. The last two frames indicate the different types of EUV flares. A large variation in Table 3 , would require that virtually all of the EUV region change emissivity. For example, the third frame represents such a high state that we consider an variation. This is in distinction to the last frame that has only a small region of the EUV annulus filled with high emissivity gas and this might be represent a short period of slightly elevated flux with a variation parameter .
The idea of associating the variability of the EUV continuum in luminous quasars to the emptying and refilling of the inner annulus formed by the optically thin surface layer of an accretion disk naturally produces dispersion in the amplitude of variation and the time frame required for significant variation. In principle, local physics or advection from larger radii can create different magnitudes of the changes in the surface layer emissivity (surface temperature) in the inner disk. Different time scales for variability naturally arise from variations in black hole mass, accretion rate and black hole spin. Based on Figures 6 - 8, the crude model should be capable of explaining the coarse features of EUV continuum variability in highly luminous quasars (Figures 2-5 and 9), the observed spread in variability magnitude and and the observed spread in the time scale for significant variation.
This article explores the variability of the EUV continuum emission in a sample of highly luminous quasars. The high bolometric luminosity of these quasars was necessary in order to insure high signal to noise absolute flux measurements of the EUV continuum in HST and SDSS observations. The end result is not a sample covering all of quasar parameter space, but only massive black holes and high Eddington rates (most of the sample has and ). The analysis utilized pairs of spectral observations from the SDSS and HST archives. Since the variability that was detected is typically less than O(1), careful considerations needed to be made in order to screen out false variation that was a mere consequence of inadequate calibration technique. The majority of this effort was consumed by the technical details of the flux calibrations in order to minimize this deleterious effect. In the end, we had a modest sample of 78 pairs of quasar observations (54 SDSS pairs and 24 HST pairs). In spite of the modest sample size, we are able to find a characteristic time scale for measurable variability. At the crudest level, the variability is larger when the sampling time between observations in the quasar rest frame is sec compared to sec (at level of statistical significance according to a Kolmogorov-Smirnov test). Based on an excess variance analysis, for time intervals sec in the quasar rest frame, of the quasars show evidence of EUV variability. Similarly, for time intervals sec in the quasar rest frame, of the quasars show evidence of EUV variability (see Figure 5 and the related discussion). The amount of variability levels off between sec and sec (1 yr). This temporal behavior is indicative of a threshold time scale for variability as opposed to a steady monotonic increase. In Section 4, our findings were compared and contrasted with similar behavior that is seen in photometry based quasar structure functions. Marginal statistical evidence was presented that the threshold time scale for significant variability is spread out possibly as wide as sec to sec.
This broad threshold was considered to be largely a consequence of the spread of black hole mass and Eddington rates. We tried to remove the black hole mass dependence by converting the time between observations to geometrized units of black hole mass, M. In these units the broad threshold appears to occur between 1800M and 2000M (see the right hand frames of Figures 2 and 4). However, this did not remove the dependence on Eddington rate. This could only be achieved by appealing to theoretical models of accretion.
We pursued the possibility that this threshold time scale was directly associated with a characteristic spatial dimension of the EUV emitting region in Section 5. This analysis cannot be done purely empirically, but requires model dependent assumptions. We chose the slim disk optically thick accretion model as our fundamental assumption. In order to characterize the black hole, we noted classical theoretical arguments and observational analysis that for a large sample of highly luminous quasars most of the central supermassive black holes should be spun up to near their maximal allowed values (Bardeen, 1970; Elvis et al., 2002). For these reasons, even though we explored a broad range of black hole spins, our fiducial value of spin was chosen to be . In Figure 6, we showed that of the EUV emission is emitted from radii, and of the EUV emission is emitted from for with parameters that are characteristic of our sample, and . Thus, the slim disk models, if they are correct, imply that the EUV continuum variability is occurring as a consequence of local physics in this region. Next, we tried to find a plausible time scale within the slim disk models that is associated with the dimensions of this region.
The most direct time scale to look for was the accreting gas infall time to the plunge region. In order to represent the inflow velocity of the optically thin surface layer responsible for the EUV emission, we chose an extremely large viscosity parameter, in the optically thin surface layer that is responsible for the EUV continuum (i.e roughly twice the value in the bulk of the disk) as indicated in 3-D modeling of slim disks. The basic logic of Section 5 can be shown in skeleton form. Our time interval is s or 1800M. The sample is dominated by powerful quasars that have central black hole masses of (see Tables 2 and 3) which is cm in geometrized units. This equates to a light travel time of s. This allows us to write the time interval in geometrized units as M. For radial inflow velocities in our fiducial model with and (based on Tables 2 and 3), the infall velocity is = 0.002c - 0.005c (see the slim disk models in Figure 7). The infall distance in a time is cm. Dividing this by to get into geometrized units indicates that . We equated this distance to the radius from which gas must fall inward to the black hole in order to produce the observed variability time scale. For our fiducial value of the spin , 75% EUV emitting gas in the slim disk models is located at . Thus, this radial infall is a viable dynamic for variability.
We tried two more analyses in order to see if the data was consistent with the infall idea. The first was simply to plot the amount of EUV flux density that emerges inside of a radius R versus the infall time from that radius. This was plotted in the top frame of Figure 8. We found that for black hole parameters that were relevant to our sample, s was associated with infall from a radius of 5M - 6M. Furthermore, 50% - 65% of the total EUV flux emerged from inside of this region. Thus, the infall time scale was consistent with the slim disk model and the observed variability. The most ambitious test was an attempt to remove the affects of the Eddington ratio in Figure 9. We computed the infall time as a function of radius in geometrized units based on slim disk models associated with the estimated black hole mass and Eddington rate for each quasar. We equated the infall time with the time between observations. The radius from which the infall time equals the time lag between observations was designated as the infall radius. In Figure 9, we found that variability started to increase significantly beyond infall radii . This was the third analysis (Figures 2-6, Figure 8 and Figure 9) that found similar dimensions for the source of the EUV variability, indicating consistency between the infall time from slim disk models and the EUV continuum variability time scale.
Finally, we speculated that the agreement between the slim disk infall time scale, the variability time scale and the location of the EUV emission in slim disk models was not coincidental. We proposed a very crude model in which the accretion of the optically thin surface layer within M of the inner edge of the accretion disk and replenishment with new gas of a different emissivity (different surface temperature) was a viable candidate physical mechanism for the EUV continuum variability in highly luminous quasars.
- Antonucci (1993) Antonucci, R.J. 1993, Annu. Rev. Astron. Astrophys., 31, 473
- Assef et al. (2011) Assef, R. 2011 ApJ 742 93
- Bardeen (1970) Bardeen, J. 1970 Nature 226 64
- Barthel (1989) Barthel, P. 1989 ApJ 336 606
- Czerny and Elvis (1987) Czerny, B. and Elvis, M. 1987, ApJ 321 305
- Dawson et al. (2013) Dawson, K. 2013 AJ 145 10
- Davis and Laor (2011) Davis, S., Laor, A. 2011, ApJ 728 98
- Dietrich et al. (2009) Dietrich, M. Smita, M. Grupe, D. and Komossa, S. 2009, ApJ 696 1908
- Eardley et al. (1978) Eardley, D., Lightman, A., Payne, D., Shapiro, S. 1978, ApJ 224 53
- Elvis et al. (2002) Elvis, M., Risaliti, G., Zamorani, G. 2002, ApJL 565 75
- Evans and Koratkar (2004) Evans, I. and Koratkar 2004, ApJS 150 73
- Filippenko (1982) Filippenko, A. 1982, PASP 94 715
- Finn et al. (2014) Finn, C. et al. 2014, MNRAS 440 3317
- Ganguly et al. (2007) Ganguly,, R. et al. 2007, ApJ 665 990
- Heap et al. (1995) Heap, S. et al. 1995, PASP 107 871
- Kittel (1961) Kittel, C. 1961, Elementary Statistical Physics. (Wiley, New York)
- Laor and Davis (2014) Laor, A., Davis, S. 2014 ApJ 428 3024
- Liu et al. (2014) Liu, T,. Wang, J-X, Yang, H., Zhu, F-F, Zhou, Y-Z, 2015 ApJL 783 106
- Lynden-Bell and Rees (1971) Lynden-Bell, D., Rees, M. 1971, Mon. Not. R. Astr. Soc. Lett. 152 461
- Macleod et al. (2012) Macleod, C. 2012, ApJ 753, 106
- Marshall et al. (1997) Marshall, H. 1997, ApJ 479, 222
- Marziani et al. (2003a) Marziani, P., Sulentic, J. W., Zamanov, R., Calvani, M., Dultzin-Hacyan, D., Bachev, R., Zwitter, T. 2003 ApJS 145 199
- Marziani et al. (2003b) Marziani, P., et al. 2003 MNRAS 345 1133
- Nandra, George and Mushotzky et al. (1997) Nandra, K., George, I. M., Mushotzky, R. F., Turner, T. J. and Yaqoob, T. 1997 ApJ 476 30
- Neira and Constantinides (1996) Neira. L. and Constantinides, A. 1996 Signal Processing 50 223
- Novikov and Thorne (1973) Novikov, I. and Thorne, K. 1973, in Black Holes: Les Astres Occlus, eds. C. de Witt and B. de Witt (Gordon and Breach, New York), 344
- Park et al. (2013) Park, D., Woo., J.-H., Denney, K,. Shin, J. 2013 ApJ 770 87
- Penna et al. (2010) Penna, R. et al 2010 MNRAS 408 752
- Punsly (2008) Punsly, B. 2008, Black Hole Gravitohydromagnetics, second edition (Springer-Verlag, New York)
- Punsly (2014) Punsly, B. 2014 ApJL 797 33
- Punsly (2015a) Punsly, B. 2015 The Formation and Disruption of Black Hole Jets, ASSL, 414 ISBN 978-3-319-10355-6. Springer International Publishing Switzerland, p. 149
- Punsly (2015b) Punsly, B. 2015 ApJ 806 47
- Punsly and Rodriguez (2013) Punsly, B., Rodriguez J. 2013, ApJ 764 173
- Punsly et al. (2009) Punsly, B., Igumenshchev, I. V., Hirose, S. 2009 ApJ 704 1065
- Punsly and Marziani (2015) Punsly, B. and Marziani 2015 MNRAS Lett 453 16.
- Punsly et al. (2016) Punsly, B.; Marziani, P.; Kharb, P/; O’Dea, C.; Vestergaard, M.. 2016 ApJ 812 79
- Sadowski et al. (2011) Sadowski, A. et al. 2011, A & A 527 A17 1
- Sadowski (2016) Sadowski, A. 2016, submitted to MNRAS http://arxiv.org/abs/1601.06785
- Shakura and Sunyaev (1973) Shakura, N., Sunyaev, R. 1973, A & A 24 337 1
- Shen and Liu (2012) Shen, Y., & Liu, X. 2012, ApJ 753 125
- Shields et al. (2003) Sheilds, G. et al. 2003, ApJ 583 124
- Soderblom and Sherbert (1997) Soderblom, D. and Sherbert, L. 1997 HST Calibration Workshop Space Telescope Science Institute, 1997 S. Casertano, et al., eds.
- Steidel and Sargent (1991) Steidel, C. and Sargent, W. 1991 ApJ 383 433
- Stevans et al. (2014) Stevans, M., Shull, M., Danforth, C., Tilton, E. 2014ApJ 794 75
- Sun and Malkan (1989) Sun, W.-H., and Malkan, M. A 1989, ApJ 346 68
- Syphers and Shull (2014) Syphers, D. and Shull, J.M. 2014, ApJ 784 428
- Szuszkiewicz et al. (1996) Szuszkiewicz, E., Malkan, A., and Abramowicz, M. A. 1996, ApJ 458 474
- Telfer et al. (2002) Telfer, R., Zheng, W., Kriss, G., Davidsen, A. 2002 ApJ 565 773
- Trakhtenbrot & Netzer (2012) Trakhtenbrot, B. and Netzer, H. 2012, MNRAS 427 3081
- Tsvetkova et al. (2010) Tsvetkova, V. et al. 2010, MNRAS 406 2764
- Welsh et al. (2011) Welsh, B., Wheatley, J., and Neil, J. 2011, A & A 527 A15
- Wilhite et al. (2005) Wilhite, B. et al. 2005, ApJ 633 638
- Yip et al. (2009) Yip, C. et al. 2009, AJ 137 5120
- Zheng et al. (1997) Zheng, W. et al. 1997 ApJ 475 469
- Zhu et al. (2012) Zhu, Y. et al. 2012, MNRAS 424 2504