© ESO, 2011

1. Introduction

Precise estimation of the time delay between components of lensed active galactic nuclei (AGN) is crucial for modeling the lensing objects. In turn, more accurate lens models give better constraints on the Hubble constant. More than 200 strong lens systems have been discovered, most of them in recent years by dedicated surveys such as the Cosmic Lens All-Sky Survey (Myers et al. 2003; Browne et al. 2003) and the Sloan Lens ACS Survey (Bolton et al. 2004). The launch of the Fermi satellite (Atwood et al. 2009) in 2008 offered the opportunity to investigate gravitational lensing phenomena with high-energy gamma rays. The observation strategy of Fermi-LAT, which surveys the whole sky in 190 min, allows a regular sampling of quasar light curves with a period of a few hours. This paper deals with the observation and estimation of time delays in strong lensing systems and not with the detection and use of microlensing phenomena such as in Torres et al. (2003).

The multiple images of a gravitational lensed AGN cannot be directly observed with high energy gamma-ray instruments such as Fermi-LAT, Swift, or ground-based Cerenkov telescopes, owing to their limited angular resolutions. The angular resolution of these instruments is at best a few arcminutes (in the case of HESS), when the typical separation of the images for quasar lensed by galaxies is a few arcseconds. This paper is concerned with estimating the time delay in strong lenses when only spatially unresolved data are available. Spectral analysis techniques and echo detection methods are investigated. A similar approach to detecting GRB lensing has been proposed by Wambsganss (1993).

Our methods of time-delay estimation have been tested on simulated light curves and on the Fermi LAT observations of the very bright radio quasar PKS 1830-211, for which the time delay was previously estimated by Lovell et al. (1998) using radio observations. The paper is organized as follows. We first give a very brief summary of properties of PKS 1830-211 and of Fermi LAT data towards this AGN. We next introduce the methods for estimating the time delay. The last section is devoted to measuring the time delay between the two compact components of PKS 1830-211.

2. The PKS 1830-211 gravitational lens system

The AGN PKS 1830-211 is a variable, bright radio source and an X-ray blazar. Its redshift has been measured to be z = 2.507 (Lidman et al. 1999). The blazar was detected in the γ-ray wavelengths with EGRET. The EGRET source was associated with the radio source by Mattox et al. (1997). The classification of PKS 1830-211 as a gravitationally lensed quasi-stellar object was first proposed by Pramesh Rao & Subrahmanyan (1988). The lensing galaxy is a face-on spiral galaxy, which was identified by Winn et al. (2002) and Courbin et al. (2002), and it is located at redshift z = 0.89 (Wiklind & Combes 1996).

PKS 1830-211 is observed in radio as an elliptical ring-like structure connecting two bright sources roughly one arcsecond distant (Jauncey et al. 1991). The compact components were separately observed by the Australia Telescope Compact Array at 8.6 GHz for 18 months. These observations and the subsequent analysis by Lovell et al. (1998) give a magnification ratio between the 2 images of 1.52 ± 0.05 and a time delay of $\mathrm{2}{\mathrm{6}}_{-5}^{\mathrm{+}\mathrm{4}}$ days. A separate measurement of a time delay of $\mathrm{2}{\mathrm{4}}_{-4}^{\mathrm{+}\mathrm{5}}$ days has been made by Wiklind & Combes (2001) using molecular absorption lines.

3. Fermi LAT data on PKS 1830-211

PKS 1830-211 has been detected by the Fermi LAT instrument with a detection significance above 41 Fermi Test Statistic (TS), equivalent to a 6σ effect (Abdo et al. 2010). The long-term light curve is presented in Fig. 1 with a two-day binning. The data analysis from this paper uses a two-day binning, which provides a sufficient photon statistic per bin with a time span per bin that is much shorter than 28 days. The data analysis was cross-checked by binning the light curve into one day and 23-h bins, with similar results. The data were taken between August 4, 2008 and October 13, 2010 and processed by the publicly available Fermi Science Tools version 9. The v9r15p2 software version and the P6_V3_DIFFUSE instrument response functions were used. The light curve was produced by aperture photometry selecting photons from a region with radius 0.5 deg around the nominal position of PKS 1830-211 and energies between 300 MeV and 300 GeV.

4. Data processing and method

4.1. Idea

If a distant source (in our case an AGN) is gravitationally lensed, the light reaches the observer by at least two different paths. We assume here that there are only two light paths. In reality, the light curves of the two images are not totally identical since (in addition to differences due to photon noise) the source can be microlensed in one of the two paths. Neglecting for the moment the background light and the differences due to microlensing, the observed flux can be decomposed into two components. One of the components is the intrinsic AGN light curve, given by f(t), with Fourier transform $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{f}\end{array}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathit{.}$ The other component has a similar time evolution to the first one, but is shifted in time with a delay a. In addition, the brightness of the second component differs by a factor b from that of the first component, so that it can be written as bf(t + a), and its transform to the Fourier space gives $\mathit{b}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{f}\end{array}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{\pi }\mathrm{i}\mathit{\nu a}}$.

The sum of two component gives g(t) = f(t) + bf(t + a), which transforms into $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{g}\end{array}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{f}\end{array}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{b}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{\pi }\mathrm{i}\mathit{\nu a}}\mathrm{)}$ in Fourier space. The power spectrum P_ν of the source is obtained by computing the square modulus of $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{g}\end{array}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}$: ${\mathit{P}}_{\mathit{\nu }}\mathrm{=}\mathrm{\left|}\begin{array}{}\mathrm{˜}\\ \mathit{g}\end{array}\mathrm{\right(}\mathit{\nu }\mathrm{)}{\mathrm{|}}^{\mathrm{2}}\mathrm{=}\mathrm{|}\begin{array}{}\mathrm{˜}\\ \mathit{f}\end{array}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}{\mathrm{|}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{b}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{b}\cos \mathrm{(}\mathrm{2}\mathit{\pi \nu a}{\mathrm{\left)}}^{\mathrm{\right)}}\mathit{.}$(1)The measured P_ν is the product of the “true” power spectrum of the source times a periodic component with a period (in the frequency domain) equal to the inverse of the relative time delay a. The microlensing of one of the components, if taken into account, gives a modulation of the amplitude of the oscillatory pattern at low frequencies. Since the typical duration of a caustic crossing microlensing event is a few months, only frequencies under 3 × 10^-7 Hz will be affected.

The usual way to measure the time delay a is to calculate the autocorrelation function of f(t). This method was investigated by Geiger & Schneider (1996). Computation of the autocorrelation of a light curve with uneven sampled data is described in Edelson & Krolik (1988). The Fermi light curve of PKS 1830-211 has very few gaps, and only one notable four-day gap. The missing data have been linearly interpolated. However, simulations with an artificial gap have shown that the results of the present paper are not affected by this gap. The autocorrelation function can be written as the sum of three terms. One of these terms models the “intrinsic” autocorrelation of the AGN, decreasing with a time constant λ. If λ is larger than a, the autocorrelation method fails, because the time delay peak merges with the intrinsic component of the AGN. Another potential problem with the autocorrelation method is the sensitivity to spurious periodicities such as the one coming from the motion and rotation of the Fermi satellite.

The periodic modulation of P_ν suggests the use of another method, based on the computation of the power spectrum of P_ν, noted D_a. This method is similar in spirit to the cepstrum analysis (Bogert et al. 1963) used in seismology and speech processing. If $\mathrm{\left|}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{f}\end{array}\mathrm{\right(}\mathit{\nu }\mathrm{)}{\mathrm{|}}^{\mathrm{2}}$ were a constant function of ν, D_a would have a peak at the time delay a. In the general case, D_a is obtained by the convolution of a Dirac function, coming from the cosine modulation, by the Fourier transform of the function: $\begin{array}{}\mathrm{˜}\\ \mathit{h}\end{array}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\mathrm{\left|}\begin{array}{}\mathrm{˜}\\ \mathit{f}\end{array}\mathrm{\right(}\mathit{\nu }\mathrm{\left)}{\mathrm{|}}^{\mathrm{2}}\mathit{W}\mathrm{\right(}\mathrm{0}\mathit{,}{\mathit{\nu }}_{\max }\mathrm{)}$(2)where W(0, ν_max) is the window in frequency of P_ν and ν_max is maximum available frequency. The Fourier transform h(a) of $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{h}\end{array}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}$ defines the width of the time delay peak in the double power spectrum D_a. For instance if $\mathrm{\left|}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{f}\end{array}\mathrm{\right(}\mathit{\nu }\mathrm{)}{\mathrm{|}}^{\mathrm{2}}\mathrm{=}{\mathrm{e}}^{\mathrm{-}\mathit{\lambda \nu }}$ and λW ≫ 1, then the time delay peak in D_a has a Lorentzian shape with a FWHM of λ/π. For a typical value of λ = 10 days, one has a FWHM of 3 days.

In the next section we describe the calculation of P_ν and D_a and illustrate the procedure with Monte Carlo simulations.

4.2. Power spectrum

To avoid problems arising from the finite length of measurements, sampling and aliasing, we use the procedure for data reduction described by Brault & White (1971) and Press et al. (2007).

We first divide the whole light curve into several segments of equal length. These segments are overlapped to obtain more segments with a sufficient number of points. After the interpolation of missing data and standard operations (oversampling, zero-padding), the power spectrum is calculated in every segment. Finally, the power spectrum is averaged over all segments. The power spectrum obtained with an artificial light curve is shown in Fig. 2.

The artificial light curve was produced by summing three simulated components. The light curve of PKS 1830-211 shown in Fig. 1 does not exhibit any easily recognized features, but it does have a rather random-like aspect. The first component was thus simulated as white noise with a Poisson distribution. It would be more realistic to use red noise instead of white noise, but the latter is sufficient for most of our purposes, such as computing D_a. The second component is obtained from the first by shifting the light curve with a 28-day time lag. The effect of differential magnification of the images has also been included. The background photon noise was taken into account by adding a third component with a Poisson distribution.

The mean number of counts per two-day bin for PKS 1830-211 is 5.42. This value was used to generate artificial light curves. The first and second components account for 80% of the simulated count rate and the rest is contributed by the Poisson noise.

4.3. Time delay determination

The methods of time delay determination use the power spectrum P_ν as described in the previous section. The simulated P_ν presented in Fig. 2 shows a very clear periodic pattern. From Eq. (1) we know that the period of the observed oscillations equals the inverse of the time delay between the images.

Our preferred approach was to calculate the double power spectrum D_a. As in Sect. 4.2, the power spectrum P_ν has to be prepared before undergoing a Fourier transform to the “time delay” domain. The low-frequency part (ν < 1/55   day^-1) of P_ν is cut off. This cut arises because of the strong power observed at low frequencies in the power spectrum of PKS 1830-211. The high-frequency part of the spectrum P_ν is also removed because the power at high frequency is low. (It goes to 0 at the Nyquist frequency.) Calculation of D_a proceeds like in Sect. 4.2, except that the P_ν data are bent to zero by multiplication with a cosine bell. The D_a distribution is estimated from five segments of the light curve. In every bin of the D_a distribution, the estimated double power spectrum is given by the average over the five segments. The errors bars on D_a are estimated from the dispersion of bin values divided by two (since there are five segments). Due to the random nature of the sampling process, some of the error bars obtained are much smaller than the typical dispersion in the D_a points. To take this into account, a small systematic error bar was added quadratically to all points. The result (with statistical error bars only) is presented in Fig. 3.

As described in Sect. 4.2, we simulated light curves with a time delay of 28 days. A peak is apparent near a time delay of 28 days on the D_a distribution shown in Fig. 3. The points just outside the peak are compatible with a flat distribution. Including also the points in the peak gives a distribution that is incompatible with a flat distribution at the 12.9 sigma level. The parameters of the peak were determined by fitting the sum of a linear function for the background plus a Gaussian function for the signal. In the case shown in Fig. 3, the time delay estimated from D_a is 28.35 ± 0.56 days.

As mentioned in Sect. 4.1, the usual approach for time delay estimation is to compute the autocorrelation of the light curve. The auto-covariance is obtained by taking the real part of the inverse Fourier transform of P_ν. The auto-covariance is normalized (divided by the value at zero time lag) to get the autocorrelation. The autocorrelation function of an artificial light curve simulated as in Sect. 4.2 is presented in Fig. 4.

A peak with a significance of roughly 16σ is present at 27.85 ± 0.14 days, however, the significance of this peak is overestimated since light curves are simulated with white noise instead of red noise. The autocorrelation function of a light curve driven by red noise is given by e ^− a/λ. For our simulated light curves, λ = 0, so that the peak is a hardly affected by the background of the AGN.

For the simulated light curves, both approaches of time delay determination give reasonable and compatible results.

5. Results and discussion

The results for real data were obtained with the same procedure as for the simulated light curves. Figure 5 shows the power spectrum P_ν calculated from the light curve of PKS 1830-211. An oscillatory pattern is clearly visible in the power spectrum. It is similar to the pattern expected from the simulations shown in Fig. 2. The autocorrelation function and the D_a distribution calculated for real data are shown in Figs. 7 and 6. A peak around 27 days is seen in both distributions. Several other peaks are present on the autocorrelation function as already noted by Geiger & Schneider (1996, see their Fig. 1). The peak around five days in the D_a distribution is likely to be an artefact of the time variation of the exposure of the LAT instrument on PKS 1830-211. Using the method described in Sect. 4.3, the significance of the peak around 27 days is found to be 1.1σ in the autocorrelation function and 4.2σ in the double power spectrum D_a. Fitting the position of the peak gives the time delay of a = 27.1 ± 0.6 days for the D_a distribution. The fit of the autocorrelation function to a Gaussian peak over an exponential background gives a = 27.1 ± 0.45 days. In both cases, the quoted error is derived from the fit.

The double power spectrum distribution obtained for PKS 1830-211 provides the first evidence of gravitational lensing phenomena in high-energy gamma rays. The evidence is still at the 4.2σ level but will likely improve by a factor of 2 over the lifetime of the Fermi satellite. Thanks to the uniform light curve sampling provided by the Fermi LAT instrument, it is not necessary to identify features on the light curve to apply Fourier transform methods. The example of PKS 1830-211 shows that the method works in spite of the low photon statistics. Possible extensions of the present work are finding multiple delays in complicated lens systems or looking for unknown lensing systems in the Fermi catalog of AGNs.

References

All Figures