Numerical investigation of lens models with substructures using the perturbative method.
Abstract
We present a statistical study of the effects induced by substructures on the deflection potential of dark matter halos in the strong lensing regime. This investigation is based on the pertubative solution around the Einstein radius (Alard 2007) in which all the information on the deflection potential is specified by only a pair of onedimensional functions on this ring.
Using direct comparison with raytracing solutions, we found that the isocontours of lensed images predicted by the pertubative solution is reproduced with a mean error on their radial extension of less than — in units of the Einstein radius, for reasonable substructure masses. It demonstrates the efficiency of the approximation to track possible signatures of substructures.
We have evaluated these two fields and studied their properties for different lens configurations modelled either through massive dark matter halos from a cosmological Nbody simulation, or via toy models of Monte Carlo distribution of substructures embedded in a triaxial Hernquist potential.
As expected, the angular power spectra of these two fields tend to have larger values for larger harmonic numbers when substructures are accounted for and they can be approximated by powerlaws, whose values are fitted as a function of the profile and the distribution of the substructures.
keywords:
methods: Gravitational lensingstrong lensing; Nbody simulations1 Introduction
The cold dark matter (CDM) paradigm (Cole et al. 2005 and references therein) has led to a successful explanation of the largescale structure in the galaxy distribution on scales 0.02 0.15h Mpc. The CDM power spectrum on these scales derived from large redshift surveys such as, for instance, the AngloAustralian 2degree Field Galaxy Redshift Survey (2dFGRS), is also consistent with the Lyman forest data in the redshift range (Croft et al. 2002; Viel et al. 2003; Viel, Haehnelt & Springel 2004).
In spite of these impressive successes, there are still discrepancies between simulations and observations on scales 1 Mpc, extensively discussed in the recent literature. We may mention the sharp central density cusp predicted by simulations in dark matter halos and confirmed by the rotation curves of low surface brightness galaxies (de Blok et al. 2001) or in bright spiral galaxies (Palunas & Williams 2000; Salucci & Burkert 2000; Gentile et al. 2004). Moreover, deep surveys (), such as the Las Campanas Infrared Survey, HST Deep Field North and Gemini Deep Deep Survey (GDDS) are revealing an excess of massive earlytype galaxies undergoing “topdown” assembly with high inferred specific star formation rates relative to predictions of the hierarchical scenario (Glazebrook et al. 2004; Cimatti, Daddi & Renzini 2006).
One problem that requires closer examination concerns the large number of sub subhalos present in simulations but not observed (Kauffmann, White & Guiderdoni 1993; Moore et al. 1999; Klypin et al. 1999). This is the case of our Galaxy or M31, although there is mounting evidence for a large number of very low mass dwarfs (Belokurov et al. 2006). However, it is still unclear whether the CDM model needs to be modified to include selfinteracting (Spergel & Steinhardt 2000) or warm dark matter (Bode, Ostriker & Turok 2001; Colín, AvilaReese & Valenzuela 2000) or whether new physical mechanisms can dispel such discrepancies with the observations. For instance, gas cooling can be partly prevented by photoionization process which may inhibit star formation in the majority of subhalos (Bullock, Kravtsov & Weinberg 2001).
This “missing satellite problem” remains an ideal framework to test cosmological models. During the past years, different methods have been employed in order to study the gravitational potential of groups or clusters of galaxies, for instance through their Xray lines emission of hot gas in the intracluster medium or through lensing considerations. However, while lensing directly probes the mass distribution in those objects, the other methods rely more often than not on strong hypotheses on the dynamical state of the gas and interactions between baryons and dark matter. For example, the gas is supposed to be in hydrostatical equilibrium in the gravitational potential well created by dark matter halo, while spherical symmetry is assumed. In this paper, we study the effects induced by substructures on the deflection potential of dark matter halos in the strong lensing regime. The presence of substructures follows from the capture of small satellites which have not yet been disrupted by tidal forces and/or suggests that the relaxation of halos is not totally finished.
Raytracing through Nbody cosmological simulations suggest that substructures should have a significant impact on the formation of giant arcs. At clusters of galaxies scales, some results indicate that lensing optical depths can be enhanced (Bartelmann, Steinmetz & Weiss 1995; Fedeli et al. 2006; Horesh et al. 2005; Meneghetti et al. 2007a) whereas some other recent studies suggest that the impact on arc occurence frequency should only be mild (Hennawi et al. 2007). On the other hand, the presence of substructures changes the properties of strongly lensed images to a point that could lead to misleading inferred cluster mass properties if not properly accounted for (Meneghetti et al. 2007b). Assuming a onetoone association of cluster subhalos in the mass range and galaxies, Natarajan, De Lucia & Springel (2007) used weak lensing techniques and found the fraction of mass in such subhalos to account for of the total cluster mass, thus in good agreement with predictions from simulations (Moore et al. 1999).
Likewise, at the scales of galaxies, strong lensing events have long been involving multiple quasars for which the impossibility of resolving the size or shape of the lensed images brought the attention toward flux ratios of conjugate images as a probe of subtructures. Depart from flux ratios expectations from a simple elliptically symmetric potential is often interpreted as a signature for local potential perturbations by substructures (Bradač et al. 2002; Dalal & Kochanek 2002; Bradač et al. 2004; Kochanek & Dalal 2004; Amara et al. 2006). It is unclear whether anomalous flux ratios actually probe “missing satellites” (Keeton, Gaudi & Petters 2003; Mao et al. 2004; Macciò et al. 2006). Due to the small source size in the case of lensed QSOs, sensitivity to microlensing events due to stars in the lens galaxy makes the interpretation less obvious. Astrometric perturbations of multiple quasars have also been considered (Chen et al. 2007) despite substantial observational limitations.
Presumably the best way out would be to consider extended sources like QSOs observed in VLBI or lensed galaxies that will be sensitive to a narrower range of scales for the pertubing potential and thus easier to interprete. New methods for inverting potential corrections that needed on top of a smooth distribution were proposed (Koopmans 2005; Suyu & Blandford 2006) but are not guaranteed to converge in all practical cases and seem to depend on the starting smooth distribution.
One interesting alternative approach is to treat all deviations from a circularly symmetrical potential as small perturbations (Alard 2007, 2008) defining the location where multiple extended images will form. Two perturbative fields, and , can then be defined to characterize deflection potential of lenses as a function of the azimuthal angle , near the Einstein radius. They respectively represent the radial and azimuthal derivative of the perturbated potential (see Eq. (10) below). Alard (2008) showed that these two fields have specific properties when one substructure of mass of the total mass is positioned near the critical lines. For instance, the ratio of their angular power spectra at harmonic number is nearly . We will investigate the detailed properties of these perturbative fields by considering more realistic lenses such as dark matter halos extracted from cosmological simulations. In order to control all the free parameters (mass fraction, and shapes of subtructures for instance) and to study their relative impact on arc formation, we will also generate different families of toy halos.
This paper is organized as follows: in section 2 we present our lensing modelling; section 3 first sketches the pertubative lens solution and applies it to our simulated lenses for validation against a ray tracing algorithm; section 4 presents our main results on the statistics of perturbations, while the last section wraps up.
2 Numerical modelling
2.1 Lens model
Halos formed in cosmological simulations tend to be centrally cuspy () and are generally not spherical, but have an triaxial shape. The triaxiality of these potential lenses is expected to increase significantly the number of arcs relative to spherical models (Oguri, Lee & Suto 2003 and references therein), and must be taken into account in numerical models. Thus, apart from dark matter halos extracted from cosmological simulations, we consider in this work typical lenses modelled by of a dark matter halo of total mass with a generalized Hernquist density profile (Hernquist 1990):
(1) 
where is the value of the scale radius, a triaxial radius defined by
(2) 
and and the minor:major and intermediate:major axis ratio respectively. We decided to use an Hernquist profile for practical reasons. However, for direct comparison with common descriptions of halos from cosmological simulation in the literature, the Hernquist profile is related to an NFW profile (Navarro, Frenk & White, 1996; 1997) with the same dark matter mass within the virial radius ^{1}^{1}1 defines the sphere within which the mean density is equal to 200 times the critical density.. Moreover, we also impose that the two profiles are identical in the inner part () which can be achieved by using relation (2) between and the NFW scale radius in Springel, Di Matteo & Hernquist (2005). By convention, we use the concentration parameter in the following to characterize the density profile of our lenses. For example, for typical lens at a redshift , we use kpc which corresponds to a NFW profile with (or equivalently kpc, kpc) and is consistent with values found in previous cosmological Nbody simulations at the specific redshift and in the framework a the CDM cosmology (Bullock et al. 2001; Dolag, Bartelmann & Perrotta 2004).
Axis ratios of each lens are randomly determined following Shaw et al. 2006: , and . These values are in good agreement with previous findings from cosmological simulations (see for instance Warren et al. 1992; Cole & Lacey 1996, Kasun & Evrad 2005). Finally, it is worth mentioning that each halo is made of particles corresponding to a mass resolution of . However, we impose a troncation at a radius of value Mpc.
2.2 Substructures model
2.2.1 Mass function
Numerical Nbody simulations show that dark matter halos contain a large number of selfbound substructures, which correspond to about 1020% of their total mass (Moore et al. 1999). In the following, the number of substructures N in the mass range – + is assumed to obey (Moore et al. 1999; Stoehr et al. 2003)
(3) 
The normalization constant is calculated by requiring the total mass in the clumps to be of the halo mass and by assuming subhalos masses in the range – M. The minimum number of particles in the substructures is about , while the more massive ones have particles.
2.2.2 Radial distribution
Substructures are distributed according to the (normalized) probability distribution , where is assumed to have an Hernquist profile of concentration , which yields the probability to find a clump at a distance within the volume element . While the abundance of subtructures in halos of different masses has recently been extensively quantified in cosmological simulation (see for instance Vale & Ostriker 2004; Kravtsov et al. 2004; van den Bosh et al. 2007), their radial distribution is less understood. However, some studies seem to suggest their radial distribution is significantly less concentrated than that of the host halo (Ghigna et al. 1998, 2000; Colín et al. 1999; Springel et al. 2001; De Lucia et al. 2004; Gao et al. 2004; Nagai & Kravtsov 2005; Macciò et al. 2006). We will use either in good agreement with those past investigations, or for comparison.
2.2.3 Density profiles and alignment
The stripping process caused by tidal forces seems to reduce the density of a clump at all radii and, in particular, in the central regions, producing a density profile with a central core (Hayashi et al. 2003). This process was further confirmed by simulations which found that the inner structure of subhalos are better described by density profiles shallower than NFW (Stoehr et al. 2003). However, other simulations seem to indicate that the central regions of clumps are wellrepresented by power law density profiles, which remain unmodified even after important tidal stripping (Kazantzidis et al. 2004a). This effect may be enhanced when star formation is taken into account since dissipation of the gas (from radiative cooling process) and subsequent star formation lead to a steeper dark matter density profile due to adiabatic contraction. To test the importance of these differences from the point of view of arcs formation, we allow for both of these possibilities: we simulate halos with clumps having a central core , where defines a core radius, or a central cusp (Hernquist profile) and study how our resulting arcs would change from one option to the other. It is worth mentioning that each subhalo concentration parameter is obtained using relation (13) in Dolag et al. (2004) within the CDM cosmology. However, to avoid spurious effects due to the lack of resolution, all subhalos represented by less than 200 particles will have a concentration parameter value corresponding to that of an halo made of exactly 200 particles (i.e ). For core profiles, we follow Hayashi et al. (2003) and take a core of size .
Finally, recent cosmological simulations suggest that subhalos tend to be more spherical than their host (Pereira et al. 2008; Knebe et al. 2008) and this effect can also be enhanced if halos are formed in simulations with gas cooling (Kazantzidis et al. 2004b). Moreover, the distribution of the major axes of substructures seems to be anisotropic, the majority of which pointing towards the center of mass of the host (Aubert, Pichon & Colombi 2004; Pereira et al. 2008). Although shapes and orientations of subhalos provide important constraints on structure formation and evolution, we have not studied their relative influence in this work. We reasonably think that modelling substructures by either triaxial shapes of spherical shapes won’t lead to any significant differences in our results.
2.3 Lens samples
Table 1 summarizes our different samples of lenses. From a statistical point of view, each sample involves one hundred realizations of halos following the above described methodology. Lenses are assumed to be at a typical redshift . They share a common total mass of and are thus described by an Hernquist profile of concentration . The sample A represents our reference catalogue in which all halos have no substructure. For each of them, we introduce a fraction of subtructures () by removing some background particles so that both the total mass and the density profile of the initial halo are conserved. These halos are classified in catalogues B and C according the definition of inner density profile (IP) of clumps. For example, each halo from samples B1 and B2 have 15% of substructures with an inner profile represented by a cusp (Hernquist profile). Their radial distribution (RD) within the halo is left as a free parameter. We consider two possibilities, as suggested by numerical simulations, and , which is the concentration parameter of the whole halo. Finally, lenses catalogues C1 and C2 have substructures represented by a core profile with and respectively.
In addition, lensing efficiency depends on the relative distance between lenses and sources. The efficiency of the lens is scaled by the critical density
(4) 
where , and are the angular diameter distances between the observer and the source, between the observer and deflecting lens and between the deflector and the source respectively. When the surface mass density in the lens exceeds the critical value, multiple imaging occurs. In order to account this effect which, for a given lensing halo, implies a different Einstein radius for a different source redshift, we consider the redshift distribution of sources taken from the COSMOS sample of faint galaxies detected in the ACS/F814W band (Leauthaud et al. 2007). It is well represented by the following expression
(5) 
with and (Gavazzi et al. 2007).
Sample  Sub F  Sub IP  Sub RD 

0%      
15%  cusp  
15%  cusp  
15%  core  
15%  core 
3 numerical validation of the perturbative solution
In this section, we take advantage of the large sample of mock lenses described in § 2 to assess the validity of the perturbative method developed in (Alard 2007). After a brief presentation of the basic idea in § 3.1, we compare the ability of this simplified procedure to reproduce multiple images lensed by complex potentials as compared to a direct raytracing method (§ 3.3) and define the validity range of the perturbative approach.
3.1 The perturbative approach
For the sake of coherence, let us sketch the motivation behind the perturbative lens method (Alard 2007) used througout this paper. The general lens equation, relating the position of an image on the lens plane to that of the source on the source plane can be written in polar coordinates as
(6) 
where is the source position, and r, and are the radial distance, radial direction and orthoradial direction respectively. Here is the projected potential. Let us now consider a lens with a projected density, , presenting circular symmetry, centered at the origin, and dense enough to reach critical density at the Einstein radius, . Under these assumptions, the image by the lens of a point source placed at the origin is a perfect ring, and equation (6) becomes:
(7) 
where the potential, , is a function of r only, and the zero subscript refers to the unperturbed solution. The basics ideas of the perturbative approach is to expand equation (7) by introducing i) small displacements of the source from the origin and ii) noncircular perturbation of the potential, which can be described by:
(8) 
where is small number: . To obtain image positions (,) by solving equation (6) directly, may prove to be analyically impossible in the general case. It is then easier to find perturbative solution by inserting equation (8) into equation (6). For convenience, we rescale the coordinate system so that the Einstein radius is equal to unity. The response to the perturbation on may then be written as
(9) 
which defines , the azimuthally dependant enveloppe of the relative deflection. Using Equation (8), the Taylor expansion of is
(10) 
where:
(11) 
3.2 Morphological effects VS astrometric distortions
As demonstrated in Alard (2008), the morphology of arcs is very sensitive to small perturbators such as substructures in the main halo. The effect of substructures on morphological features like for instance the size of an image is typically much larger than pure astrometric distortions. Indeed, astrometric distortions are only of the order of the substructure field, which in many case is smaller than the PSF size of the instrument used, and represents then an observational challenge to measure. Note that a perturbative theory of astrometric distorsions was already considered by Kochanek et al. (2001) and Yoo et al. (2005; 2006). Note also that the former perturbative approach is limited to astrometric effects, the relevant theory does not describe image formation, and thus cannot predict morphological effects. When image morphology is considered, effects are an order of magnitude larger than astrometric effects. For instance, let us consider the examples displayed in Fig.1 and Fig. 2 of Alard (2008). Whereas a giant arc is obtained from an unperturbed elliptical lens cusp caustic (Fig. 1), the introduction near the Einstein radius of a substructure of only 1 % of the main halo mass breaks the arcs into 3 subimages (Fig. 2). The detection of such effects should not require a particularly good resolution, since the amplitude of the effect is a fraction of the arc size, which is typically several time the PSF size.
Specifically, let us consider the following example: we take a subcritical configuration for an elliptical lens and evaluate the modification of the image size due to the perturbation by a substructure field. Let us assume that the ellipticity of the lens is aligned with the axis system and both source and the substructure are placed on the Xaxis. In this configuration, the size of the central image will be perturbed by the substructure field, and this will be the observable effect. The main interest of this example is the simplicity of the calculations (linear local description of the field) and the simple description of the effect (reduction of image size). Let us now define some useful quantities: the local slope of the unperturbed field , the angular size of the unperturbed image and the source size . The substructure parameters are its mass (in unit of the main halo mass) and its distance to the Einstein ring. For circular sources, the local slope, image size and source radius are related by the following relation in the unperturbed case:
(14) 
The slope perturbation introduced by the substructure is (See Alard 2008, Sec. 2):
(15) 
The substructure field modifies the images size according to:
(16) 
which gives:
(17) 
For an elliptical isothermal potential with ellipticity parameter , the field reads (see Alard 2007, 2008):
(18) 
Thus, which is the field derivative in 0 is: leading to:
(19) 
The former equation evaluates the angular perturbation of the image size by the substructure. To obtain the corresponding perturbation on the image length , we have to multiply by the Einstein Radius:
(20) 
Note that astrometric effects are identical to the effects of the substructure on the field for circular sources. Thus, the astrometric effect is of the order (See Alard 2008, Sec. 2):
(21) 
Consequently, the ratio between morphological and astrometric effects here is:
(22) 
For arcs, the source size (source diameter= ), and the parameter have the same typical scale, which gives: , and ; then:
(23) 
This means that morphological effects are 10 times larger than astrometric effects. This point is critical, since it really makes the effect of substructure observable. To illustrate this latter point, let’s consider some numerical values for typical galaxies. For a Milky way like galaxy, we have a mass . The Einstein radius in arcsec is given by:
(24) 
which gives the size of the image perturbation, :
(25) 
For a perturbator with 0.5 % mass of the galaxy, taking typical scales, , , and , we obtain:
(26) 
Such effects should be within the reach of spatial instruments such as DUNE or SNAP. On the contrary, astrometric effects are typically 10 times smaller and should be then much more difficult to detect. Consequently morphological effects are definitely our best hopes to detect substructures.
Clearly, other configurations will also allow the measurement of morphological effect of substructures, though in general the calculations will be a little bit more complicated, but the effects will more or less be of the same order, since they are related to the amplitude of the perturbation of the field by the substructure.
3.3 Reconstruction of images
One interesting feature of the perturbative method is to provide a framework for the reconstructions of images. By first defining an elliptical source centered on position , with a characteristic size , ellipticity , and inclination of the main axis such that,
(28)  
(30)  
(31) 
one can express the equations of the image contours using Eq. (3.1): (see Alard (2007) for more details)
(33)  
(35)  
(36) 
This equation corresponds to Eq. (15) in Alard (2007). The functional is defined to take into account the effect of the translation of the source by the vector .
(37) 
As emphazised in Alard (2007), the image contours are only governed by the two fields, and , which contain all the information on the deflection potential at this order in the perturbation. For instance, the two first terms in the bracket of Eq. (33) give informations on the mean position of the two contour lines while the last term provides informations on the image’s width along the radial direction as well as a condition for image formation. Therefore, the characterization of these two fields represents a simple and efficient way to track possible signatures of the deflection potential induced by substructures, so long as the perturbative framework holds, as we will illustrate below. Alard (2007) already implemented the method with a lens described by a NFW profile yielding an analytical solution for the projected potential profile. In this section, we illustrate and validate the method while considering more complicated and realistic situations. In particular, we use lenses either from cosmological simulations or from toy models presented in section 2.
For direct comparison between arc reconstructions predicted by the perturbative method and theoritical ones, we use a raytracing method. Part of our investigations indeed makes use of the Smooth Particle Lensing technique (SPL), described in details in Aubert, Amara & Metcalf (2007) and summarized in this section. SPL has been developed to compute the gravitational lensing signal produced by an arbitrary distribution of particles, such as the ones provided by numerical simulations. It describes particles as individual light deflectors where their surface density is arbitrarily chosen to be 2D Gaussian. This choice makes it possible to compute the analytical corresponding deflection potential, given by:
(38) 
where , is the mass of the particle, its extent and the critical density. From the deflection potential, expressions for the deflection angles , the shear components and the convergence can be easily recovered (see Aubert, Amara & Metcalf (2007) for more details). Knowing the lensing properties of a single particle, one can recover the full signal at a given ray’s position on the sky by summing the contributions of all the individual deflectors:
(39)  
(40) 
where, e.g. is the contribution of the ith particle to the shear at ray’s position . This summations are performed efficiently by means of 2DTree based algorithm, in the spirit of Nbody calculations. The tree calculations are restricted to monopolar approximations where an opening angle of is found to give results accurate at the percent level on analytical models. Finally, Aubert, Amara & Metcalf (2007) found that an adaptative resolution (i.e. an adaptative extent for the particles) provides a significant improvement in the calculations in terms of accuracy. For this reason, the smoothing depends on the rays location: particles shrink in high density regions in order to increase the resolution while they expand in lowdensity regions, smoothing the signal in undersampled areas. For all simulations, we use rays within a square of size , an opening angle of 0.7 and (where is the number of particles over which the smoothing is applied).
3.3.1 Lenses from the toy model
We present in this section characteristic examples of arc reconstruction. Three lenses , and belonging to samples A, B2 and C2 respectively are considered. They have a common mass, density profile, axis ratios and random orientation in 3D space. They only differ via the presence or not of substructures as well as via the inner density profile of substructures: has no substructure whereas has substructures with a central cusp while a core describe the inner density profiles of substructures in . In the present case, the source is at a redshift , has an elliptical contour with and a radius , characterizing a Gaussian luminosity profile (see figure 1).
For the arc reconstructions presented below, we shall consider 3 different radii , and corresponding to 3 specific isophotal contours defined by , and (see figure 1).
In Figure (2), we show the projected mass density of lenses and near the Einstein radius, the image’s solution obtained from raytracing and the contours predicted by the perturbative method when the source is placed at the origin. When no substructure is considered, both projected density and potential are nearly elliptical. As expected, we obtain four distinct arcs in a cross configuration. The predicted arcs reconstruction are in good agreement with the numerical solution obtained via the raytracing algorithm. From a theoretical point of view, it is easy to show that the functions and are proportional to and respectively. These functionnal form are recovered in our experiment and shown in Figure (3).
When substructures are present, the shape of images is significantly altered. First, we notice that the positions of substructure tend to break the ellipticity of the halo center. Thus, it is not surprising that the shape of the image is approaching a ring in that case. Moreover, it is interesting to see that the position of one substructure (at the top left of the figure) is exactly at the Einstein radius. This produces an alteration of the luminosity, while the effects is more violent when substructures present a cups profile. This effect can be clearly seen when comparing the perturbative fields and relative to the lens in Figure (3). For instance, we can see two clear bumps in the evolution of . The second one () is produced by substructures in the lower right part and induced an alteration of the luminosity again.
To estimate the systematic error between the theoretical contours provided by the raytracing and those predicted from equation (33), we use a simple procedure with a low computational cost. First, each predicted contour is divided into a sample of N points. Each of them is defined by polar coordinates (, ) which coincide with a luminosity value of the image corresponding to an unique radius in the source frame. By using relation (33), we then compute which gives the image contour radius of the isophotal contour of the source in Einstein radius unit. By defining the radial distance of point i, the mean error (in Einstein radius unit) is then computed by .
For illustration, we have estimated the mean error reached for the lenses and using three luminosity contours source (see fig. 1). For , the mean errors are respectively 0.67%, 0.71%, 0.95% of the Einstein radius for isophotes , and respectively, while we obtain 0.74%, 0.86% and 1.04% for lens and for same luminosities. We have also studied how the mean error evolves for random positions of the source inside an area limited by the caustic lines. To do that, we have used the lens and have studied 100 realizations with different impact parameters. We found for isophotes equal to .
3.3.2 Lenses from cosmological simulation halos
In this section, lenses are modelled by dark matter halos extracted from a cosmological simulation of the Projet HORIZON^{2}^{2}2http://www.projethorizon.fr/. The simulation was run with Gadget2 (Springel 2005) for a CDM universe with , , , km/s/Mpc, in a periodic box of 20 Mpc. We use particles corresponding to a mass resolution of and a spatial resolution of 2 kpc (physical). Initial conditions has been generated from the MPgrafic code (Prunet et al. 2008), a parallel (MPI) version of Grafic (Bertschinger 2001). In this simulation, we selected two regions. In the first one, the lens is a typical halo of total mass at redshift . The source is at , assumed to be elliptical () with and placed near a caustic in order to obtain a giant arc. Fig. (4) shows the projected density of the lens near the Einstein radius; both the raytracing solution and the predicted contours by the perturbative method are shown. Here again, the three different contours are well reconstructed since the error are 0.76%, 0.83% and 0.91% for isophotes , and respectively. For illustration, the angular variation of the perturbative fields are also representated in the figure (4).
The second example is a lens modelled by another halo from the same Nbody simulation. Its total mass is at . This is an extreme case since a significant number of substructures are still falling toward the center of the host halo which suggests that the dynamical relaxation is still operating. This strongly affects the potential and the perturbative fields (see figure 5). However, mean errors remain small of the order of 1.11%, 1.20% and 1.26% for isophotes , and respectively, which proves the accuracy of the method to deal with more complex systems. We may reasonably think that lenses in our different samples at tend to be more relaxed that the present configuration, and, consequently, the mean errors in our arc reconstruction should be less pronounced.
4 Fourier series expansion
4.1 Motivation
Given Eq. (33), it is straightforward, for a given lensed image and an assumed underlying spherical lense and elliptical source, to invert it for and as
(41)  
(42) 
where , and are given by Eq. (33). This inversion formula depends explicitly on the source parameters, which are unknown. However using Eqs ( 41, 42) it is possible to compute the two functions and for each couple of parameters . The proper solution corresponding to the true parameter has minimal properties. Consider for instance a circular solution (), if the inversion formula is used with additional Fourier terms with order will appear in the inversion formulae. Thus, it is clear that minimizing the power in higher order Fourier modes is a criteria that will allow to select the best solution when exploring the plane . This criteria has also a very interesting property, considering that power at order usually reveal the presence of substructures (see Table 3 for instance), the solution with minimum power at higher order is also the one that puts the more robust constraint on the presence of substructure. Thus the elliptical inversion can be performed by exploring the plane in a given parameter range, computing the corresponding fields and and their Fourier expansion, and selecting the solution with minimum power at . For non elliptical sources, one can use the general inversion method presented in Alard (2008). This inversion method remaps the images to the source plane using local fields models (basically the scale of the images). The solution is selected by requiring maximum similaritiy of the images in the source plane. Image similarity is evaluated by comparing the image moments up to order . Provided the number of image moments equations exceed the number of model parameters, the system is closed and has a definite solution. Note that the local models may be replaced with general Fourier expansion in the interval , but in this case, the additional constraint that no image are formed in dark areas must be implemented (Diego et al. 2005). In the perturbative approach this requirement can be reduced to in dark areas, where is the radius of the smallest circular contour that contains the source. We may therefore assume for now that observational data may be inverted, and that an observationnal survey of arcs should provide us with a statistical distribution of the perturbatives fields. Hence we may use our different samples of halos presented in paragraph 2.3 in order to measure the relative influence on arc formation of the different free parameters such as the inner profile of substructures or their radial distribution within the host halo. To conduct this general analysis the fields will be represented by Fourier models, due to the direct correspondance between Fourier models of the fields and the multipolar expansion of the potential at (Alard 2008).
4.2 Results
The angular functions and can be characterized by their Fourier expansion:
(43)  
(44)  
(45) 
where , correspond to associated power spectra. We have derived the multipole expansion of and for each halo of the different catalogues and we focus in the following on the mean amplitudes and obtained.
Tables (2) and (3) respectively summarise the seven first orders of the power spectrum of and for the 3 lenses , and .
Lens  1  2  3  4  5  6  7 

0.07  4.21  0.02  0.20  0.04  0.07  0.03  
1.62  3.80  0.42  0.18  0.29  0.20  0.33  
1.38  2.86  0.18  0.20  0.10  0.11  0.11 
Lens  1  2  3  4  5  6  7 

0.08  8.17  0.04  0.39  0.02  0.08  0.03  
1.14  4.12  0.32  1.50  0.28  0.59  0.24  
1.54  5.36  0.20  0.74  0.07  0.14  0.18 
When substructures are absent, both harmonic power spectra of and are dominated by the second order mode, which is characteristic of a projected elliptical potential. The situation is totally different when substructures are taken into account. First, we notice that first mode () increase for lenses and . This is due to the fact that we kept the same definition of the mass center between the three lenses. The random position of subtructures generates a non zero impact parameter which affect the first order mode according the relation (37). Moreover, since substructures tend to break the ellipticity of the halo center in the present case, one expects that the second mode decreases. However, the most interesting feature is that modes corresponding to increase when substructures are present.
4.3 Practical limitations
A fraction of the error is produced by the raytracing simulation as well as the limitation of the considered resolution. Let us therefore consider a toy halo with an isothermal profile:
(46) 
where is evaluated so that the mass enclosed inside a radius kpc is . An isothermal profile is appropriate to estimate systematic error since it leads to exact solution with the perturbative method. Here again, we have evaluated the quantity by considering 100 realizations with different impact parameters. Sources have circular contour with and 15 millions particles have been used. For isophotes , which is supposed to have the higher error values, we have obtained % . Thus, in the following, we will consider that both raytracing method and the resolution limitation induce to a mean error of 0.3 % in contours reconstructions. Moreover, as we shall see in section 4, and can be characterized by their multipole expansion and their associated power spectrum (see equations 45). In Figure (6), we plot the amplitudes and as a function of n derived from the present experiments. These values for the amplitude () correspond to a noise that we have to take into account below. For this reason, we put a confident limit to .
In Aubert et al. 2007, the influence of smoothing , number of particles and opening angles have been extensively investigated on softened isothermal spheres. The number of particles and the models used here are similar and the parameters used in the current study can be considered as the most appropriate considering these previous tests (, , ). For instance a reanalysis of the Aubert et al. tests case present an average error in the deflection angle of 0.30.4 % at the Einstein radius. Considering that the profiles and the number of particles are similar, the current error estimation is in good agreement.
In the same studies, the (inverse) magnification reconstruction was also previously tested and the critical lines fluctuates around their theoretical location because of Poisson noise and rayshooting artefacts. By increasing the number of particles and by means of adaptive smoothing, these errors can be limited. Again, the reanalysis of these tests cases allows an estimation of the error on magnification of close to the Einstein’s radius. The inverse magnification is obtained from the joint calculation of the convergence and the shear through the same rayshooting technique. Hence, the error estimation does not rely on some propagation procedure but on the effective calculation and thereby includes Poisson sampling effects and rayshoot errors.
To finish, it’s interesting to determine how the error of 0.3 % in arc reconstruction is reflected in error in the image length. Let’s consider again the configuration studied in Sect. 3.2 and taking equations (14) and (17), one obtains:
(47) 
Taking with the same typical value, , and we have:
(48) 
As in the previous calculation, the actual size of the perturbation due to the error is obtained by direct multiplication with the Einstein radius. We take the same numerical value, , thus:
(49) 
Given the amplitude of morphological effects, this noise is not a source of concern, but had we considered astrometric effects, this noise would become be a real problem.
4.4 Statistics
Figures (7) and (8) show respectively the variations of the mean amplitudes and as a function of derived from our different samples of lenses. As mentioned above, we put a confident limit () and we exclude in our calculation all amplitudes below. First, when substructures are disregarded, power spectra are dominated by the second harmonic which just reflects the fact that our simulated lenses have a mean ellipticity. However, values of the fourth orders appear to be not negligible. This is probably due to the 2D projection which can lead to boxy projected densities. On the other hand, odd orders are negligible, as expected. When substructures are present, we can notice that the amplitudes of first harmonic, , have rather high values. As emphasized before, this is simply due to the fact the position of the lenses center is assumed to be the same for all objects. In fact, subtructures tend to modify the position of the mass center or, equivalently, tend to generate a non zero impact parameter which affect first harmonic coeficient according to the relation (37). However, the most interesting and important result is the presence of a tail in the power spectra which clearly suggests that the amplitude of high order harmonics () are not negligible anymore. Note that these effects are clearer when the density profile of substructures is modelled by a cusp as was the case for sample . Note that the power spectra of high orders harmonic () can be fitted by powerlaws:
(50) 
where the relevant fitting parameters are shown in the Table 4.
Having estimated the statistical contribution of substructures to the perturbative fields, and we may now address the problem of defining observational signatures of substructures. The observational effects are of two types: (i) an effect on the position of images which is controlled by the field , and (ii) effects on the image morphology, for instance its size in the orthoradial direction which depends on the structure of the field . An interesting point is that effects of type (i) are directly related to , while type (ii) effets is not related directly to but to its derivative (at least for images of small extension). Indeed, let us consider a source with circular contour; provided the image is small enough, a local linear expansion of the field will be sufficient to estimate the image morphology. We make the following field model:
(51) 
where is a position angle which should be close to the image center. At the edge of the images, we have:
(52) 
By defining , it follows that:
(53) 
Thus, scales like the orthoradial image size and can be directly related to observational quantities. The main difference with effects of types (i), which are related to the field rather than the field derivative, is that the derivative introduce heavy weights on the higher orders of the Fourier serie expansion of the field, where the substructure contribution is dominant. Indeed, the derivative of the Fourier serie introduce a factor at order , which translates in a factor on the components of the power spectrum. Thus the image morphology (and in particular its orthoradial extension for smaller images) will be much more sensitive to substructure than the average image position. To illustrate this, we have plotted in Figure (9) the amplitudes for lenses samples , , and . Moreover, to study the contribution of order , it is convenient to define the following quantities:
(55)  
(57)  
(58) 
Table 4 summarize mean values of relative to each catalogue of lenses. As expected, has higher values when substructure have a cusp profile. The total contribution of high order vary between and according the model of lenses used, which is quite significant.
Sample  

          
3.57  1.84  2.21  1.69  0.113 0.057  
7.33  2.07  6.65  2.21  0.179 0.097  
0.74  0.950  0.546  0.845  0.083 0.040  
2.08  1.575  1.462  1.468  0.095 0.049 
In closing, for the sample , mean errors between predicted contours from the perturbative method and raytracing solution are , and for isophotes , and respectively.
5 Conclusions
The structuration of matter on galactic scales remains a privileged framework to test cosmological models. In particular, the large number of dark matter subhalos predicted by the CDM cosmology is still a matter of debate. In this paper, potential signatures of substructures in the strong lensing regime were considered. This investigation makes use of the perturbative solution presented by Alard (2007, 2008), in which small deviations from the “perfect ring” configuration are treated as perturbations. In this framework, all informations on the deflection potential are contained in two onedimensionnal fields, and which are related to the radial expansion of the perturbed aspherical potential. The analysis of the properties of these two fields via their harmonic decomposition represents a simple and efficient way to track back possible signatures of subtructures. For instance, the perturbative method offer a simple and clear explanation of the multiplicity of the images. As shown in Alard (2008), the unperturbed image is a long arc for a cusp caustic configuration, but when a substructure is introduced near the arc, this latter is broken in 3 images. This effect is also illustrated in the present work in the example shown in the figure (2). Basically the breaking of the image and the change in image mutiplicity is due to the perturbative terms introduced by the substructure on the field . In the case of a circular source with radius , the image is broken when .
In this paper, lenses were modelled either by dark matter halos extracted from cosmological simulations or via toy models. The advantage of toy models is to reach a higher resolution and to allow us to study the influence of free parameters such as the inner profiles of substructures which are expected to play a central role here. We have first estimated the accuracy of the perturbative predictions by comparing the mean error between predicted contours images and theoretical ones derived via raytracing. We found that in general, the relative mean error is of the Einstein radius () when different impact parameters configurations and different source redshifts are considered. We also found that both resolution limitation and raytracing procedure lead to a systematic error of . Furthermore, although the accuracy of this approach for elliptical lenses is demonstrated in Alard (2007), we have checked and shown in the present study that the method works for the range of ellipticities derived from realistic 3D dark matter halos. This implies that the numerical evaluation of the coefficient of the perturbed potential, and at the Einstein radius is accurate enough to carry a statistical investigation.
We have generated several mock catalogues of lenses in which all objects have a total mass of and are at redshift . This value is motivated by comparison with observational surveys, and is close to where the strong lensing efficiency of clusters is the largest for sources (Li et al. 2005). The first catalogue represents our reference sample since all lenses are modelled by dark matter halos without substructures. In the other ones, substructures are described by either a cusp profile or a core profile. Their radial distribution is also a free parameter and we have used and . Our statistical investigation involves a MonteCarlo draw: the ellipticity of host halos, position of substructures, sources redshift for instance are randomly derived according to specific distributions. We found that the harmonic power spectra of and tend to develop a tail towards the large harmonics when substructure are accounted for. This effect is more pronounced when substructures have a cusp profile.
Several improvements of the present investigation are envisioned since the ultimate goal of the method is to provide a clear estimate of the amount of substructures in observations.

Statistically, the properties of the crosspower spectrum of and will be instructive. A clear characterisation of the covariance of these fields observed in Fig. 3 along with the result of Alard (2008) upon which, as high multipole order , about the same power is contained in each of the fields, would allow us to reduce the dimensionality of the problem and perhaps only consider either or .

In this vein, it remains to be shown to which extent the statistical properties of and can be approximated as Gaussian random fields. If so, the realization of mock giant arcs would be greatly simplified. Concerning this point, the number of experiments of the present work must be increased to provide a clear diagnostic.

A natural extension of this work would be to consider more realistic lenses while taking into account the dynamics of substructures inside the host halo, its connection to cosmology via the expected statistical distribution of substructures (see e.g. Pichon & Aubert 2006), as well as star formation mechanisms. Indeed, stripping proccesses caused by tidal forces may lead to more complex structures. For instance, the study of some merger events in the phasespace (radial velocity versus radial distance) reveals the formation of structures quite similar to caustics generated in secondary infall models of halo formation (Peirani & de Freitas Pacheco 2007). On the other hand, cooling processes and subsequent star formation may lead to steeper dark matter profile due to adiabatic contraction. Then, such processes should have an impact in the amplitude of high orders of our study.

Obviously a perturbation of the simple circularly symmetrical case needs not lie in the same lens plane as the main deflecting halo. Uncorrelated halos superposed along the line of sight to the background source may well introduce perturbations on top of substructures belonging to the main halo. This will contribute to some additional shot noise background in the power spectra of and , which has to be quantified and subtracted off using raytracing through large simulated volumes. This work is beyond the scope of the present analysis.

On the path to a possible inversion yielding and from observed arcs shapes and locations, several unknowns left on the rhs of Eq. (41) have to be controlled and will need to be fitted for in a nonlinear way in order to attempt a reconstruction of fields and . In addition, we only have considered a simple representation of the background source. A lumpier background source will translate into a less regular arc with some small scale signature in the observable quantities such as . However, the replication of these internal fluctuations along the arcs and, possibly, in the counter image, as well as the information contained in the various isophotes could allow us to reconstruct and directly from the observations. In this respect, Alard (2008) provides a general inversion method when two circular sources are considered for instance.
To conclude, the upcoming generation of high spatial resolution instruments dedicated to cosmology (e.g. JWST, DUNE, SNAP, ALMA) will provide us with an unprecedented number of giant arcs at all scales. The large samples expected will make standard lens modellings untractable and require the development of new methods able to capture the most relevant source of constraints for cosmology. In this respect, the perturbative method we present here may turn out to be a promising research line.
6 Acknowledgements
S. P. acknowledges the financial support through a ANR grant. This work made use of the resources available within the framework of the horizon collaboration: http://www.projethorizon.fr. It is a pleasure to thanks T. Sousbie, K. Benabed, S. Colombi, B. Fort and G. Lavaux for interesting conversations. We thank the referee for his useful comments that helped to improve the text of this paper. We would also like to thank D. Munro for freely distributing his Yorick programming language (available at http://yorick.sourceforge.net/) which was used during the course of this work.
References
 (2006) Amara A., Metcalf R. B., Cox T. J., Ostriker J. P., 2006, MNRAS, 367, 1367
 (2007) Alard, C. 2007, MNRAS, 382, L58
 (2008) Alard, C. 2008, MNRAS, 388, 375
 (2007) Aubert, D., Amara, A., & Metcalf, R. B. 2007, MNRAS, 376, 113
 (2004) Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376
 (1995) Bartelmann, M., Steinmetz, M., & Weiss, A. 1995, A&A, 297, 1
 (2007) Belokurov, V., et al. 2007, ApJ, 654, 897
 (2001) Bertschinger, E. 2001, ApJS, 137, 1
 (2001) Bode, P., Ostriker, J. P., & Turok, N. 2001, ApJ, 556, 93
 (2002) Bradač M., Schneider P., Steinmetz M., Lombardi M., King L. J., Porcas R., 2002, A&A, 388, 373
 (2004) Bradač M., Schneider P., Lombardi M., Steinmetz M., Koopmans L. V. E., Navarro J. F., 2004, A&A, 423, 797
 (2001) Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2001, ApJ, 548, 33
 (2001) Bullock, J. S., Kolatt, T. S., Sigad et al., 2001, MNRAS, 321, 559
 (2006) Cimatti, A., Daddi, E., & Renzini, A. 2006, A&A, 453, L29
 (2007) Chen J., Rozo E., Dalal N., Taylor J. E., 2007, ApJ, 659, 52
 (1996) Cole, S., & Lacey, C. 1996, MNRAS, 281, 716
 (2005) Cole, S., et al. 2005, MNRAS, 362, 505
 (1999) Colín, P., Klypin, A. A., Kravtsov, A. V., & Khokhlov, A. M. 1999, ApJ, 523, 32
 (2000) Colín, P., AvilaReese, V., & Valenzuela, O. 2000, ApJ, 542, 622
 (2002) Croft, R. A. C., Weinberg, D. H., Bolte, M., Burles, S., Hernquist, L., Katz, N., Kirkman, D., & Tytler, D. 2002, ApJ, 581, 20
 (2002) Dalal N., Kochanek C. S., 2002, ApJ, 572, 25
 (2001) de Blok W. J. G., McGaugh S. S., Bosma A. & Rubin V. C., 2001, ApJ 552, L23
 (2004) De Lucia, G., Kauffmann, G., Springel, V., White, S. D. M., Lanzoni, B., Stoehr, F., Tormen, G., & Yoshida, N. 2004, MNRAS, 348, 333
 () Diego, J.M.,Protopapas, P., Sandvik, H.B., Tegmark, M., 2005, MNRAS, 360, 477
 (2004) Dolag, K., Bartelmann, M., Perrotta et al., 2004, A&A, 416, 853
 (2006) Fedeli, C., Meneghetti, M., Bartelmann, M., Dolag, K., & Moscardini, L. 2006, A&A, 447, 419
 (2004) Gao, L., De Lucia, G., White, S. D. M., & Jenkins, A. 2004, MNRAS, 352, L1
 (2007) Gavazzi, R., Treu, T., Rhodes, J. D., Koopmans, L. V. E., Bolton, A. S., Burles, S., Massey, R. J., & Moustakas, L. A. 2007, ApJ, 667, 176
 (2004) Gentile, G., Salucci, P., Klein, U., Vergani, D., & Kalberla, P. 2004, MNRAS, 351, 903
 (1998) Ghigna, S., Moore, B., Governato, F., Lake, G., Quinn, T., & Stadel, J. 1998, MNRAS, 300, 146
 (2000) Ghigna, S., Moore, B., Governato, F., Lake, G., Quinn, T., & Stadel, J. 2000, ApJ, 544, 616
 (2004) Glazebrook K. et al, 2004, Nat 430, 181, astroph/0401037
 (2003) Hayashi, E., Navarro, J. F., Taylor, J. E., Stadel, J., & Quinn, T. 2003, Apj, 584, 541
 (2007) Hennawi, J. F., Dalal, N., Bode, P., & Ostriker, J. P. 2007, ApJ, 654, 714
 (1990) Hernquist, L. 1990, ApJ, 356, 359
 (2005) Horesh, A., Ofek, E. O., Maoz, D., Bartelmann, M., Meneghetti, M., & Rix, H.W. 2005, ApJ, 633, 768
 (2005) Kasun, S. F., & Evrard, A. E. 2005, ApJ, 629, 781
 (1993) Kauffmann G., White S. D. M. & Guiderdoni B., 1993, MNRAS 264, 201
 (2004a) Kazantzidis, S., Mayer, L., Mastropietro, C., Diemand, J., Stadel, J., & Moore, B. 2004, ApJ, 608, 663
 (2004b) Kazantzidis, S., Kravtsov, A. V., Zentner, A. R., Allgood, B., Nagai, D., & Moore, B. 2004, ApJL, 611, L73
 (2003) Keeton C. R., Gaudi B. S., Petters A. O., 2003, ApJ, 598, 138
 (1999) Klypin A., Kravtsov A. V., Valenzuela O. & Prada F., 1999, ApJ 522, 82
 (2008) Knebe, A., Draganova, N., Power, C., Yepes, G., Hoffman, Y., Gottloeber, S., & Gibson, B. K. 2008, astroph/0802.1917
 () Kochanek, C. S., Keeton, C. R., & McLeod, B. A. 2001, ApJ, 547, 50
 (2004) Kochanek C. S., Dalal N., 2004, ApJ, 610, 69
 (2005) Koopmans L. V. E., 2005, MNRAS, 363, 1136
 (2004) Kravtsov, A. V., Berlind, A. A., Wechsler, R. H., Klypin, A. A., Gottlöber, S., Allgood, B., & Primack, J. R. 2004, ApJ, 609, 35
 (2007) Leauthaud, A., et al. 2007, ApJS, 172, 219
 (2005) Li, G.L., Mao, S., Jing, Y. P., Bartelmann, M., Kang, X., & Meneghetti, M. 2005, ApJ, 635, 795
 (2006) Macciò A. V., Moore B., Stadel J., Diemand J., 2006, MNRAS, 366, 1529
 (2004) Mao S., Jing Y., Ostriker J. P., Weller J., 2004, ApJ, 604, L5
 (2007) Meneghetti M., Argazzi R., Pace F., Moscardini L., Dolag K., Bartelmann M., Li G., Oguri M., 2007a, A&A, 461, 25
 (2007) Meneghetti M., Bartelmann M., Jenkins A., Frenk C., 2007b, MNRAS, 381, 171
 (1999a) Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J. & Tozzi P., 1999, ApJ 524, L19
 (2005) Nagai, D., & Kravtsov, A. V. 2005, ApJ, 618, 557
 (2007) Natarajan P., De Lucia G., Springel V., 2007, MNRAS, 376, 180
 (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
 (199)) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
 (2003) Oguri, M., Lee, J., & Suto, Y. 2003, ApJ, 599, 7
 (2000) Palunas P. & Williams T.B., 2000, AJ 120, 2884
 (2007) Peirani, S., & de Freitas Pacheco, J. A. 2007, astroph/0701292
 (2008) Pereira, M. J., Bryan, G. L., & Gill, S. P. D. 2008, ApJ, 672, 825
 (2006) Pichon, C., & Aubert, D. 2006, 368, 1657
 (2008) Prunet, S., Pichon, C., Aubert, D., Pogosyan, D., Teyssier, R., & Gottloeber, S. 2008, astroph/0804.3536
 (2000) Salucci, P., & Burkert, A. 2000, ApJl, 537, L9
 (2006) Shaw, L. D., Weller, J., Ostriker, J. P., & Bode, P. 2006, ApJ, 646, 815
 (2000) Spergel, D. N., & Steinhardt, P. J. 2000, Physical Review Letters, 84, 3760
 (2001) Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726
 (2005) Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776
 (2005) Springel, V. 2005, MNRAS, 364, 1105
 (2003) Stoehr, F., White, S. D. M., Springel, V., Tormen, G., & Yoshida, N. 2003, MNRAS, 345, 1313
 (2006) Suyu S. H., Blandford R. D., 2006, MNRAS, 366, 39
 (2004) Vale A., Ostriker J. P., 2004, MNRAS, 353, 189
 (2007) van den Bosch, F. C., et al. 2007, MNRAS, 376, 841
 (03) Viel M., Matarrese S., Theuns T., Munshi D. & Wang Y., 2003, MNRAS 340, L47
 (04) Viel M., Haehnelt M.G. & Springel V., 2004, MNRAS 354, 684
 (1992) Warren, M. S., Quinn, P. J., Salmon, J. K., & Zurek, W. H. 1992, ApJ, 399, 405
 () Yoo, J., Kochanek, C. S., Falco, E. E., & McLeod, B. A. 2005, ApJ, 626, 51
 () Yoo, J., Kochanek, C. S., Falco, E. E., & McLeod, B. A. 2006, ApJ, 642, 22