IDENTIFYING ANOMALIES IN GRAVITATIONAL LENS TIME DELAYS

Gravitational lensing has become a valuable probe of dark matter substructure in distant galaxies (see Section B.8 of Kochanek et al. 2006b, and references therein). A growing body of evidence suggests that anomalous flux ratios, which are observed in many four-image quasar lenses (Metcalf & Zhao 2002; Keeton et al. 2003, 2005), can be explained if a few percent of the projected mass within each lens galaxy's Einstein angle (which sets the spatial scale for lensing) is contained in cold dark matter (CDM) "clumps" (Mao & Schneider 1998; Metcalf & Madau 2001; Chiba 2002; Dalal & Kochanek 2002; Kochanek & Dalal 2004). It is then natural to ask whether it is possible to use lensing to measure properties beyond the mean substructure mass fraction. Of particular interest is the clump mass function, because the discrepancy between the observed number of Local Group satellite galaxies and the theoretically predicted abundance of dark matter clumps varies strongly with mass (e.g., Klypin et al. 1999; Moore et al. 1999; Strigari et al. 2007). Unfortunately, it is difficult to constrain the mass function with anomalous flux ratios, because there is a degeneracy such that a small clump near a lensed image (in projected distance) can produce the same flux perturbation as a large clump farther away (see Dalal & Kochanek 2002).

One possible solution to this problem is to measure flux ratios at multiple wavelengths. Quasar emission regions are believed to have different sizes at different wavelengths. Roughly speaking, only clumps with Einstein radii (Einstein angles translated into units of length) larger than the size of the source can produce flux ratio anomalies (e.g., Dobler & Keeton 2006); clumps with small Einstein radii act collectively as a smooth mass component. Since the Einstein radius depends on mass, the source size effectively translates into a mass threshold. The net effect is that multi-wavelength observations provide a way to study clumps with different mass ranges.⁴ For example, continuum (Woźniak et al. 2000) and broad-line (Richards et al. 2004; Keeton et al. 2006) optical emission regions in quasars can be perturbed by objects of stellar mass and larger. Radio emission regions, by contrast, are large enough that only objects with masses ≳10⁶M_☉ are important (Dobler & Keeton 2006). If a flux-ratio anomaly is observed at optical wavelengths but not at radio wavelengths, it is likely that microlensing by stars in the lens galaxy is the culprit. To avoid contamination by microlensing, Dalal & Kochanek (2002) focused on radio anomalies in their study. Infrared emission regions are intermediate in size, so observations in this band provide a way to probe objects with mass scales between stars and clumps. Recent mid-IR observations with the Spitzer (Poindexter et al. 2007a, 2007b) and Subaru (Chiba et al. 2005; Minezaki et al. 2009) telescopes have achieved the spatial resolution necessary for strong lensing studies. In particular, Chiba et al. (2005) found evidence for subhalos of ∼10⁴M_☉ in the lens system B1422+231, and possibly in PG 1115+080 as well. In addition, Minezaki et al. (2009) recently found evidence of a clump with mass >10⁵M_☉ in the lens MG 0414+0534. A related way to study the substructure mass function is to compare continuum and broad-line emission in optical lenses (Moustakas & Metcalf 2003), which originate from the central accretion disk and an extended distribution of fast-moving clouds, respectively. These two regions differ in length scale by a factor of ∼100, whence their utility for substructure lensing. This technique has been employed to study the four-image lenses HE 0435 − 1223 (Wisotzki et al. 2003), SDSS J0924+0219 (Keeton et al. 2006; Eigenbrod et al. 2006), SDSS J1004+4112 (Richards et al. 2004), RX J1131 − 1231 (Sluse et al. 2007), and Q2237+0305 (Metcalf et al. 2004; Eigenbrod et al. 2008), along with several two-image lenses (Burud et al. 2002; Inada et al. 2005a, 2006; Sluse et al. 2008).

While the method of probing substructure with differential source size effects is promising, this approach does have some limitations. For one thing, the mass threshold defined by the source size and Einstein radius does not represent a sharp cutoff (see Dobler & Keeton 2006), so it is difficult to obtain tight constraints on the clump mass. In addition, it may not be possible to measure flux ratios at enough wavelengths to sample the mass function well. Finally, lensing constraints on the mass function will only be as accurate as the mapping of wavelength to source size. While the model of Shakura & Sunyaev (1973) has been quite successful in explaining the observed properties of accretion disks, and gives a relation between wavelength and source size, it does not apply to the full range of wavelengths relevant to strong lensing.

A second possible approach to constrain the substructure mass function is to use flux ratios in combination with other lens observables. Although they are smaller in amplitude than flux ratio perturbations, substructure effects on both image positions (Chen et al. 2007) and time delays (Keeton & Moustakas 2009) should be detectable. Combining different observables will be valuable because they depend on the lens potential in different ways: time delays, image positions, and flux ratios depend on the potential and its first and second derivatives, respectively. The different observables, in other words, contain different information about the substructure population. For example, Keeton & Moustakas (2009) showed that time delays are sensitive to the slope and dynamic range of the substructure mass function, which suggests that precise time delay measurements may provide a way to constrain these properties. Before we can use this method, we need to determine whether observed time delays differ from the predictions for smooth mass models in a way that may indicate the presence of substructure. Developing a method to identify such "time delay anomalies" is the focus of this paper.

In particular, we wish to understand the range of time delays that can be produced by reasonable smooth models. We can then classify any outlier as "anomalous" and use it as evidence of complexity in the lens potential. We must be careful when interpreting anomalies, however, because dark matter substructure may not be the only relevant source of complexity in the lens galaxy. Stars also constitute complex structure that is important when interpreting flux ratios (due to microlensing), but they have essentially no effect on time delays (Keeton & Moustakas 2009). This means that time delays should not depend on wavelength, so we do not need to distinguish between radio and optical time delays.⁵ Finally, extinction by dust in the lens galaxy can perturb flux ratios, but it will not have any effect on time delays, which are measured through flux variability alone and do not depend on color information.

The other main source of complexity we need to consider is the environment of the lens galaxy. Many lens galaxies lie in groups or clusters of galaxies (e.g., Momcheva et al. 2006; Auger et al. 2007). To lowest order, such an environment contributes a tidal shear to the lens potential (e.g., Keeton et al. 1997), and we therefore include shear in our analysis, but there may be higher-order terms that are non-negligible. Extreme examples of this situation include the lens B1608+656, which has two galaxies inside the Einstein angle (Koopmans & Fassnacht 1999; Koopmans et al. 2003; Suyu et al. 2009), and the lens SDSS J1004+4112, which is produced by a cluster of galaxies (Inada et al. 2003; Oguri et al. 2004). In general, time delay anomalies in close pairs of images potentially provide the strongest evidence of dark matter substructure, because environmental effects are fairly large-scale and should not produce dramatic differences between images separated by a distance much less than the Einstein angle. With time delay anomalies in image pairs that have larger separations, by contrast, we will need to take more care to consider environmental effects as well as substructure.

Our approach follows that of Keeton et al. (2003, 2005), who employed analytic flux-ratio relations that are generic for all lenses with fold and cusp configurations produced by smooth mass models. To be more specific, a fold lens contains a bright pair of images whose fluxes F_A and F_B should satisfy the relation

$$\begin{equation} R_{\rm fold}\equiv \frac{F_A - F_B}{F_A + F_B} \approx A_{\rm fold}\, d_1, \end{equation} \tag{ 1 }$$

where d₁ is the image separation and A_fold depends on properties of the lens potential. A cusp lens contains a triplet of bright images whose fluxes should satisfy the relation

$$\begin{equation} R_{\rm cusp}\equiv \frac{F_A - F_B + F_C}{F_A + F_B + F_C} \approx A_{\rm cusp}\, d_1^{2}, \end{equation} \tag{ 2 }$$

where d₁ is the distance between the closest pair of images and A_cusp depends on properties of the lens potential (see, e.g., Congdon et al. 2008). (Note that R_fold and R_cusp vanish in the limit d₁ → 0.) If these relations are violated, we may conclude that the mass distribution of the lens galaxy cannot be smooth and must contain additional structure, most likely in the form of CDM substructure or a complex lens environment. Determining in practice that a given lens is anomalous requires some care. Equations (1) and (2) show that R_fold and R_cusp increase with the distance between the images, so a non-zero value for one of these quantities does not automatically imply a flux ratio anomaly. Making such an identification would require knowledge of A_fold and A_cusp, which depend on the exact form of the lens potential. Since these quantities are not directly observable, Keeton et al. (2003, 2005) performed Monte Carlo simulations to determine distributions for R_fold and R_cusp at fixed d₁ using parameter values appropriate for a realistic population of lens galaxies. With this information, it is possible to determine the probability that an observed lens is anomalous and hence contains small-scale structure.

Oguri (2007) employed a method similar to that of Keeton et al. (2003, 2005) to study time delays, although his emphasis was different from ours. While both studies are based on time delay distributions, Oguri's emphasis was on determining the Hubble constant, whereas we focus on small-scale structure in lens galaxies. In fact, we explicitly remove the Hubble constant from our analysis by working with scaled time delays and time delay ratios. Another difference is that we concentrate our attention mainly on cusp and fold lenses, which are most relevant to the questions we seek to answer.

This paper is organized as follows. In Section 2 we consider how the differential time delay depends on the parameters of the lens potential, and on the shape of the associated (astroid) caustic. The results we present there motivate our Monte Carlo simulations, which we describe in Section 3. We apply our formalism to the 25 known four-image lenses in Section 4. Finally, we present our conclusions in Section 5.

In Congdon et al. (2008), we showed that the time delay between a close pair of images in a fold lens scales with the cube of the image separation, d₁, i.e., Δτ/τ₀ ≈ |h|d³₁/2 (see Equation (25) of Congdon et al. 2008, and pp. 190–191 of Schneider et al. 1992), where τ₀ is a cosmology-dependent scale factor (see Equation (6) below). The coefficient h comes from a Taylor expansion of the lens potential, ψ, and can be written as a particular third derivative: h ≡ ψ₂₂₂/6, where the subscript "2" indicates differentiation with respect to the coordinate in the image plane that corresponds to the direction perpendicular to the caustic in the source plane, and is evaluated at the point on the critical curve that serves as the coordinate origin (see Section 3.2 of Congdon et al. 2008). In this section, we study how the time delay for a fold pair depends on the lens potential and the distance between the fold point and the nearest cusp point. This is equivalent to studying the variation of h along the caustic, since, for fixed d₁, the time delay is given solely in terms of this coefficient. Because h < 0, we find it more convenient to work with its absolute value.

Most lens galaxies are of early type and have density profiles close to isothermal, i.e., the three-dimensional density scales with radius as ρ ∝ r⁻² (e.g., Treu & Koopmans 2004; Rusin & Kochanek 2005; Treu et al. 2006), so we compute |h| for a singular isothermal ellipsoid (SIE) lens. To determine an appropriate value for the ellipticity parameter e ≡ 1 − q, where q is the minor-to-major axis ratio, we turn to the observed galaxy samples of Bender et al. (1989), Jørgensen et al. (1995), and Saglia et al. (1993). These samples have mean ellipticities and dispersions of $(\bar{e}, \sigma _e) = (0.28, 0.15), (0.31, 0.18)$, and (0.30, 0.16), respectively. Note that these values measure the distribution of light rather than mass, so it is possible that the dark matter halo in which the galaxy presumably resides is rounder or flatter than the observed isophotes.

Figure 1 shows |h| at various points on the caustic, for ellipticities of 0.1, 0.3, and 0.5. We see that |h| remains roughly constant along the caustic for points away from the cusps. This suggests that lenses whose fold pairs have comparable separations will have similar time delays as well, at least for galaxies with similar ellipticities. Testing this prediction will require large samples of fold lenses, for which both the differential time delay between the fold pair and the ellipticity of the lens galaxy are known. For points near a cusp, |h| decreases rapidly and vanishes right at the cusp point (as it must; see Section 2 of Congdon et al. 2008). Finally, we see that |h| depends on the size of the caustic. It is not yet clear whether this merely reflects a simple correlation between |h| and e or is indicative of a more subtle relationship between caustic size and the time delay for a fold pair.

Zoom In Zoom Out Reset image size

In addition to ellipticity, Bender et al. (1989) and Saglia et al. (1993) find that many early-type galaxies have small departures from elliptical symmetry. "Disky" or "boxy" isophotes can be represented by an m = 4 multipole mode, with coefficient a₄ > 0 for disky isophotes and a₄ < 0 for boxy isophotes. In terms of the (two-dimensional) lens potential, the m = 4 mode for an isothermal galaxy can be written as (e.g., Keeton 2001)

$$\begin{equation} \psi (r,\theta)=-A_4 r \cos {4\theta }{.} \end{equation} \tag{ 3 }$$

The coefficient a₄ is defined in terms of surface brightness rather than potential, so we must convert from A₄ to a₄:

$$\begin{equation} a_4 = 15A_4 \sqrt{1-\epsilon } \,, \end{equation} \tag{ 4 }$$

where ≡ (1 − q²)/(1 + q²). For typical values of a₄ = -0.01, 0, 0.01, and 0.02 (e.g., Bender et al. 1989) and ellipticity e = 0.3, we find maximum |h|-values of 0.07, 0.07, 0.08, and 0.09, respectively. We therefore conclude that m = 4 multipole terms have fairly small effects on the time delay for a fold pair, regardless of whether the isophote is disky or boxy. We nevertheless include m = 4 modes in our simulations (Section 3) given that this is both simple and observationally motivated.

Since many lens galaxies lie within groups or clusters (e.g., Momcheva et al. 2006; Auger et al. 2007), a nonzero tidal shear is common. Using numerical simulations and semianalytic models of galaxy formation, Holder & Schechter (2003) find that shear can be described by a lognormal distribution with median γ = 0.05 and dispersion σ_γ = 0.2 dex.⁶ For such shear amplitudes, the caustic structure is qualitatively similar to what is seen in Figure 1, so we do not show it explicitly. We return to shear when we consider realistic lens populations in the following section.

For any given lens galaxy, we could estimate h (or more generally the smooth model time delays) directly, by fitting a lens model and inferring the lens potential. However, we would like to avoid any dependence on models and modeling as much as possible. We therefore elect to examine the full probability distributions for time delays given a realistic population of lens galaxies, using Monte Carlo methods to construct the time delay distributions.

We follow the approach of Keeton et al. (2003, 2005) and create lens galaxies with model parameters drawn from the galaxy samples of Bender et al. (1989), Jørgensen et al. (1995), and Saglia et al. (1993). Although the mean ellipticities and dispersions are similar for the three samples, the underlying galaxy populations are not the same. Jørgensen et al. (1995) and Saglia et al. (1993) include galaxies within clusters, while Bender et al. (1989) use bright, nearby galaxies to construct their sample. In addition, Bender et al. (1989) and Saglia et al. (1993) report a₄ measurements, while Jørgensen et al. (1995) do not. In order to determine the extent to which our conclusions depend on the input data, we carry out the numerical method described below separately for each of the three galaxy samples. To model the environment of a lens galaxy, we add tidal shear based on the simulations of Holder & Schechter (2003) that we discussed at the end of the previous section. For the sample of Jørgensen et al. (1995) we pick 2000 ellipticities from their observed distribution, and assign a random shear to each. For the samples of Bender et al. (1989) and Saglia et al. (1993), we use the actual (e, a₄) pairs for the observed galaxies (87 in the Bender et al. 1989 sample, and 54 in the Saglia et al. 1993 sample) in order to retain any correlation between the parameters; and for each galaxy we use 100 different realizations of the shear.

For each model lens potential, we use an updated version of the GRAVLENS software⁷ (Keeton 2001) to solve the lens equation numerically and obtain the image positions and time delays for a large set of random source positions that are uniformly distributed in the four-image region. Using a uniform distribution has several consequences. First, it means that each model lens potential is automatically weighted by its four-image cross section, which seems like the proper statistical approach.⁸ Second, it means that we neglect magnification bias. While lensing magnification bias is quite important when comparing statistical samples of four-image and two-image lenses, it is less dramatic within a sample of four-image lenses, and still less so within the subset of four-image lenses that have a particular image configuration. Generally speaking, including magnification bias would give more weight to sources that lie closer to the lensing caustic, which tend to produce shorter time delays, so it would shift the time delay distributions we derive to somewhat shorter values. The effect is not strong, however, and we checked that it does not change any of our conclusions about time delay anomalies. All told, our simulations based on the galaxy samples from Bender et al. (1989), Jørgensen et al. (1995), and Saglia et al. (1993) contain 1,267,555, 2,205,515, and 851,261 mock four-image lens systems, respectively.

The time delay of an image at angular position θ relative to an unlensed light ray from the true source with position β is given by

$$\begin{equation} \tau ({\mbox{\boldmath $\theta $}}) = {\tau _0} \left[ \frac{1}{2} |{\mbox{\boldmath $\theta $}}-{\mbox{\boldmath $\beta $}}|^2 - \psi ({\mbox{\boldmath $\theta $}}) \right], \end{equation} \tag{ 5 }$$

where the time scale is

$$\begin{equation} \tau _0 = \frac{1+z_L}{c} \frac{D_L D_S}{D_{\rm LS}}. \end{equation} \tag{ 6 }$$

Here D_L, D_S, and D_LS are the angular-diameter distances from the observer to lens, observer to source, and lens to source, respectively. For our mock lenses we focus on the dimensionless, scaled time delay

$$\begin{equation} \hat{\tau }\equiv \frac{\tau }{\tau _0 \theta _E^2}\,, \end{equation} \tag{ 7 }$$

where θ_E is the Einstein angle.⁹ The advantage of working with the scaled time delay is that our analysis of the mock lenses is independent of cosmology, and of the lens and source redshifts. It does mean that we need to convert all observed time delays from physical to dimensionless units before we can compare them to our mock lens catalog (see Section 4). For convenience, we quote distances (e.g., d₁ and d₂ below) in units of the Einstein angle as well.

Our analytic time delay relation (Equation (25) of Congdon et al. 2008) specifically applies to fold image pairs, so they are the most obvious targets for anomaly searches. However, we do not actually use the analytic relation to make predictions for our mock lenses, so it is reasonable to consider non-fold image pairs, as well as cusp and cross lenses, in our comparisons between observed lenses and Monte Carlo simulations. There are six distinct image pairs in a four-image lens. We are interested in the subset of pairs that contain one image of positive parity (which lies at a minimum of the time delay surface) and one of negative parity (which lies at a saddle point). In this way, we consider only adjacent images. The effect is to make our analysis local¹⁰ rather than global, which is appropriate for studying small-scale structure. There are four such mixed-parity image pairs in each four-image lens.

We characterize each image pair by the separation between the two images, d₁. We also use the distance to the next-nearest image, d₂,¹¹ because it helps specify the lens morphology: a fold lens has d₁ ≪ d₂ ∼ 1; a cusp lens has d₁ ∼ d₂ ≪ 1; a cross lens has d₁ ∼ d₂ ∼ 1. Another reason to use d₂ is that it allows us to take ratios of time delays (see below). While d₁ and d₂ technically define an image triplet, we will use the term "pair" to refer to the two images defined by d₁, since this is the image pair of primary interest.

We seek to understand the full range of time delays that can be produced by smooth lens models that are consistent with the positions of an observed image pair. We consider a mock lens image pair to "match" an observed pair if two conditions are satisfied. First, we require that the d₁ and d₂ values of the mock and observed pairs agree within a tolerance of ±0.05. This criterion is conservative in the sense that we are using minimal information from quantities that are observable and local; using additional information would only narrow the range of smooth models that are considered to be consistent with an observed image pair. The distance-based matching criterion does have one limitation, which is illustrated in Figure 2: it may allow not only mock pairs with the same parities as the observed images, but also pairs with the opposite parities. We could avoid this problem by adding information about the fourth image. We do not want to do that, however, because the fourth image is always "far" from the pair in question (i.e., 1–2 Einstein angles away); including it would make our analysis more sensitive to global properties of the lens potential and countermand our goal of using a local analysis to search for substructure. Instead, we add the image parities to the matching criterion.¹² Specifically, we insist that the parities of the three images defined by d₁ and d₂ are the same for an observed lens and its simulated counterpart. This is straightforward to do since the image parities are known for the mock lenses, and they can be determined unambiguously for observed lenses (by measuring time delays or just analyzing the image configuration; see, e.g., Saha & Williams 2003). Since the two images in each pair we consider have opposite parities (by construction), the parity cut represents only one additional bit of information.

Zoom In Zoom Out Reset image size

For each image pair we compute the differential scaled time delay between the two images, which we call Δt₁. For the second differential time delay we use the image pair characterized by the distance d₂, and we call this Δt₂. We then adopt the following sign convention. By construction, any image pair we consider will consist of a minimum image (positive parity, labeled M₁) and a saddle image (negative parity, labeled S₁). If the next-nearest image is a minimum (labeled M₂), we have a triplet consisting of a saddle flanked by two minima. We define the two differential time delays to be

$$\begin{equation} \Delta t_1 = \hat{\tau }(M_1) - \hat{\tau }(S_1) \quad {\rm and} \quad \Delta t_2 = \hat{\tau }(M_2) - \hat{\tau }(S_1). \end{equation} \tag{ 8 }$$

Since saddles have larger time delays than minima, both Δt₁ and Δt₂ are negative. In the case that the next-nearest image is a saddle (labeled S₂), we have a triplet consisting of a minimum flanked by two saddles. We then define the differential time delays to be

$$\begin{equation} \Delta t_1 = \hat{\tau }(S_1) - \hat{\tau }(M_1) \quad {\rm and} \quad \Delta t_2 = \hat{\tau }(S_2) - \hat{\tau }(M_1). \end{equation} \tag{ 9 }$$

Now both Δt₁ and Δt₂ are positive.

It is also useful to compute the ratio of the differential time delays for two image pairs. Time delay ratios are attractive because they do not depend on cosmology (τ₀ factors out), and they are largely immune to the radial-profile degeneracy (e.g., Kochanek 2002; Keeton & Moustakas 2009). Our sign convention ensures that Δt₁/Δt₂ is always positive since Δt₁ and Δt₂ have the same sign.

To reprise, we define the two time delays relative to the "middle" of the three images we use to determine d₁ and d₂. This means the scaled time delays can be either positive or negative, but the time delay ratio is always positive. Note that regardless of the sign, when labeling image pairs we always list the minimum first.

From the set of mock lenses that match an observed image pair, we construct histograms of the scaled time delay and the time delay ratio. These represent the range of values that can be produced by realistic smooth mass distributions. We can then use these predicted distributions to assess whether observed time delays are or are not consistent with lensing by a smooth mass distribution.

The number of observed four-image lenses has been steadily increasing in recent years. We focus on the 25 lens systems for which the lensed images are point-like (see Table 1). This restriction prevents us from using the system Q0047 − 2808 (Warren et al. 1996, 1999), as well as lenses from the SLACS survey¹³ (Bolton et al. 2006, 2008). However, this condition is necessary to ensure that the lensed source is compact, so that the variability timescale is short enough to make time delay measurements practical.

As noted above (see Equation (7)), we compute scaled time delays for our mock lenses, so we must scale each observed time delay by τ₀θ²_E in order to compare it with the mock lens sample. The Einstein angles and lens and source redshifts are listed in Table 1. We compute the angular diameter distances in τ₀ assuming cosmological parameters Ω_M = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹. The benefit of working with scaled time delays in our mock lens analysis is that we can present results for all observed lenses, regardless of whether any of their time delays are currently known. As more time delays are measured, it will be a simple matter to compare them with our compilation of predicted scaled time delays. In addition, it is possible to consider different cosmological models simply by recomputing the scale factor τ₀ and retranslating between observed and scaled time delays.

4.1. Time Delay Histograms

To get a sense of the range of time delays that can be produced by smooth lens models, we show histograms of the scaled time delay (Figures 3–6) and time delay ratio (Figures 7–10) for each image pair in the 25 known four-image lenses. There are many lenses whose time delays have not yet been measured, but the histograms for those cases are still pedagogically useful and provide a way to predict what the time delay should be if the lens in question is not anomalous. In this subsection, we discuss the general features of the histograms. In the following subsections, we analyze the lenses with known time delays and make predictions for the remainder.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We divide the observed lenses into three groups: folds (Figures 3, 4, 7, and 8), cusps (Figures 5 and 9), and crosses (Figures 6 and 10). We assign one of these canonical lens morphologies to each observed four-image lens according to the classification scheme of Keeton et al. (2003, 2005). For each lens, we define d*₁ to be the smallest value of the pairwise image separations. The lenses in each figure are arranged in rows such that d*₁ increases from top to bottom. Within each row, the panels are arranged such that d₁ increases from left to right. As a check for systematic effects, we plot histograms of the time delays and time delay ratios produced by Monte Carlo simulations using galaxy samples from Bender et al. (1989), Jørgensen et al. (1995), and Saglia et al. (1993; see Section 3).

We first consider histograms of the scaled time delays, which are shown in Figures 3–6. A negative time delay indicates that the middle image in the triplet defined by d₁ and d₂ has negative parity (Equation (8)), while a positive time delay indicates a middle image with positive parity (Equation (9)). The histograms are not highly sensitive to the input galaxy sample, especially in the tails (which are most important for identifying anomalies). The histograms are somewhat asymmetric with a longer tail on the side away from zero, which presumably reflects the fact that the histograms are bounded by zero by construction. Many of the histograms have a large width relative to the mean/median, which is not too surprising because we have been quite generous in matching mock lenses to observed lenses on the basis of minimal information (just d₁, d₂, and parity). This is consistent with our goal of being conservative in identifying time delay anomalies.

Next, we consider histograms of the time delay ratios, which are shown in Figures 7–10. These histograms are skewed to the right, presumably because the time delay ratios are positive and hence the histograms are bounded by zero on the left but unbounded on the right. The overall structure of the time delay ratio histograms does not vary much from one lens to another, or between image pairs of a given lens. This suggests that conclusions drawn from time delay ratios are not terribly sensitive to the lens morphology (fold, cusp, or cross), which may prove quite useful.

4.2. Identifying Time Delay Anomalies in Observed Lenses

4.3. Predictions for the Remaining Lenses

We have introduced a new method to use gravitational lens time delays to detect complex structure in the lens potential. The complexity may be associated with CDM substructure, in which case time delays offer the chance to learn more about the substructure population than is possible with lens flux ratios; or it may be associated with the lens environment, such as a group or cluster of galaxies surrounding the lens. The basic approach is to determine the range of time delays that can be produced by reasonable smooth lens models, so that we can identify outliers as being anomalous. To get a sense of how this could work, we first studied the dependence of the time delay between the close pair of images in a fold lens on the position of the source and the form of the lens potential. For a source near a fold point, we have found that the time delay remains approximately constant as the source moves along the caustic. For a lens modeled by an elliptical galaxy with m = 4 multipole perturbations and tidal shear, the time delay increases with ellipticity and shear, but is not very sensitive to m = 4 modes.

Using Monte Carlo simulations, we then constructed distributions of the time delays in four-image lenses. This approach can handle fold, cusp, and cross lenses, which comprise the three canonical four-image lens morphologies. By constructing a catalog of mock lenses based on observed populations of elliptical galaxies, we computed the range of time delays and time delay ratios that would be expected for a smooth lens potential (i.e., one with ellipticity, shear, and m = 4 multipoles). By comparing observed time delays with the predicted ranges, we have found time delay anomalies in the systems RX J0911+0551, RX J1131 − 1231, B1422+231, B1608+656, and WFI 2033− 4723. It is unlikely that these anomalies can be explained by errors in our assumed values of the Hubble constant or the slope of the density profile: the Hubble constant would have to be unreasonably high or the density profile surprisingly shallow in order to explain some of the apparent anomalies (we specifically discussed RX J0911+0551 and WFI 2033 − 4723), and neither possibility could explain RX J1131 − 1231. It is possible to further reduce sensitivity to the Hubble constant and density profile by working with time delay ratios, although in general we have found that time delay ratios have less power to reveal anomalies. Part of the problem is that there are fewer time delay ratios known than time delays themselves, and a large uncertainty in a given time delay will lead to a correspondingly large uncertainty in the ratio, making definitive conclusions difficult.

In general, anomalies between close pairs of images in fold and cusp lenses should provide the cleanest evidence of substructure. The cusp lens RX J1131 − 1231 contains such anomalies, which is consistent with conclusions based on more detailed modeling by Morgan et al. (2006) and Keeton & Moustakas (2009). The cusp lenses B1422+231 and RX J0911+0551 show evidence of time delay anomalies in the cusp triplet, but the large uncertainties in the measured time delays prevent firm conclusions at present. The only fold pair that is clearly anomalous is in B1608+656, but the peculiar nature of this system (with two lens galaxies inside the Einstein angle) makes it difficult to draw a definitive conclusion about substructure from our analysis. Among larger-separation image pairs we have found time delay anomalies in the fold lens WFI 2033 − 4723 and the cusp lens RX J0911+0551. Both lenses show evidence of a complex environment, but in order to distinguish between that and substructure as the origin of the time delay anomalies it will be necessary to precisely measure the time delays between the close images (the fold pair in WFI 2033 − 4723, and the cusp triplet in RX J0911+0551).

In the hope that the sample of observed precision time delays will continue to grow, we have predicted the time delays for all mixed-parity image pairs in all 25 known four-image lenses. As new time delays are measured, it will be a simple matter to compare them with our predicted confidence intervals to determine whether they are anomalous. If a lens galaxy contains complex structure, time delays should help reveal it. Flux ratios provide a powerful way to find small-scale structure, but they are not unique in this; time delays hold great promise for contributing to our understanding of the role played by dark matter in the universe, especially when combined with other lensing observables.

We thank Leonidas Moustakas and Ross Fadely for helpful conversations. We also thank the referee, Olaf Wucknitz, for several helpful suggestions that improved the manuscript. A.B.C. would also like to thank Mark Eichenlaub and Dennis Lam for their input. Part of this work was funded by NSF grant AST-0747311. A.B.C. is currently supported by an appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, administered by Oak Ridge Associated Universities through a contract with NASA.