THE EFFECTS OF HALO-TO-HALO VARIATION ON SUBSTRUCTURE LENSING

In the standard cold dark matter (ΛCDM) cosmology, structure forms hierarchically. Small dark matter (DM) halos form first and merge to form larger objects. The process of merging and dynamical evolution leaves remnants of smaller halos (called subhalos) within larger host halos. We refer to the population of subhalos as substructure. The presence of substructure in galaxy-sized DM halos is an important test of CDM. While the distribution and number of galaxies and clusters of galaxies are in agreement with the generic predictions of CDM, at sub-galactic scales numerical simulations predict an abundance of DM substructure which has not been detected by the presence of baryonic matter (stars and gas), a situation known as the "missing satellites problem" (Klypin et al. 1999; Moore et al. 1999; Diemand et al. 2007; Springel et al. 2008). Detecting DM subhalos and putting constraints on the substructure mass function is a key step toward the identification of the DM particle (see, e.g., Moustakas et al. 2009).

One of the most promising avenues to measuring substructure in galaxies is via strong gravitational lensing. As light from a distant object passes an intervening galaxy, it may be bent sufficiently for multiple images of the source to be observed. Strong lens systems with quasar sources and galaxy-sized lenses have been shown to be sensitive to satellite galaxy-sized DM subhalos in the lensing halo³ (Mao & Schneider 1998; Metcalf & Madau 2001; Metcalf & Zhao 2002; Dalal & Kochanek 2002; Chiba 2002; Kochanek & Dalal 2004). Substructure in the lensing halo may be observed as perturbations to the image magnifications (flux ratio anomalies), image positions (e.g., Metcalf & Madau 2001; Chen et al. 2007), and/or time delays in systems with time variable sources (Keeton & Moustakas 2009).

The most well-explored approach to substructure lensing is the study of anomalous flux ratios in multiply imaged systems. Anomalous flux ratios are sensitive to the projected mass fraction in substructure, f_sub, at projected radii where lensed images are present, ∼3% of the projected virial radius of the lens halo. Results, however, have not been definitive. Dalal & Kochanek (2002) find that f_sub = 2% with 0.6% < f_sub < 7% within 90% confidence intervals. This is marginally consistent with the results of numerical simulations. Metcalf & Amara (2010) also find a result which is generally consistent with numerical simulations.

Several studies have shown that f_sub is small (≲ 1%) but with large variance. For example, Zentner & Bullock (2003) used semi-analytic simulations of a large number of halos to argue that substructure mass fractions at projected radii of 3% of halo virial radii should be f_sub ≈ 0.6% on average and f_sub ≲ 1% for 95% of halos. Mao et al. (2004) found f_sub ≲ 0.5% using 12 halos at galaxy, group, and cluster masses. Xu et al. (2009) found 0.01% ⩽ f_sub ⩽ 0.7%, using six halos from the Aquarius simulation. Using 26 simulated halos, Zentner (2006) showed that f_sub < 1% along most lines of sight, but that there is a large variance due, in part, to the orientation of the halo with respect to the line of sight.

Other studies have suggested that lensing observations are inconsistent with expectations from numerical simulations. Using simple models for substructure, Chen (2009) found that observations support the presence of more substructure than expected from simulations. Macciò et al. (2006) and Amara et al. (2006) each compared observations to a single simulated halo and found more and greater magnification perturbations than can be accounted for in their simulated halos. In addition, the results of Dalal & Kochanek (2002) and Metcalf & Amara (2010) (which are consistent with simulations) use substructure models that make it possible that their results underestimate the value of f_sub.

Future, more definitive, results will require better characterizations of sources of error, more observational data, and better predictions for CDM substructure. In this work we focus on this last element. More specifically, as the initial density field is unknown, predictions for substructure populations in numerical simulations are specific to a particular realization of the simulated halo. Halo-to-halo variation, which should be present, is essential in determining the expected variance about each measurement. Quantifying this variance can allow more accurate error assessments and potentially help bring observational results into better agreement with theoretical predictions. In addition, understanding the variance of the substructure population may provide valuable insight into previously unexplored biases in the lensing samples.

In this work, we characterize the halo-to-halo variation of substructure populations in galaxy-sized DM halos and investigate how substructure properties correlate with other properties of DM halos. We study the effects of halo-to-halo variation to mock lensing observations, focusing on a simple statistic that quantifies violations of the cusp relation. We employ a semi-analytic method to generate a large number of substructure populations. We describe this method in Section 2 and describe basic properties of halo substructure in Section 3. We use our substructure models to generate strong-lensing simulations in Section 4. We analyze our simulations and compare the results to observations in Section 5, we discuss caveats in Section 6, and we conclude in Section 7. Throughout this paper, we work in a flat cosmological model with Ω_dm = 0.228, Ω_b = 0.0227, h = 0.71, and σ₈ = 0.81 as favored by the five-year Wilkinson Microwave Anisotropy Probe results (Komatsu et al. 2009).

In order to generate a large statistical sample of host halos, we use a semi-analytic technique (Zentner et al. 2005). This technique allows us to approximate the hierarchical assembly of a DM halo analytically and thus avoid the prohibitive computational cost associated with numerous, high-resolution numerical simulations (see also Taylor & Babul 2001, 2004).

The semi-analytic approach employs a variation of the extended Press–Schechter formalism to generate merger histories of host DM halos (Bond et al. 1991; Lacey & Cole 1993; Somerville & Kolatt 1999; Zentner 2007) and an approximate analytic approach to the dynamical evolution of subhalos in an evolving host halo. The details of this method are given in Zentner et al. (2005), including a number of comparisons to self-consistent numerical simulations. This technique has been tested for subhalos of mass M > 10⁻⁴M_h, where M_h is the mass of the host halo, and produces subhalo abundances and spatial distributions in good agreement with N-body simulations. Though this approach is approximate, it allows us to compute a large number of realizations of a host halo at a subhalo resolution which is inaccessible at present with state-of-the-art numerical simulations. In addition, this procedure can be used to produce physically motivated extrapolations to numerical simulations (see Koushiappas & Zentner 2006; Koushiappas et al. 2004, 2010).

We generate 200 realizations of a host halo of mass M_h = 1.26 × 10¹² h⁻¹ M_☉. We track the accretion and dynamical evolution of subhalos of mass greater than 10⁻⁵M_h. The mean properties of substructure in the sample of realizations is in agreement with numerical simulations (see Figure 1 in Koushiappas et al. 2010). The average number of subhalos with V_max > 4 km s⁻¹ within the inner 200 kpc of the halo is 2481 with a 68 percentile range of [1964–3007], which is consistent with the 2469 subhalos found within the same radius of a numerical simulation of a similar-size halo (Kuhlen et al. 2008). In addition, the scatter in the number of subhalos from halo-to-halo variation is consistent with a Poisson scatter added in quadrature to an intrinsic scatter of about 20% (as derived from recent numerical simulations; Boylan-Kolchin et al. 2010).

Each realization represents a possible subhalo population within the host halo. The distribution, structure, and properties of the subhalo population differ between realizations due to the statistical nature of the initial density field and the complexity of any individual subhalo's orbit. The output of each realization contains the concentration of the host and a list of all subhalos, including their evolved structural parameters—total bound mass, scale radius, and tidal radius—and their locations within the host (see Zentner et al. 2005). The individual merger history of each host influences the abundance and properties of its substructure population.

We construct a two-dimensional projection of the subhalo distribution by projecting the three-dimensional density along an arbitrary direction. The three-dimensional density profile of each subhalo is described by a Navarro–Frenk–White profile (NFW; Navarro et al. 1997), ρ(r) = ρ_s/[x(1 + x)²], where x = r/r_s, r_s is the scale radius, and r is the three-dimensional distance from the halo center. Each subhalo has a projected distance from the center of the halo, which we call R. For simplicity, we define a dimensionless projected distance for each subhalo as $\tilde{R} = R / R_{\rm h}$, where R_h is the virial radius of the host halo. This two-dimensional projection of the subhalo population is used in the mock lensing simulations we describe below. We also utilize a dimensionless subhalo mass $\tilde{M}$ defined such that the mass of a subhalo is $M= \tilde{M} M_{\rm h}$, where M_h is the mass of the host halo.

We can glean information about the variation in the distribution of substructure by calculating the subhalo mass fraction within a projected radius. The projected mass fraction within the inner $\tilde{R}$ fraction of the virial radius of the host in a given subhalo mass bin is

$$\begin{equation} f_{\mathrm{sub}}({<}\tilde{R}) = \frac{1}{\tilde{M}_{\rm h}({<}\tilde{R})} \sum _i \tilde{M}_i(\tilde{R}_i < \tilde{R}), \end{equation} \tag{ 1 }$$

where

$$\begin{equation} \tilde{M}_h({<}\tilde{R}) = \frac{c^2}{\mathcal {F}(c)} \int _0^{\tilde{R}} \tilde{R}^{\prime }\int _0^{{\ell _{\mathrm{max}}}} \frac{1}{\tilde{r} [ 1 + \tilde{r} ]^2} \, \,d \tilde{\ell } \,\, d\tilde{R}^{\prime } \end{equation} \tag{ 2 }$$

is the projected mass within the inner $\tilde{R}$ fraction of the virial mass of the host halo. In Equation (1), the sum over subhalo masses is over all subhalos i whose projected radial position $\tilde{R}_i$ is less than the projected host position $\tilde{R}$, independent of the projected size of the subhalo. In Equation (2), $\tilde{r} = c r / R_{\rm h} = \sqrt{\vphantom{A^A}\smash{\hbox{${ \tilde{R}^{\prime 2} + \tilde{\ell }^2}$}}}$, the upper integral in the line-of-sight element is $\tilde{\ell }_{\mathrm{max}} = \sqrt{\vphantom{A^A}\smash{\hbox{${1 - \tilde{R}^{\prime ^2}}$}}}$, and $\mathcal {F}(c) = \ln (1 + c) - c /(1 + c)$, where c is the concentration of the host halo. Note that values of f_sub ⩾ 1 are possible. This may be the result of two effects: $\tilde{M}_{\mathrm{h}} ({<} \tilde{R})$ from Equation (1) only refers to the smooth mass of the host halo and not the total mass, and we use the full subhalo mass even when most of the subhalo's bound mass falls outside the relevant region.

In the top panel of Figure 1 we show the projected cumulative substructure mass fraction, f_sub, in the complete sample of 200 realizations as a function of radius. The median value of f_sub at $\tilde{R}=1$ is 9% for subhalos with mass greater than 10⁻⁵M_h, with a 95 percentile in the range [2–17]%. Note that the spread about the median increases as $\tilde{R}$ increases. This is due to the fact that the outer regions of halos are dominated by recently accreted subhalos, while the subhalo population of the inner regions is old and highly evolved. The recent accretion of subhalos is an outcome of the individual accretion history of each host halo, and the large spread about the median in the outer regions mirrors the range of recent accretion histories.

Zoom In Zoom Out Reset image size

Relevant to lensing studies is the projected substructure mass fraction within the inner 3% of the virial radius of the host halo. The bottom panel in Figure 1 shows f_sub in the inner regions. At $\tilde{R} = 0.03$, the median value of f_sub is 0.25% with a 95 percentile in the range [0–1]%, meaning that only 2.5% of the realizations have a projected substructure mass fraction greater than 1% in the inner 3% of the virial mass of the host halos.

In order to explore the relationship between the host halo properties and the substructure lensing effects, we distribute realizations into three groups. In the first group, which we call C₀, we include host halos with a concentration c (see Equation (2)) within the 68% range of the mean concentration of the complete sample, 〈c〉 = 9.7. This includes host halos with concentrations in the range 6.7 ⩽ c ⩽ 13.4. We then collect the upper 16% of the concentration range (c > 13.4) in group C₊ and the lower 16% range (c < 6.7) in group C₋. As they form earlier, subhalos present in high-concentration hosts are expected to be fewer and with more concentrated subhalo density profiles relative to subhalos present in low-concentration hosts. The inverse relation between substructure abundance and halo concentration has also been seen in the simulations of Ishiyama et al. (2008) and Ishiyama et al. (2009), where they measure the abundance of subhalos around DM halos of ∼10¹² M_☉ and larger. Grouping the hosts (and their substructure populations) allows us to investigate the effects of host concentration on the expected lensing signal from its substructure.

The cumulative substructure mass fraction for halos with mass greater than 10⁻⁵ the mass of the host is also dependent on the concentration of the host halo. In the top panel of Figure 2 we show the median and 68 and 95 percentiles of the distribution of mass fractions found in our sample. Host halos with concentration in the 68 percentile about the median concentration of all realizations (middle figure) have a substructure mass fraction which is of order ∼10%, and a 68 percentile of ∼[5–20]% (the 95 percentile range is ∼[3–35]%). Higher concentration hosts have a smaller fraction (of order ∼6% and a 68 percentile range of ∼[3–10]%, and a 95 percentile range of ∼[2–25]%), while low-concentration hosts have a higher substructure mass fraction (median of ∼22%, and a 68 percentile range of [9–50]% and a 95 percentile of ∼[4–55]%). Systems of low concentration (and subsequently large values of f_sub) are systems which typically contain recent large mergers. Such systems may appear as interacting galaxy pairs or groups in observations.

Zoom In Zoom Out Reset image size

In the bottom panel of Figure 2 we show the projected mass fraction within the inner 3% of the projected virial radius (most relevant to lensing studies), for a number of mass bins.⁴ We find that the abundance of subhalos with masses greater than ∼ a few times 10⁻⁴ the mass of the host is negligibly small: even though the upper 95 percentile of the mass fraction increases with mass in the panels, the median and the upper 68 percentile drop to zero (as does the number abundance) as mass increases. Therefore, large subhalos contribute little to projected mass fractions on average, although, in the rare cases in which a large subhalo is present, the projected mass fraction is correspondingly large. This deficit of large subhalos is independent of the host halo concentration. However, on smaller scales (masses of order $\tilde{M} \sim 10^{-4}\hbox{ to }10^{-5}$), the effects of the host concentration are prominent. For low-concentration host halos, the median mass fraction is f_sub ≈ 0.3% with an upper 68 percentile of f_sub ≈ 0.8% (95 percentile is f_sub ≈ 1.3%), within the projected inner 3% of the virial radius. In contrast, the high-concentration host halos have a median of f_sub ≈ 0.035% and with an upper 95 percentile of less than 0.2%.

The origin of these differences between low- and high-concentration halos is the same as alluded to earlier. High-concentration halos generally form earlier than low-concentration counterparts, undergoing the majority of their mergers at higher redshift. As a consequence, the merging subhalos evolve for a relatively longer period within the dense environments of their hosts in high-concentration halos. This allows more time for orbital decay and mass loss. The increased propensity to lose mass in high-concentration halo is compounded by the fact that tidal forces are moderately stronger in high-concentration host halos at fixed, total bound mass. These effects, working together, render the survival rate of mass bound to halo substructure within high-concentration hosts smaller than in low-concentration hosts. The opposite is true of substructure in low-concentration host halos, so low-concentration host halos tend to have higher substructure mass fractions. These concentration-dependent effects are consistent with those found and discussed in Zentner et al. (2005) and Wechsler et al. (2006).

In Figure 3, we summarize these effects by showing the projected substructure mass fraction as a function of radius for the three concentration groups of host halos and three representative subhalo mass bins. It is apparent that host halos with low concentrations have a higher overall projected substructure mass fraction than halos with high concentration (see also the top panel of Figure 2). In addition, massive subhalos (subhalos with mass greater than 10⁻⁴ the mass of the host) tend to a median substructure mass fraction close to zero. However, the upper 68 percentile range is increasing relative to the 68 percentile of lower mass subhalos. For example, for masses ∼10⁻⁴ of the mass of the host, the upper 68 percentile implies that up to 16% of the low-concentration halos have substructure mass fractions which are greater than 3% (with 2% of the low-concentration halos having a fraction greater than ∼7.5%). The corresponding fractions for subhalos with mass ∼10⁻⁵ the mass of the host are 1.2% and 3% (for the top 16 percentile and 2 percentile of hosts, respectively).

Zoom In Zoom Out Reset image size

In addition, in high-concentration hosts, the 68 percentile range increases with projected radius. The inner regions of these halos have evolved for a considerable amount of time, while the outer parts reflect more recent accretion. On the other hand, in host halos with low concentrations, there has been, on average, less time for halo substructure to evolve. Moreover, recent accretion of subhalos with either highly elongated orbits or with orbits that decay rapidly due to dynamical friction can deposit large amounts of substructure at either small or large halo-centric radii. The consequence of this is that the range of substructure mass fractions extends to much higher values at all mass bins in low-concentration hosts.

We create mock lens systems with characteristics typical of observed lens systems and include the substructure populations derived in the semi-analytic catalogs (see the previous section). We choose a lensing galaxy at z = 0.5 and a source quasar at z = 2 and model the lensing halo as a singular isothermal ellipsoid with projected surface density,

$$\begin{equation} \kappa (\xi) = \frac{\Sigma (\xi)}{\Sigma _{\rm crit}} = \frac{b}{2 \xi }, \end{equation} \tag{ 3 }$$

where ξ is the elliptical coordinate, ξ² = x² + y²/q², and Σ_crit is the critical surface density for lensing. In the case of a spherical lens, q = 1 and b is the Einstein radius of the model. In addition to the smooth lens model, an external shear, γ, is applied to account for nearby structure, such as a group of galaxies. θ_γ is the angle between the shear and the major axis of the smooth lens model.

The host halo in the semi-analytic simulations is based on a Milky Way sized DM halo. It differs from lensing halos in a few respects. First of all, it has no baryons. The density profile of the Milky Way in DM only simulations matches an NFW profile. This profile is shallower than isothermal in the inner portions and steeper than isothermal in the outer portions. Real galaxies, however, contain baryons, and observations have shown that lensing halos (baryons + dark matter) follow a roughly isothermal profile near the halo center, where strong lensing is sensitive to substructure (Treu et al. 2006). In addition to the differences in density profile, we might expect that models with baryons would have less substructure than our semi-analytic simulations, as the presence of the host galaxy makes tidal forces stronger. On the other hand, this effect may be small as subhalos that are projected close to the halo center may be located at considerably larger physical distances.

The host halo for the substructure catalogs also differs with respect to typical lensing halos in a few, more subtle ways: the total mass, the formation, and merger history, and the epoch of observation. The host halo is a Milky Way sized object, while lensing halos are mostly those which contain large, elliptical galaxies. We might expect that realistic large, lensing halos might have somewhat more substructure than Milky Way sized galaxies, as larger halos form later (Gao et al. 2004; Giocoli et al. 2010). Gao et al. (2004), Zentner et al. (2005), and van den Bosch et al. (2005) suggest that the number of substructures scales with host mass as ${\it dN}/d\ln \tilde{M} \propto M_{\rm h}^{0.08}$, which would correspond to only a 15% increase in substructure compared to our simulations. A similar effect is possible as we use substructure catalogs at z = 0 instead of at z ∼ [0.5–1] where lensing galaxies are found, since larger fractions of substructure might be expected at earlier epochs (Madau et al. 2008). These effects, however, are likely to be significantly smaller than the effect of halo-to-halo variation, as shown by Xu et al. (2009) and Zentner et al. (2005).

Overall, the effect of baryons on halo structure remains an important question, and the optimal method of placing DM substructure in a lensing halo remains unclear. It might be possible to construct a sophisticated model that would connect the properties of the host halo to those of the lens halo, but it is neither relevant nor necessary for a general study, as is this work. Further, treating the host halo and the lensing halo rather independently is logical because strong lenses have significant surface density contributions from baryonic components of galaxies and because lensing halos are described in the data well by isothermal models with a narrow range of dispersion (Treu et al. 2006). In the absence of a model for these effects and all the effects discussed above and in order to be conservative in our calculations, we preserve (1) the projected distance from the center of the halo (i.e., $\tilde{R}$ from the previous section) and (2) the fractional mass of substructure from our semi-analytic model in our lensing calculations. In addition, the host halo mass, M_h = 1.8 × 10¹² M_☉, is scaled to that of a more typical lensing halo, M_lens = 10¹³ M_☉. Thus, the mass of each subhalo is M_{lens, sub} = (M_lens/M_h)M_sub and all scale lengths are x_{lens, sub} = (M_lens/M_h)^1/3x_sub, where x denotes either the tidal radius or the scale radius, and both $\tilde{R}_{\rm lens} = \tilde{R}_{\rm }$ and f_{lens, sub} = f_sub.

Using this model, the Einstein radius of a singular isothermal sphere (SIS) lensing halo at z = 0.5 can be estimated from Equation (3), b = 0586. This is generally consistent with observed lenses which have Einstein radii ∼1'' (e.g., Browne et al. 2003). It is also a reasonable value for any likely lensing halo. The Einstein radius of any SIS can be calculated from b = 4π(σ/c)²D_ls/D_os, where D_ls is the angular diameter distance from the lens to the source, D_os is the angular diameter distance from the observer to the source, and σ is the one-dimensional velocity dispersion. We expect lensing halos to have σ = [200–350] km s⁻¹ (Treu et al. 2006), which corresponds with b = [0.366–1.120]''. In Table 1, we show the model parameters for the lens halo and environment. We refer to these as the macromodel parameters.

Table 1. Simulation Parameters

Name	b	q	γ	θγ (rad)
M1	0586	0.7	0.2	π/3
M2	0586	0.7	0.1	π/6
M3	0586	0.9	0.2	π/3

Download table as:  ASCIITypeset image

We model subhalos as truncated NFW profiles. There is no closed form to calculate the image position deflections from NFW profiles with a sharp cutoff at the truncation radius. We use one of the modified, truncated forms provided in Baltz et al. (2009), where

$$\begin{equation} \rho (r)=\frac{\rho _s}{x(1+x)^2} \frac{\tau ^2}{\tau ^2+x^2}, \end{equation} \tag{ 4 }$$

where τ = r_t/r_s and r_t is the truncation radius. This form falls off as fast as r⁻⁵ beyond the truncation radius and the relevant lensing quantities can be calculated simply.

Using the combination of the three different macromodels and the set of 200 substructure realizations, we perform lensing simulations using an inverse ray-tracing code. The position of the source, $\vec{\beta }$, and the images, $\vec{\theta }$, are related by the lens equation,

$$\begin{equation} \vec{\beta } = \vec{\theta } - \nabla \phi (\vec{\theta }), \end{equation} \tag{ 5 }$$

where ϕ is the lensing potential. In the case of circular symmetry in the lens, ∇ϕ_R∝M(< R), where M(< R) is the mass of the lens enclosed by a cylinder of radius R. The lens equation is multi-valued and several image positions may correspond to a single source position. Thus, it is solved most simply by establishing a grid of image positions and calculating the source position for each point on the grid. For our simulations, the resolution on the image plane is 0.469 mas.

Every set of three points in the image plane defines a triangle and corresponds to a triangle in the source plane. The magnification of an image can be calculated by the ratio of the area of the image plane triangle to the area of the source plane triangle. A common assumption—one also used in this work—is that the angular size of the source is much smaller than that of the DM substructure. We approximate the source by a point. However, for high magnification images, the area of the triangle on the source plane may be smaller than observed source sizes, as discussed in Section 6.

For each macromodel and each substructure realization, we sample the source plane 50,000 times, creating ∼10, 000 mock four-image lens systems. We confine our investigation to four image systems, since those have a sufficient number of observables and have been proved to be most useful for substructure studies.

A parameter that can be used to characterize deviations from a simple, smooth lens is the signed sum of the image fluxes from any three adjacent images in a four-image system. We work with a scaled version of such a sum,

$$\begin{equation} R_{\rm cusp} \equiv \frac{| \mu _1 + \mu _2 + \mu _3|}{|\mu _1| + |\mu _2| + |\mu _3|} = \frac{| F_1 + F_2 + F_3|}{|F_1| + |F_2| + |F_3|}, \end{equation} \tag{ 6 }$$

where μ_i is the magnification and F_i is the observed flux, and both are signed quantities which indicate the image parities. In practice, determining the image parities is unnecessary; in a triplet of images, the outer two images have one image parity, and the middle image has the opposite parity. The ideal cusp relation has R_cusp = 0 in cases where the lens potential is smooth and the opening angle Δθ spanned by the images from a vertex at the lens center approaches zero. An illustrative example of a mock lens system, providing a clear definition of Δθ, is shown in Figure 4. In practice, we take the sum in Equation (6) over the three images that yield the smallest opening angle. The cusp relation is the measure used in many anomalous flux ratio studies (e.g., Macciò et al. 2006; Amara et al. 2006; Metcalf & Amara 2010).

Zoom In Zoom Out Reset image size

Table 2 lists R_cusp values for systems with a single lens galaxy and four distinct, point-like images, observed in the radio or mid-infrared (noted by "IR"). This excludes B1608+656 (Koopmans & Fassnacht 1999; Fassnacht et al. 2002), B1127+385 (Koopmans et al. 1999), and B1359+154 (Myers et al. 1999; Rusin et al. 2001) which have multiple lensing galaxies. Limiting ourselves to distinct images also excludes MG 2016+112 (Garrett et al. 1994, 1996; Koopmans et al. 2002), MG 0751+2716 (Lehar et al. 1997; Tonry & Kochanek 1999), and B1938+666 (King et al. 1997; Tonry & Kochanek 2000), in which merging pairs or arcs are seen. We specifically include the lens systems which have observable satellite galaxies near the lens, B2045+265 and MG 0414+053, which each have one. Not all of the systems in Table 2 are nearly cusp catastrophes—i.e., a cusp configuration where Δθ → 0—but in all cases R_cusp is calculated from the three images with the smallest opening angle.

Radio data have been traditionally used in lensing studies of substructure as it is free of microlensing contamination from stars in lensing galaxies. This is because the size of the quasar source in the radio (∼[1–10] pc) is larger than the scale size of perturbations from stars (e.g., Mao & Schneider 1998) and from subsolar mass DM halos (Chen & Koushiappas 2010). The size of the source relative to the perturbations from DM substructure, on the other hand, is very small. Therefore, in our simulations we approximate the source as a point. Mid-infrared observations are probably the best flux measurements for substructure lensing, as the flux measurements are additionally free from extinction effects, as well as from microlensing effects (Minezaki et al. 2004; Chiba et al. 2005). Only five lens systems, however, have been observed in the mid-IR (see Table 2).

Figure 5 shows the values of R_cusp for macromodels M1 (blue dots), M2 (red dots), and M3 (green dots) without including the effects of substructure. The observed R_cusp values from Table 2 are shown as black crosses and circles, where mid-IR values are used in place of radio values where both are available. Small opening angles correspond to small cusp relation values, as is expected. At larger opening angles, for each model the values of R_cusp bifurcate, corresponding to source positions along the minor axis (smaller R_cusp values) or the major axis (larger R_cusp values) of the lensing halo. As expected, values of R_cusp depend on the macromodel parameters (Keeton et al. 2003). At the smallest of opening angles, R_cusp → 0, regardless of macromodel. At opening angles, ≲ 90°, the observed values of R_cusp clearly deviate from what may be expected from smooth models.

Zoom In Zoom Out Reset image size

The effect of substructure on the cusp relation is shown in Figure 6. The no substructure results from Figure 5 are reproduced here in red and orange shaded areas. The borders of the colored areas are contours of equal probability, and the red and orange patches encompass the areas where 99% and 100% of the cusp relation values are found. These areas are well contained within a small portion of the R_cusp–Δθ parameter space. Including the effect of substructure gives values which are shown by the dark gray and light gray shaded areas. The dark gray area, which encompasses the area where 99% of the values are found, overlaps significantly with the region in which the no substructure values lie. And where the substructure-included values do extend beyond this area, they mostly lie fairly close to the no substructure values. In these cases, the deviations from the no substructure case are not due to small perturbations to individual images but perturbations that effectively change the macromodel of the lens systems—e.g., shifting the lens center or changing the overall ellipticity of the lens. The light gray region shows where the final 1% of R_cusp values lie for the substructure-included systems. Here, some large deviations from the no substructure case can be seen. In particular, some small number of lenses show R_cusp ∼ 1 at opening angles less than 90°. So large deviations from the cusp relation due to the influence of substructure are possible but are infrequent.

Zoom In Zoom Out Reset image size

We bin the cusp relation values by opening angle and focus on three bins that contain observed cusp systems: 30° ⩽ Δθ ⩽ 40°, 70° ⩽ Δθ ⩽ 80°, and 100° ⩽ Δθ ⩽ 110°. When comparing our mock lens systems to observed R_cusp values, we find that the probability of high R_cusp values remains small. A summary is presented in Table 3. Here we use all 200 substructure realizations and estimate the cumulative probability of observing R_cusp values that are greater than provided by smooth models, P(R_cusp > RS), where RS is the R_cusp value in the no substructure case where P(R_cusp > RS|nosubhalos) = 0. For mock lens systems using the substructure realizations, P(R_cusp > RS) ≈ 10⁻². In addition, the cumulative probability of observing R_cusp values that are greater than the observed values, P(R_cusp > RO) is estimated, where RO is the largest R_cusp value observed for a given bin of opening angle. For mock lens systems using the substructure realizations, P(R_cusp > RO) ≈ 10⁻³.

In Figure 7, we show the distribution of R_cusp values in a given bin of opening angle. Here, the distribution of values for mock lenses with no substructure cluster near small R_cusp values (shown as cyan bars). On the other hand, the observed values are found at large R_cusp values (cross-hatched, magenta bars). The distribution of values for mock lenses with substructure are shown in red bars and black-lined bars. Substructure induces high-R_cusp tails which extend to (and, in some cases, past) the observed values of the cusp relation. In all three opening angle bins, these tails contain only a small fraction of the lens systems (a few percent). Interestingly, the high-R_cusp tails extend to larger values of R_cusp for smaller opening angles than for larger opening angles. The substructure realizations are split by the concentration of the host galaxy, as discussed in Section 2, plotting only the realizations that constitute the smallest 16% of host concentrations (red bars) and the realizations that make up the largest 16% (black-lined bars). We see that the correlation between host concentration and amount of substructure persists into observable lensing results. The hosts with the smallest concentrations have the largest amount of substructure and larger values of R_cusp when compared to the hosts with the largest concentrations. This is best seen in the bins of 70°–80° and 100°–110°, where the number of mock lens systems number in the thousands. In the smallest opening angle bin of 30°–40°, the number of mock lens systems is less than 20 per substructure realization and differences based on halo concentration are more difficult to discern.

Zoom In Zoom Out Reset image size

We test which parameterizations of the substructure content of lensing galaxies are most sensitive to anomalous fluxes and are most useful for comparing observed lens systems to the results of simulations. First, we explicitly tailor the definition of f_sub from Section 2, Equation (1) to match our macromodel conditions:

$$\begin{equation} f_{\rm sub} ({<}2b) = \frac{1}{M_{\mathrm{lens}} (R< 2b)} \sum _i M_{i} (R_{i}< 2b), \end{equation} \tag{ 7 }$$

where M_i(R_i < 2b) is a subhalo with mass M that is found to be projected within twice the Einstein radius of the lens b (within which all the images using macromodels M1, M2, and M3 have been found) and M_lens(R < 2b) is the mass of the lensing galaxy within 2b. In this equation, we have the freedom to define what it means for a subhalo to be projected within R = 2b. As discussed previously, the substructure in our catalogs have density profiles which are described by a truncated NFW profile, which has three parameters: the concentration, the bound mass, and the tidal radius. We could choose to include any subhalo which is projected such that some portion of the area within R = 2b falls within the tidal radius of the subhalo. This parameterization is supported by Rozo et al. (2006), who suggest that the most important DM halos for magnification perturbations are those within a tidal radius of an image. Another possibility is that a single subhalo could be projected such that some portion of the area within R = 2b falls within the scale radius of the subhalo. A final possibility is that a single subhalo could be projected such that the center of the subhalo falls within the area defined by R = 2b. This is the parameterization used in previous studies, when analyzing the outputs of DM simulations (Mao et al. 2004; Xu et al. 2009). The list of different parameterizations, then, is as follows: (1) include subhalos whose tidal radii fall within 2b, (2) include subhalos whose scale radii fall within that area, and (3) include subhalos whose halo centers fall within that area.

In Figure 8, the substructure realizations are split into low f_sub and high f_sub bins using the three criteria. Note that some realizations have no subhalos which satisfy the criterion, i.e., f_sub = 0: 4 out of 200 for case 1, 45 out of 200 for case 2, and 58 out of 200 for case 3. In addition, note that values of f_sub ⩾ 1 are possible. This may be the result of two effects: M_lens(R < 2b) from Equation (7) only refers to the smooth mass of the host halo and not the total mass, and we use the full subhalo mass even when most of the subhalo's bound mass falls outside the relevant region. Using all three criteria, we can see that substructure realizations with high f_sub values also have larger and more frequent violations of the cusp relation, when compared to low f_sub realizations. The differences between low and high values of f_sub are greatest for values of f_sub calculated via halo centers (case 3). Here, the median f_sub value is f_sub = 0.13%, and substructure realizations with f_sub greater than the median have an appreciable but small (∼10⁻²) fraction of large R_cusp values while those with f_sub smaller than the median may have no large cusp relation values.

Zoom In Zoom Out Reset image size

The value of the cusp relation for every mock lens system is dependent on the relative positions of the subhalos and the source. The placement of the subhalos and the source is stochastic—depending upon the merger history of the lens halo, the direction in which the halo is projected, and the sampling of the source plane. Thus, the results of Figure 8 could be due to just a few substructure realizations that always give high-R_cusp values, while the majority of substructure realizations give small R_cusp values regardless of f_sub. We address this possibility in Figure 9. For every bin of f_sub, we measure the fraction of substructure realizations that have at least one large cusp relation value among their set of mock lenses. A small fraction of the realizations with small values of f_sub have large R_cusp values, while a large fraction of the realizations with large values of f_sub have large R_cusp values. This behavior is most pronounced for a f_sub definition that includes only subhalos whose center falls near the image positions, even when considering that, using this definition, the bins with f_sub > 0.05 contain only a few realizations each. Thus, f_sub is reliably correlated with R_cusp values.

Zoom In Zoom Out Reset image size

The parameterization of substructure in lensing halos which is most often used in studies of anomalous flux ratios is f_sub. Other possibilities for parameterizing substructure exist: the number of subhalos and the average subhalo mass.⁵ We test these possibilities in Figure 10. Here, subhalos whose halo centers fall within R = 2b are counted and the mean subhalo mass measured. The substructure realizations with a larger number of subhalos also have larger and more frequent violations of the cusp relation, when compared to realizations with a smaller number of subhalos. In addition, the substructure realizations with a larger average subhalo mass have larger and more frequent violations of the cusp relation, when compared to realizations with a smaller average subhalo mass. However, when compared to Figure 8, neither the number of subhalos nor the mean subhalo mass is as discriminating as f_sub. For example, for the observed value of R_cusp = 0.25 (for opening angle bin 70°–80°), Figure 10 suggests that the mean subhalo mass and the number of subhalos could be large or small. On the other hand, Figure 8 is clear in showing that f_sub must be large. This result is consistent with that of Metcalf & Amara (2010), who find that the frequency of flux anomalies is a function of f_sub and not of the number or average mass separately.

Zoom In Zoom Out Reset image size

In this section, we discuss some specific caveats to our calculation that we are aware of and should be taken into account in the interpretation of these results.

In Section 4, we discussed how the simulations of substructure differ between the semi-analytic host halo and the typical lensing halo and the accommodations we make in order to incorporate the substructure catalogs into our mock lenses. We make simple adjustments—adopting the projected distance distribution and projected mass fraction from the simulations while scaling the mass of the host galaxy from M_h = 1.8 × 10¹² M_☉ to M_lens = 10¹³ M_☉. The scaleability of the halo and subhalo masses is in line with the results of DM numerical simulations (Gao et al. 2004; Giocoli et al. 2010). The host halo in the semi-analytic realizations also differs with respect to typical lensing halos in more subtle ways: the formation and merger history and the epoch of observation (Madau et al. 2008). While these differences are not likely to significantly effect the distance distribution or mass fraction of substructure (e.g., Xu et al. 2009, and as discussed in Section 4), they have unquantified effects on the internal structure of the subhalos—their scale and tidal radii. The concentration of a subhalo, however, is a slowly varying function of mass and time, so that biases in this analysis induced by these effects are likely to be small compared to natural halo-to-halo variation.

With regards to finite source effects, in this work we use the point-source approximation, which—as stated previously—could bias the calculated magnification of high magnification images. For example, Metcalf & Amara (2010) suggest that finite source effects could have large effects on the magnifications of images. For lens systems with substructure and total magnifications above μ = 10, the fractional deviation in magnification for a point source and for a 1 pc source approaches 1. For smooth lens halos, the same holds true for magnifications above μ = 20. The bias in magnifications most often manifests itself as an overestimate of the magnification of an image when a point-like source is used in place of a finite source.

However, tests using our code fail to find significant biases due to the point-like source assumption. This difference can be accounted for by two effects. First of all, we consider only the cusp relation, which is the difference between magnifications of images, mitigating the effects of any finite source effect. Second, Metcalf & Amara (2010) use a model of substructure where the surface density profile is a power law Σ∝R^−α, and α = 1 or α = 0.5. This profile is significantly steeper than the NFW profile used in this work and is likely to cause bigger magnification perturbations than the subhalos in our work.

Finally, with regards to magnification bias, the distribution of opening angle values in our mock lenses, as shown in Figures 5 and 6, is not representative of the distribution expected in observed lens samples. For example, our mock data set has fewer lenses with small opening angles than expected observationally. This is because the distribution of opening angles in our mock samples is a result of a combination of the macromodel parameters and the manner in which the source plane is sampled in order to create the mock lens systems. We sample the source plane uniformly, a method which does not account for the magnification bias of observations—the fact that brighter systems are more likely to be observed. Following Keeton & Zabludoff (2004), it can be shown that the probability density in the source plane is related to the luminosity function of sources, such that

$$\begin{equation} p({\bf u})d{\bf u} = \mu ({\bf u})^{\nu -1} d{\bf u}, \end{equation} \tag{ 8 }$$

when the luminosity function of sources has a power-law form dN/dS∝S^−ν, where S is the source flux and where $\bf {u}$ is the position of the source. Thus if ν = 2, the probability of observing a lens system with source position $\bf {u}$ is proportional to the total magnification of the system, μ. Studies of the Cosmic Lens All-Sky Survey (Browne et al. 2003) have suggested that ν = 2.1 (Rusin & Tegmark 2001; Chae 2003). The total magnification of a system is related to the opening angle of a cusp so that the smallest opening angles are found in the systems with the largest magnifications. In our observational sample, one of ten in the observational set has an opening angle less than 40°. The mock lens systems, in contrast, have very few systems (<1%) with opening angles less than 40°. Accounting for magnification bias would increase the number of systems with opening angles less than 40° to a few percent. Given the effect of magnification bias, we negate its effect by making comparisons of observed lens systems to our mock lens systems using only data that are binned in opening angles.

DM substructure is a generic prediction of the CDM model and probes of substructure are a key test of the nature of the DM. Strong gravitational lensing can be used to detect and measure DM subhalos, as they cause perturbations to the image magnifications, positions, and time delays in lens systems with galaxy lenses. Putting robust constraints on these observational results require robust predictions for CDM substructure that account for halo-to-halo variation. Current predictions from numerical simulations have used a limited number of simulated halos. We investigate the question of halo-to-halo variation in substructure using semi-analytical models and test how it affects perturbations to the magnifications of images, as measured by the cusp relation parameter, R_cusp. We find the following results.

• 1.  
  The total projected substructure mass fraction for halos with mass greater than 10⁻⁵ of the mass of the host is ∼9% with a 95 percentile in the range [2–17]%. Low-concentration halos have a higher projected mass fraction (∼20%), with a significantly larger 95 percentile range. High-concentration halos have a lower mass fraction (∼5%) and a much smaller 95 percentile range.
• 2.  
  We find that the mass fraction within the inner 3% of the projected radius is f_sub = 0.25%, with a 95 percentile in the range of f_sub = [0–1]% for the complete sample of 200 realizations. For low-concentration host halos, the median substructure mass fraction is f_sub ∼ 0.5% with a 68 percentile in the range of f_sub = 0.2%–0.8%. These results are generally consistent with the results of Xu et al. (2009), who find 0.01% ⩽ f_sub ⩽ 0.7%, using six halos from the Aquarius simulation, as well as with all previously published limits on f_sub. The mass fraction of host halos with higher concentrations is lower, with f_sub ⩽ 0.1% for host halos with concentrations in the upper 16% of the distribution.
• 3.  
  Accounting for halo-to-halo variation, the observed values of R_cusp are greater than the values predicted by smooth lens models and larger than 99% of the values predicted by substructure models. Given halo-to-halo variation, the probability of observing cusp relation values greater than supported by smooth models is ∼10⁻² and the probability of observing cusp relation values greater than the largest currently observed values is ∼10⁻³.
• 4.  
  Substructure fractions are correlated with the concentrations of the host halos. Lower-concentration host halos have larger substructure fractions so that lower concentration host halos give rise to more substructure and more violations of the cusp relation.
• 5.  
  For lensing studies of substructure, the substructure content of lensing halos is best parameterized by a definition of the substructure mass fraction, f_sub, which includes only the subhalos whose halo centers fall within the area in which images are expected to be found (twice the Einstein radius in our models, ∼3% of the virial radius of the lens halo, generically) and which excludes more distant subhalos that fall within a scale or tidal radius of a lensed image.
• 6.  
  f_sub is a better predictor of R_cusp when compared to the number of subhalos or the mean subhalo mass.

We attempt to isolate the effect of halo-to-halo variation in substructure populations from macromodel effects and, thus, do not model any particular observed lens system. Our efforts can still be compared qualitatively to previous attempts to constrain the substructure mass fraction. For example, Dalal & Kochanek (2002) find that f_sub = 0.02 with 0.6% < f_sub < 7% within 90% confidence intervals. This result is marginally consistent with predictions of f_sub in previous works and in this paper, while Metcalf & Amara (2010) find rough consistency with such predictions. Both of these works use substructure models with density profiles that are steeper than the NFW profile used in this work. Steeper profiles are more effective at causing magnification perturbations. Thus, it is possible that both of these works are underestimating the f_sub that is supported by lensing data. Other works have also suggested that lensing observations are inconsistent with expectations from numerical simulations (Macciò et al. 2006; Amara et al. 2006; Chen 2009). Xu et al. (2009) suggest that, using the six halos in the Aquarius simulation, the probability of the observed values of R_cusp is ∼10⁻³.

In summary, we find that when accounting for halo-to-halo variation the probability of observing R_cusp values that match observations may be as large as ∼10⁻², but the tension between observations and numerical simulations persist. This tension may be eased by assuming larger values of the substructure mass fraction. Possible solutions, then, are (1) additional substructures which are unresolved in current numerical and semi-analytic simulations and (2) biases which cause the observational sample to prefer lensing halos with high substructure mass fraction or, relatedly, low concentration. Both these possibilities, however, are speculative.

While this work only probes the effect of subhalos with masses ∼10⁸ M_☉ or larger, others have studied smaller substructures. Using a ΛCDM model, Macciò & Miranda (2006) investigated the role of subhalos with masses 10⁵–10⁷ M_☉ and found them to be insufficient in creating large cusp relation violations. The Aquarius simulation used by Xu et al. (2009) resolves subhalos down to 10⁵ h⁻¹ M_☉. The cumulative effect of smaller and larger subhalos produces their result of a probability that observations match expectations from ΛCDM of ≈10⁻³.

Currently, there are no known observational biases which would prefer high substructure or low-concentration lensing halos. If the observational bias from the concentration of lensing halos is considered alone, high-concentration halos should be preferred (Mandelbaum et al. 2009). In addition, the relation between observational biases, substructure content, and other characteristics of lenses may be complicated. For example, triaxial halos projected along the long axis are better strong lenses in general and may have more substructure projected near images (Zentner 2006), but triaxial halos along the middle axis are better at creating four-image lenses (Rozo et al. 2007). Accounting for all the biases due to lensing halo characteristics is difficult and requires future investigation.

S.M.K. is funded by NSF grant PHY-0969853 and by Brown University. A.R.Z. is funded by the National Science Foundation under grant PHY 0968888.