A&A 506, 609621 (2009)
Perturbative reconstruction of a gravitational lens: when mass does not follow light
C. Alard
Institut d'Astrophysique de Paris, 98bis boulevard Arago, 75014 Paris, France
Received 7 March 2009 / Accepted 7 July 2009
Abstract
Aims. The structure and potential of a complex gravitational
lens is reconstructed using the perturbative method presented in Alard
(2007, MNRAS, 382, L58; 2008, MNRAS, 388, 375). This lens is composed
of 6 galaxies belonging to a small group.
Methods. The lens inversion is reduced to the problem of
reconstructing nondegenerate quantities: the 2 fields of the
perturbative theory of strong gravitational lenses. Since in the
perturbative theory the circular source solution is analytical, the
general properties of the perturbative solution can be inferred
directly from the data. As a consequence, the reconstruction of the
perturbative fields is not affected by degeneracy, and finding the best
solution is only a matter of numerical refinement.
Results. The local shape of the potential and density of the
lens are inferred from the perturbative solution, revealing the
existence of an independent dark component that does not follow light.
Conclusions. The most likely explanation is that the particular
shape of the dark halo is due to the merging of cold dark matter halos.
This is a new result illustrating the structure of dark halos on the
scale of galaxies.
Key words: gravitational lensing  dark matter
1 Introduction
In 1986 Lynds and Petrosian reported the discovery of giant arcs in Abell 370 and 224202, with a length of 100 Kpc and luminosities comparable to giant elliptical galaxies, but with bluer colors. Bohdan Paczynski (1987) interpreted these features, he suggested that the arcs may be the result of the clusters acting as gravitational lenses. The Paczynski interpretation was later confirmed by Grossman & Narayan (1988), and initiated a series of very fruitful developments on gravitational lenses. The lens inversion problem was tackled by Kochanek et al. (1989) using the ring cycle method. Other methods have been developed to solve the lens inversion problem, such as maximum entropy methods by Wallington et al. (1996), methods based on Fourier series expansion (Trotter et al. 2000), LensCLEAN (Wucknitz 2004), semi linear inversion methods (Warren & Dye 2003), Bayesian methods (Brewer & Lewis 2006), and potential reconstruction method (Suyu et al. 2009). The lens inversion methods may be classified in 3 types: (i) model dependent reconstructions; (ii) potential reconstruction on a grid; (iii) expansion of the potential in Fourier series. All methods have to consider the degeneracy of modeling, for instance Wayth et al. (2005) found that 6 different models were consistent with the observations for the lens ER 00472808. The problem of degeneracy in lens inversion is widespread. The degeneracy issue is related to the fact that the constraints on the potential are local. The value of the potential is constrained only in areas where images of the source are formed. Additional constraints comes from the fact that no images are formed in dark areas (Diego et al. 2005). Since there are no constraints on the potential in other areas, without additional assumptions, the problem is fully degenerate, and an infinite number of models are possible. One possibility is to explore the family of models that is consistent with the observations (Saha & Read 2009) and to try to find the common features between the models (nondegenerate quantities). However, in such an approach it is difficult to explore the full range of possible models. Another possibility to reduce the degeneracy is to use a local representation of the potential, like the methods (ii) and (iii). However, the grid methods or Fourier methods require many free parameters, and it is not clear whether the degeneracy problem has been improved. To solve the degeneracy problem it is essential to reduce the number of parameters, and to describe the potential in terms of fundamental, nondegenerate quantities. Since giant arcs resemble a part of a circle, it is natural to consider that arcs corresponds to a small break of the circular symmetry. In such an hypothesis, the fundamental properties of the arcs should be related to the quantities breaking the circular symmetry. These ideas were first explored by Blandford & Kovner (1988) who demonstrated that for a small impact parameter (with respect to the size of the critical circle) and small deviations from circularity, the equation for caustics depends only on the orthoradial component of the deflection field at the critical circle. Blandford & Kovner (1988) also presented a geometrical method to reconstruct the images of the lens. The same concepts were investigated later by Alard (2007) who derived an analytical equation to describe the formation of images by the lens. The main advantage of the analytical description is that it directly relates the morphological properties of the arcs to the potential expansion. This relation between the potential and the observational features suggests that the description of the potential is nondegenerate, and also simplifies the lens inversion problem. The reduction of the degeneracy problem in this approach is well illustrated by the fact that an infinite number of models corresponds to a given perturbative representation. Lens inversion in this perturbative framework was already tested in Alard (2008) using numerical simulations. It was found that the perturbative approach allows an accurate reconstruction of the potential for lenses with substructures. The concepts developed in this approach are very different from some recent perturbative numerical schemes (see for instance Vegetti & Koopmans 2009). In this type of work there is no attempt at reformulating the lensinversion problem in terms of fundamental quantities; the main goal is to derive an efficient numerical scheme to solve the ordinary lens equation. Peirani et al. (2008) showed that the perturbative approximation is accurate for realistic lens models. The numerical accuracy is of about 1% in units of the critical radius. The lens model of Peirani et al. (2008) were extracted from cold dark matter simulations (Horizon project), and some of them have a lot of structure (merging). This paper is organized as follows: the first part will present the perturbative method, and its basic properties. In an effort to make this work selfcontained, a synthesis of the other useful results on the pertubative method is presented in the appendices. This section also has the advantage of presenting all the equations essential to lens inversion. The second part of the paper will present the application of the perturbative method to the inversion of an interesting case of gravitational lens. This paper also will be an occasion to extend the technique of perturbative inversion and to relate the perturbative fields to the geometry of the lens.
2 The perturbative approach in gravitational lensing
In the perturbative approach (Blandford & Kovner 1988; Alard 2007),
gravitational arcs are interpreted as a perturbation of the perfect
ring situation. A perfect ring is formed when a point source is
situated at the center of a circularly symmetric lens. In the perfect
ring situation, the central point has an infinite number of images
located on the critical circle. There are 2 different kinds of
perturbations that breaks the symmetry of the system: the
noncircularity of the lens, and the source impact parameter.
When such noncircular perturbations are introduced, the symmetry of
the system is lost and the ring is replaced with a series of individual
images. However, as illustrated in Alard (2007) whatever the angular position
of these perturbed images, there is always an unperturbed image of the
central point on the critical circle at the same angular position. As a
consequence there is only a radial displacement
between the perturbed and unperturbed images. The unperturbed situation is represented by a circular lens with potential ,
and a point source with null impact parameter. In the perturbed situation, the source has an impact parameter
,
and the lens is perturbed by the noncircular potential
.
Both perturbations
and
are assumed to be of the same order :
Note that for convenience the unit of the coordinate is the critical radius. Thus by definition in this coordinate system the critical radius is situated at r = 1. As a consequence, in the continuation of this work all distances and their associated quantities (errors,..) will be expressed in units of the critical radius. In response to this perturbation, the radial position of the image will be shifted by an amount with respect to the unperturbed point on the critical circle. The new radial position of the image is:
The response to the perturbation can be evaluated by solving the lens equation in the perturbative regime. The lens equation reads:
The last step to solve the lens equation is to expand the potential is series of :
with:
By inserting Eqs. (1), (2), and (4) in the lens equation (Eq. (3)), it is straightforward to obtain an equation relating the response to the perturbation and :
and
This equation corresponds to Eq. (8) in Alard (2007).
Considering that the source has a mean impact parameter
,
the position in the source plane may be rewritten:
.
Assuming Cartesian coordinates (x_{0},y_{0})for the vector
,
Eq. (6) reads:
For a circular source with radius R_{0}, the perturbative response takes the simple following form (Alard 2007, Eq. (12)):
Equation (6) depends on . However this variable can be eliminated from Eq. (6) by renormalizing the fields: , and the source plane coordinates, (mass sheet degeneracy). These renormalized variables will be adopted in the rest of this work. The renormalization is equivalent to solving Eqs. (6) and (8) for . The variable will reappear when the renormalized quantities are converted to the original quantities.
3 Preprocessing of the HST data
The gravitational lens SL2SJ021400532 (Cabanac et al. 2007, SL2S public domain) was observed by HST in 3 spectral domains, F475W, F606W, F814W, with an exposure time of 400 s. The main difficulty in the preprocessing of the images is to remove the large number of cosmic present in the images. Cosmic are identified using a wavelet approach based on an estimation of the local scale in the image. Once a cosmic has been identified at a given position, the local area around this position is investigated to remove any other pixels connected to the cosmic. This procedure is based on routines from the ISIS package (Alard 2000). Note that since the number of cosmic is large in the images, the local density of cosmic sometime become too large in some areas. When too many nearby cosmic are detected within a small area, the area is flagged, and will be considered cautiously, or even rejected. Once the 3 images are cosmic cleaned, a recentering to a common grid is performed using the position of the bright objects. The last step in the processing of the images is to remove the background. There is some diffuse light coming from the galaxy halos at the level of the arcs, and the proper removal of this diffuse light is especially important. The diffuse background is removed in several steps. The first step is to build a general background model. To build this model, the image is divided in square areas of 32 pixel, the local background is estimated within each area using a cleaning method. The background value in each local area is then interpolated to compute a general background value at each position in the image. Once this background model has been subtracted from the image, the residuals around the arcs are investigated in order to control any systematic errors. A tiny systematic background residual is detected; this residual depends on the distance to the 3 main galaxies center of gravity. This background residual is related to an imperfect subtraction of a diffuse halo centered around the 3 main galaxies. An elliptical component centered at the center of gravity of the 3 galaxies is fitted to this residual. The best elliptical parameters are estimated by minimizing the deviations to elliptical symmetry. The mean radial profile of the residual is estimated, and finally the relevant elliptical component is subtracted. Note that this last background correction is small, it is only a few percent of the initial background subtraction. Thus, the 3 cosmic cleaned, recentered and background subtracted images are stacked, to produce the final reference image of the arc system. A color illustration of the arc system is provided in Fig. 1.
Figure 1: Color image for SL2SJ021400532 reconstructed from 3 noise filtered HST images in 3 bands. 

Open with DEXTER 
4 Lens inversion
This work presents the first successful model of the images formed by the lens SL2SJ021400532. This system was investigated by Limousin et al. (2009), who tried to fit this system using an elliptical lens model. They found that it was not possible to reproduce this system using such simple elliptical models. Limousin et al. (2009) were also confused by the color of the last image (right side in Fig. 1) which is slightly different to the color of the other images. They concluded that the identification of the images formed by this lens was uncertain, and as a consequence renounced modelling of this system. We will see in the continuation of this work that this lens requires models much more complex than elliptical, and that such models can be built naturally using the perturbative method presented in this paper. The development of this perturbative model also will naturally solve the problem of the different color of the last image. We will see that the first part of this image has the same color as the other images, while the remaining part of the image has a different color. This is a direct prediction of the perturbative model, which offers an elegant solution to the problems encountered by Limousin et al. (2009). We now turn to the analysis of this lens using the perturbative method. The lens SL2SJ021408053532 forms a system of 4 images of the source. The first 2 images are situated inside the large arc in the upper left (they appear as 2 identical mirror images). The 2 other image, are at the bottom and the right side of the frame (see Fig. 1). In the perturbative approach the problem of reconstructing this lens is equivalent to the reconstruction of the perturbative fields, and the next sections will concentrate on this issue. Obviously, the reconstruction problem is reduced to the reconstruction of the first order perturbative field only if the first order perturbative theory applies. The requirement for this approximation to work is that the gradient of the potential can be linearized in the vicinity of the critical circle (Alard 2007). For single dark matter halos numerical experiments shows that this local constraint on the gradient is usually met with good accuracy (Alard 2007; Peirani et al. 2008). It is only in the case where a local minima of the potential is present close to the critical circle that this approximation may not work. But such situations are very rare because they require that for instance a substructure falls very close to the critical line. However Alard (2008) showed that the general treatment of substructure in this approximation is correct even when the substructure is quite close to the critical circle (up to one tenth of the critical circle). Additionally, when an unlikely situation happens, such as a substructure lying on the critical line, the consequences are easily observable: a few additional images are formed near the perturbator. Thus such situation is easy to recognize, and is not apparent in the present data. In any case the best test of the first order approximation is to perform a successful reconstruction of the data, which means that in practise second order terms can be neglected. This will be the approach favored in this paper. Before the field reconstruction problem can be tackled, a necessary step is to evaluate the size of the critical circle, which will become by convention the coordinate system unit.
4.1 Estimation of the critical circle
The critical circle is estimated by fitting a circle to the mean position of the images. The center and radius of the circle are adjusted. The nonlinear adjustment procedure starts from a circle centered on the small galaxy at the center of the image. The initial estimate for the circle radius is the mean distance of the images to the center. A non linear optimization shows that the best center is close to the initial guess and that the optimal radius is close to 152 pixels; with a pixel size of the critical circle radius is . This value should be compared to the radius of found by Limousin et al. for the same lens. Note that the final result will not depend upon the particular choice of a given critical circle. Taking another circle close to this one would change a little the estimation of the perturbative fields, but the total background plus perturbation would remain the same.
4.2 General properties of the solution
The general properties of the solution can be derived from the properties of the circular source solution. To estimate the circular source solution, the outer contour of the source will be approximated with a circle. This approximation of the source may not be very accurate; the typical error on the field estimation will be of the order of the source outer contour deviations from circularity. However, this approximation is sufficient to obtain a first estimate of the solution and to derive its general properties. For circular sources the perturbative fields are directly related to the data (see Appendices A, and Eq. (8)). This direct relation between the fields and the data is possible only in areas of the lens plane where images of the source are present. In dark areas, the field will have to be reconstructed by interpolating the values in bright areas. The field will be approximated with the mean radial position of the images (with respect to the critical circle). The reconstruction of is simple and not ambiguous, but the reconstruction of requires more work. On the other hand is better constrained since no images should be formed in dark areas ( ). Furthermore, the field is of particular interest since it is directly related to the potential isocontour (Eq. (C.1)). The derivation of the general properties of require some specific guess of the local behavior of this field in the region of image formation. Images form in minima of and the example presented in Appendix A indicates that the behavior of near the minima is of 2 types. For small images has a linear behavior, while for larger images (caustics) the field behaves like a higher order polynomial. In the case of this particular lens system, we have 4 images of comparable size. The typical image size is quite small and it is reasonable to assume that the field behavior in the image is linear. Of course other models could be investigated, with for instance order 2 behavior of the field, but obviously these models are much more complex and would require higher order Fourier expansions (see Fig. 2). Additionally, adopting a higher order description (lower curve in Fig. 2) may introduce false features like local peaks. The locally linear model (upper curve in Fig. 2) is the simplest description consistent with the general properties of ; it has the smallest possible number of minima and maxima. Thus adopting the locally linear model guarantees that no unnecessary features are introduced in the model. On the other hand, in reality this description may be too simple. The validity of this model can only be tested by a direct application to the data, and an evaluation of the quality of the fit. A higher order description should be adopted only if the locally linear model does not fit the data properly. As a consequence, the locally linear model will be adopted in the continuation of this work. Once the general shape of the solution is known, everything is now just a matter of improving the numerical accuracy. The process of building the solution will proceed in several steps. First the local solution will be refined numerically, then some interpolation of the solution in dark areas will be computed. A Fourier expansion will be fitted to this piecewise polynomial solution. And finally the adjustment of the Fourier expansion to the data will be optimized using a non linear minimization procedure. Note that since only the square of enters Eq. (8), the sign of is not known, as a consequence two solutions with different sign will have to be explored. Hopefully, when the source is not circular this degeneracy of the sign can be broken.
Figure 2: Two example of possible shapes for the field . The black dots represents the mean image position. The red curve represents the local behavior of the field near the image, while the black curve is a smooth extrapolation. The upper curve present a solution where the field has local linear behavior near the images. In the lower curve a higher order behavior is presented. In dark areas the field must be larger that the source radius to avoid forming images (see Eq. (8)). Note also that the field must also have zero sum. 

Open with DEXTER 
4.3 Refinement of the local solution
The study of the general properties of indicates that the simplest model of this field is locally linear in the neighborhood of the images. This local linear model will be refined in this section. The same local linear approximation will be used for the field . As illustrated in the former section, the estimation will start from the circular source guess. For a circular source and a linear field model, the field is zero at the center of the image. This property will be assumed in the refinement of the local solution. The slope of the linear model for is estimated by solving Eq. (8). Note that solving Eq. (8) requires an estimation of the source radius R_{0}. Since the image maximum half radial width is R_{0}, an estimation of R_{0} can be obtained by averaging the maximum radial width of the 4 images. By adopting this linear local model of , taking the mean radial position as an estimation of linear parameters, and using the perturbative lens equation (Eq. (6)), it is possible to fold the 4 images to the source plane. Once folded to the source plane, the 4 images must be similar and their comparison is a direct test of the quality of the lens reconstruction. With the initial guess for the field a, correlation between the 4 folded images is effectively observed: the mean crosscorrelation between the images is 0.6. This cross correlation can be improved quite significantly by nonlinear optimization. The crosscorrelation between images in the source plane is maximized using the Nelder & Mead (1965) simplex method. The parameter space of the maximization is the 4 local slopes of the field, and also the 4 local slopes of the field . The simplex maximization brings the mean crosscorrelation between images to 0.8 which is a good indication that the linear approximation is already a good model. This maximization is also the occasion to build a first realistic model of the source, and to go beyond the circular source approximation. The result for the field is presented in Fig. (3). The local linear approximation is the dark blue line segments in Fig. (3).
Figure 3: An estimate of the field. The blue line segments represent the local linear guess at the image location, while the light blue curves are the second order interpolations between images. The red line is an adjustment of a 6th order Fourier series to the blue curve. 

Open with DEXTER 
4.4 Interpolating the local linear solution in dark areas
Since there are no constraints on the field in dark areas, the fit of this field in the full angular range is straightforward. Due to the constraint that no images are formed in dark areas, the situation is more complicated for the field . The local linear solution obtained in the former section needs to to be completed with second order polynomials to fill the gaps between the images. The condition that no images are formed in dark areas (Diego et al. 2005) is equivalent to the condition that the functional is larger in dark areas than in areas of image formation. Another constraint is that the sum of the field must be zero. And finally the field must be continuous and as smooth as possible. The first constraint is a constraint on the sign of the curvature of the local polynomial, or equivalently on the sign of the second order coefficient of the polynomial. The second constraint reduces the number of unknowns in the polynomial coefficients. And the last constraint settles the value of the remaining coefficients. Note that the constraints on zero and first order polynomial coefficients can be solved using linear methods, which reduce the number of unknown and leaves us with variables depending only on the second order coefficients. The remaining equations depending on the second order coefficients have to be solved with a constraint on the sign of the coefficients. A simple grid method was adopted to solve the relevant equations. The final solution is presented in Fig. 3 (light blue curve).
4.5 Final fitting of the fields
The final fitting of the fields will be conducted by reconstructing the images of the source and comparing these images with the HST data by chisquare estimation. This procedure will take into account the convolution of the images with the PSF, which was ignored in the former reconstructions. The HST PSF is estimated using the Tiny Tim software (Krist 1995). Note that the chisquare value between model and data can be interpreted as a comparison between the images. In the model, all images come from the same source, thus the chisquare directly measures the image similarity in the lens plane. In practice, the fitting procedure will work in the following way: start from a numerical guess of the field, reconstruct the source, then reconstruct the images of the source, and finally estimate the chisquare with respect to the HST data. An essential component in this procedure is the source reconstruction method, which is illustrated in the next section.
4.5.1 Source reconstruction
At this point the simplest method to estimate the source would be to use the Warren & Dye (2003) method. Unfortunately the source is very complex, and this method would require using too many basis functions, which in practice is untractable. The following method has been preferred: images are coadded in the source plane using an algorithm that preserves the geometrical properties of the images elements. The basic idea of this geometrical method is that pixels in the lens plane are divided in 2 triangles. These triangles are transported to the source plane and are assumed to remain triangles. The triangles in the source plane are projected on a grid and the flux is coadded. The optimal size of the grid is estimated by numerical experiments. The noise level in the source image is estimated by looking at a specific area, and only the significant pixels in the source grid are kept. The images of the source in the lens plane are computed using the same geometrical method, but this time in the opposite direction. The final step is to compare the images of the source with the data. The drawback of this method is that the pixels are affected by the convolution with the PSF, and thus the source reconstruction is not perfect. However, it is easy to correct this effect by using the following procedure: first ignore the effect of the PSF on source reconstruction, find a solution, and then correct the effect using a deconvolution procedure very similar to the VanCittert method (Van Cittert 1931). This deconvolution of the source starts from the solution obtained ignoring the PSF convolution. The initial guess is convolved with the PSF and a difference image is formed by subtracting the data. This difference image is taken as a new input to the source reconstruction procedure in order to produce a first correction of the initial guess of the source. This procedure is further iterated until only noise is present in the difference image. In practice this method requires to control the noise in the difference image. In this work a 3 level was adopted; tests were performed with other thresholds, 2 and 4, with no significant change in the results. The lower cuts require more iterations and provide somewhat more noisy source reconstruction. However, since the areas where effective corrections are needed are quite small (these areas correspond to the bright sharp features in the images), the source reconstruction is not affected much. It is also important to note that these PSF corrections have very little effect on the reconstruction of the perturbative fields. By comparing the fit obtained by ignoring the PSF and correcting for the PSF, it was found that the difference is only a small fraction of the noise in the field reconstruction. Thus for this particular lens this procedure of PSF deconvolution is important only for the reconstruction of the source. This is due to the fact that for this particular lens the thickness of the arc system is quite large with respect to the PSF size, and as a consequence the PSF convolution has only a limited effect.
4.5.2 Fitting of a Fourier model
The practical implementation of the fitting procedure presented above now will be described. The first step is to choose a starting estimate. It could be the piecewise polynomial model presented in Sect. 4.6, but since the Fourier expansion of the fields is directly related to the multipole expansion of the potential (see Eqs. (B.1) (B.2) and (B.3)), it is more convenient to use Fourier series. The initial guess is obtained by fitting a 6th order Fourier series to the piecewise model of , and for by fitting a Fourier series to the mean radial position of the images. Starting from this initial guess, the source and image reconstruction procedures are iterated using the simplex method, in order to minimize the chisquare between the model and data. This procedure converges to the minimum of the chisquare in the parameter space. This minimum corresponds to the optimal Fourier coefficients. Note that as pointed out in Sect. 4.2, the sign of the initial guess is unknown and 2 solutions with different signs have been explored. Fortunately in this case the source is not circular and the 2 solutions have different chisquare, which solves the problem of the degeneracy of the solution sign. The final solution, the source reconstruction and the image reconstruction are presented in Fig. 5, and the field solution is presented in Fig. 4. The fact that the source is complex and has a number of small features (see Fig. 7) can be used to estimate the accuracy of the modeling. There are 3 bright features in the source that are visible in all images (red, yellow, blue in Figs. 6 and 7). The position of these features in the HST data and model has been compared for each image. To obtain the best possible accuracy in the comparison of the positions, the following procedure was used: a small image of each feature is extracted in the model frame, and this small image is crosscorrelated with the HST data at various positions around the feature. The crosscorrelation map obtained using this procedure is fitted with a polynomial to locate the map maxima, to derive a direct estimate of the shift between the features. The residuals are somewhat smaller for the first 2 images. In particular a small discrepancy is visible in the bright part of the third image (lower left image); the distance between the two images of the bright features in the source are closer in the model. The discrepancy is of about 2.4 pixels, which is about 1.6% of the critical radius. On average the mismatch between model and observations is 0.9% of . This estimation is very close to the results obtained by Peirani et al. (2008) who found that for realistic lens models the mean errors on the reconstruction of images using the perturbative approximation was about 1% of . Note that two other features are visible in the source (see Fig. 7); due to variable levels of amplification, noise, and blurring by other nearby features, the images of these features are not always clearly visible in the HST data. However, when they are visible in the images, their position and the correspondence between model and data are indicated in Fig. 6. Note that the caustic system associated with the solution is complex and includes a number of sharp features (see Fig. 9). The complexity of this caustic system is not easily visible in the images; this is due to the fact that the source is large with respect to the small features in the caustic system.
Figure 4: Best solution for the fields: red curve, , and black curve . The blue dashes indicates areas where images are present. 

Open with DEXTER 
Figure 5: Results from image reconstruction: HST data ( upper right), model ( upper left), data and model subtraction ( lower left), source ( lowerright). Note the reproduction of the source details in the images of the model. An image of the source has been superimposed on the model ( upper left), to illustrate the scaling. 

Open with DEXTER 
Figure 6: Correspondence between the features in the model ( right side) and the HST data ( left side). The color code of the relevant source feature is presented in Fig. 7. The green ellipses indicate areas where cosmic cleaning was not satisfactory. There are a very large number of cosmic in the images, and sometimes several of them falls nearly in the same area, making cosmic cleaning very difficult. 

Open with DEXTER 
4.6 Noise
The model and the data are compared in a small region around the arcs with total number of pixels N. The Poisson weighted difference R_{i} between the model at pixel i, M_{i} and the HST data D_{i} is very close to a Gaussian (see Fig. 8). Considering that the model has parameters, 25 parameters of the Fourier expansion of the fields, and the approximately 450 pixels of the source, the chisquare is, . Changing a little the size of the area, either reducing it by moving closer to the center of the arcs, or enlarging the area does not significantly change the chisquare value.
4.6.1 Errors on the Fourier coefficients due to noise in the data
The errors on the estimation of the Fourier coefficients
are evaluated by local linearization methods near the minima (see, for
instance Press et al. 2007). The intensity M_{i} is a function of the
Fourier parameters p_{n},
M_{i}=M_{i}(p_{n}). Assuming that p_{n} is very close to the solution and adding a small shift to the parameters
,
the intensity variation can be linearized.
The parameters can be estimated by a linear leastsquare fit from Eq. (9), and thus the relevant errors are also estimated by taking the diagonal elements of A = B^{1}, the normal leastsquare matrix invert B, . This error estimation indicates that the errors on the Fourier coefficients due to the Poisson noise are .
4.6.2 Errors on the reconstruction of the perturbative fields
In the former section it was shown that the accuracy in the
reconstruction of the images was close to 0.9% of .
This error
is directly related to the errors on the reconstruction of the
perturbative fields. To estimate the amplitude of the errors on the fields
it is particularly useful to consider circular sources.
For a circular source there is a direct relation between the errors on
the fields and the errors on image reconstruction.
Eq. (8) shows that for circular sources,
is the mean
radial position of the image. As a consequence, the errors on
are of the same order as the errors on the image position, which were
estimated to be of the order of a percent. The field
is related by Eq. (8) to the radial width W of the image:
As a consequence, the error on W and the error on are related:
Assuming the statistics of the variable W and are the same and differ only by a scale factor, the above equation shows that this scale factor must be unity. As a consequence, the errors on are of the same order as the errors on the image radial width W, which are similar to the errors on image reconstruction: . The corresponding errors on the Fourier coefficients of the field expansion of order m can be estimated by analyzing the variance of the expansion:
Note that here and . Taking the square of F_{i} and averaging over
Figure 7: The source features visible in the images (see Fig. 7). The 3 brightest features are outlined by red, yellow, and blue ellipses. These features are used to estimate the accuracy of the reconstruction. The part of the source at the lower left of the image is effectively visible only in one of the images (see Sect. 5 for more explanations). 

Open with DEXTER 
Provided it is assumed that the variance of , the variance of F_{i} reads
Note that if the Fourier expansion has no constant term, the mean value of F_{i} is zero ( ). In practice there is an additional constant term for the field f_{1}, however considering that the number of Fourier parameters is quite large (12 parameters), the variance due to this additional parameter will be neglected. As a consequence for both fields the errors on the Fourier coefficients is . This error is about 10 times larger than the error due to the Poissonian noise in the image, which is clearly negligible. As a consequence, in the continuation of this work, only the noise due to the perturbative approximation will be considered.
4.6.3 Effects of the errors on the potential and density estimation
Some interesting properties of the lens are directly related to the
perturbative fields as illustrated in Eqs. (C.1) and (D.2). The potential isocontour depends on f_{0}, and
it will be assumed that the variance of f_{0} is identical to the variance of
.
The variance on
can be estimated using the method presented in Sect. 4.6.2. The expansion for
read:
Thus
Using Eq. (10) the associated variance of is
As a consequence, the average error on the estimation of is . In the same spirit the errors on the density isocontours (Eq. (D.2)) are
In practice the expansion in the above formula is conducted to order , with in order to limit the amplitude of the noise. Some numerical experiments shows that for there is no significant changes to the isocontour shape, while on the other hand the noise grows quickly with . Thus taking , the resulting error is .
Figure 8: Histogram of the difference image between the model and data. The black curve is the histogram of the pixels in the difference image. For each pixel the difference has been weighted by the local Poisson noise expectation. The red curve is the theoretical Gaussian expectation for the normalized residuals. The closeness between the black and red curves indicates that the result of the fit is very close to optimal. The blue curve is the histogram of normalized residuals of the data without subtracting the model. This curve is plotted here for reference; it illustrates the effect of subtracting the model from the data. 

Open with DEXTER 
5 Prediction of image color
Figure 9 indicates that the last image (lower right position in Fig. 1) does not have the typical ``hatshaped'' color diagram observed for the 3 other images. But this feature is actually a prediction of the perturbative model. The model predicts that the first half of the last image comes from the same source area as the other images. The other half of the last image corresponds to an area of the source that has no counterpart in the large arcs, and only a tiny counterpart in the lower left image. This situation is due to the crossing of the caustic line by a faint part of the source which is visible in Fig. 9 in the lower left of sourcecaustic diagram. The part of the source situated inside the caustics corresponds to the first half of the last image, while the other half of this image corresponds to the faint part of the source outside the caustics. As a consequence, the first half of the last image should have the same color diagram as the other images, and this is exactly what is observed (see Fig. 9). The remaining part should have the same color as the faint parts of the source, and this is also what Fig. 9 illustrates. The color variations in the source are entirely related to the different colors of the bright and faint parts as illustrated in Fig. 9. The accurate prediction of the observed image colors is a confirmation of the model.
6 Shape of the dark component
The shape of the potential isocontour is given by Eq. (C.1). However, in practice, only is known, and Eq. (7) shows that the first order coefficients of depend also on the impact parameters. But these first order terms are related to the centering of the potential, not to its shape, for instance the ellipticity is related to order 2 Fourier terms. It is however possible to obtain some partial information about the centering of the potential by using Eq. (B.3) the first order moments within the unit circle do not depend on the impact parameters (they cancel out in the calculations). The numerical value of the inner moments is quite small, of the order of . The outer moments have to be of the same order, and thus the centering terms are not large. Note that equations have been normalized by the factor which depends on specific properties of the background. In practice acts as an unknown scaling factor on and thus on the potential isocontours. However, an important point is that this scaling degeneracy does not affect the shape of the potential, and for realistic halos is close to unity (Peirani et al. 2008). The potential isocontours are presented in Fig. 10.
Figure 9: The 3 color diagrams obtained by taking the sequence of raw flux in the 3 HST bands, as well as the sourcecaustic system ( lower right). The upper left presents the color digram of the last image (red), and of the 3 other images. The upper right presents the color diagram for the bright parts of the first 3 images (red curves), while the blue curves represent the faint parts of the source. The lower left plot, compares the color diagram of the first half of the last image (black curve), to the mean diagram of the 3 other images (blue curve), and also compares the color diagram of the second half of the last image (green curve), to the mean color diagram of faint source parts in the 3 other images (red curve). The last figure presents the source and the system of caustics in the source plane. 

Open with DEXTER 
Figure 10: The upper plot presents the potential isocontours, the lower plot presents the density isocontours. In the upper plot, the red curves represent the two 3contours of the potential isocontours. The dashed curve is the unit circle. The blue curve represents the local deviation from circularity of the potential. The distance between the blue curve and the unit circle is proportional to the deviation from circularity. The lower plot presents the two 3 isocontours of the density. In each plot the small elliptical contours represent the galaxies in the lens. The small black contours represent the galaxies. 

Open with DEXTER 
6.1 Mass to light relation
The unperturbed potential isocontour is a circle with radius r=1when a local perturbation is introduced, some distortions to the circular isocontours will be visible in the neighborhood of the perturbator. As a consequence, provided that the potential is generated by the visible matter, the distortions should correlate with the position of the galaxies. Since the center of the potential is unknown, the center of the unperturbed unit circle is also unknown, and will have to be adjusted. The adjustment of this local circle is performed in the neighborhood of points on the potential isocontour. The angular range of the circle adjustment is . The residual of this local unit circle adjustment is presented in Fig. 10. The local deviations from circularity are not consistent with the positions of the galaxies. The opposite is visible; the maximum deviations from circularity are observed in the lower part of the figure where no nearby galaxies are present. The average deviation from circularity of the potential in this area is 0.05 . According to the mean error on the potential estimated in Sect. 4.6.3, this is a 10 deviation. This result demonstrates the need for a dark component that does not follow light. This is an interesting point for theories that try to avoid dark matter by modifying gravity. A natural solution is to consider that the dark component was formed by merging of cold dark matter halos.
6.2 Inner and outer contribution to the potential
Note also that Eq. (B.3) relates the Fourier coefficients of the fields and to the inner (a_{n},b_{n}) and outer terms (c_{n}, d_{n}) of the multipole expansion. As a consequence, the expansion of can be separated into inner contributions (which is equivalent to taking c_{n} = d_{n} = 0 in the original expansion), and outer terms (a_{n}=b_{n}=0). A numerical estimation shows that about 75% of the total power spectrum of the field comes from inner terms, which demonstrates clearly that the observed effects are dominated by the inner contributions. Since the deviation from circularity of the potential in areas with no light are much above the noise (10), the result that mass does not follow cannot be affected by the 25% contribution of the outer terms.
6.3 Density
To estimate the density of the lens it is necessary to make
some assumptions. Making a specific assumption about the
background allows us to make an estimation of
and also
break the masssheet degeneracy. For nearly isothermal lens distributions, the density isocontours are given by Eq. (D.2).
Note that in this equation the terms corresponding to the
impact parameters cancel out
To optimize the signal to noise ratio, the isocontours presented in Fig. 10 are reconstructed by using an order 4 Fourier expansion (see Sect. 4.6.3 for details). Note that a method to reconstruct the isocontours of a nearly isothermal lens was also studied by Evans & Witt (2003). However, their description of the potential is less general since it contains only one angular functional, it is suitable only for potentials with constant isophotes.
7 Discussion
This paper describes the structure of the dark matter envelope of a small group of galaxies. It is shown that the distribution of the dark mass does not follow light (Fig. 10). This is direct evidence for an independent dark component at the scale of galaxies. Such new results about a system with a mass much lower than a typical cluster (for instance the Bullet Cluster) are essential. First because there has always been a problem with missing baryonic mass in clusters (for instance X ray gas); second because structures like the Bullet Cluster contains many galaxies, and their distances are not known with good enough accuracy to avoid some degeneracy (Clowe et al. 2006). The galaxies in this small group have very much the same colors, and as a consequence are probably very nearly at the same distance. Interpretating these results is not ambiguous and we are probably observing a merging event between two cold dark matter halos. Since the dark matter halos are much more extended than the bright cores of the galaxies, the interaction between the dark halos is already quite strong, while at the same time the stellar distribution is much less affected. However, by looking carefully at the outer part of the luminous halo's of the galaxies, some indications of interaction are visible. In particular, the lack of asymmetry in the outer contours of the light distribution indicates that the small galaxy at the center of the group has probably been stripped of its outer parts. This issue will be analyzed in more detail in a forthcoming paper. The analysis of this particular lens makes it clear that the combination of the perturbative method and good HST data open interesting new possibilities. Using this approach it should be possible to probe the structure of the dark component at the mass scale of galaxies for a large number of systems. Statistical results about the geometry and structure of the dark halos should be derived, and will probably offer the possibility to evaluate the amount of substructure in the halos (Alard 2008). It is clear also that the massluminosity relation could be analyzed for a number of other systems, and that the discovery of systems similar to this one would constrain even more a modification of gravity.
AcknowledgementsThis work is based on HST data, credited to STScI and prepared for NASA under Contract NAS 526555. I would like to thank François Sevre, Leslie Sage, J. P. Beaulieu, S. Mao, R. Blandford, and P. Alard for helping with this paper, and the referee for helping to improve the presentation of this paper.
Appendix A: Physical meaning of the perturbative approximation
For circular sources, Eq. (8) shows that
is the mean radial position of the image. The width of the image is
.
The condition for image formation is
,
thus images are formed near the minimas
of
.
The minimas of
are of two different
types: linear behavior of
near
the minima (small images), or a higher order polynomial behavior:
caustics. To illustrate these properties, let's take the simple case of
an elliptical potential, with ellipticity parameter :
Let's now develop this potential to first order in :
As a consequence, , (see Eq. (4)) and since the background potential must be critical at r=1, . Using Eqs. (5) and (7) the perturbative fields reads:
Placing the source on the Xaxis (y_{0} = 0), at an impact parameter Eq. (A.2) one obtains at : , n=1..3, which corresponds to a broad minima of near (see Fig. A.1). This minimum corresponds to the formation of an arc. This is a case of a cusp caustic, and behaves like an order 3 polynomial near the singularity. The counterimage in the opposite direction is much smaller and corresponds to a narrow minimum of . For this small image for n = 1 only, and the behavior of this field is linear near the image.
Figure A.1: Plot of the function ( and y_{0}=0 in Eq. (A.2)). Images are formed when . The nature of the image is related to the behavior of near its minima, the arc corresponds to a large minima, with cancellation of the first and second order derivatives of . For the small image, the behavior of near zero is linear, and corresponds to a much smaller minima of . 

Open with DEXTER 
Appendix B: Relation to multipole expansion
The expansion of the perturbative potential at a given radial position r
reads (Kochanek 1991):
The perturbative theory requires the evaluation of the potential gradient on the critical circle. Since the coordinate system has been renormalized so that the critical radius is situated at r=1, all the evaluation of quantities related to the potential will be performed at r=1. The coefficients (a_{n},b_{n},c_{n},d_{n}) are related to the projected density of the lens by the following formula:
Using Eqs. (B.1) and (B.2), and noting that: the fields and :
As a consequence, there is a direct relation between the Fourier expansion of the perturbative fields and the multipole expansion of the potential at the critical circle. This relation also demonstrates that the Fourier expansion of the perturbative fields is not equivalent to the Fourier expansion of the potential by Trotter et al. (2000). The expansion proposed by Trotter requires a number of additional coefficients with respect to the multipole expansion. On the other hand, Eqs. (B.1)(B.3) indicate that the perturbative expansion is directly related to the multipole expansion without needing any additional coefficients. It is also clear from Eqs. (B.1) and (B.3) that since the Fourier expansion of elliptical potentials is approximately of order 2 (see Eq. (A.1)), the corresponding Fourier expansion will also be of order 2. Higher order terms in the Fourier expansion of the perturbative fields will thus be related to the skewness, boxyness, etc. , of the potential. The first order terms are related to the centering of the potential.
Appendix C: Geometrical interpretation
The perturbative fields are directly related to the geometrical
structure of the potential. The potential isocontour near the
critical circle is easy to infer from Eq. (D.1). The
isocontour equation reads:
Searching for the isocontour near the critical circle and developing the former equation to first order in :
When no perturbation is introduced (f_{0}=0), the isocontour must be reduced to the critical circle ( ), as a consequence the right side of the equation must be zero too, and the former isocontour equation reduces to:
The potential isocontour depend only on f_{0}, while the isocontour variation depends also on f_{1}. Considering a contour with radial position r_{i}, defined by , the equation of a nearby contour is , or equivalently . There is no isocontour variation provided that is constant, which is equivalent to being constant. To first order in :
Considering that Eq. (C.2) depends only on f_{1} for purely isothermal backgrounds, and that CDM halos are quite close to isothermal in general, the isocontour variation is dominated by f_{1}.
Appendix D: Relation to the geometry of the density
The two perturbative fields are related to the potential by the following
expansion (Eq. (4)):
Note that since the expansion presented in Eq. (D.1) is of order 2 in . This expansion reduces to order 1 when introduced in the lens equation. By inserting Eq. (D.1) in the Poisson equation, taking , and developing to first order in , the density reads
The isocontour near the critical circle is defined by the equation . In the absence of a perturbation (f_{i}=0) the isocontour must reduce to the critical circle ( ). As a consequence the equation of the isocontour is:
For a nearly isothermal distribution, the higher order derivatives (n > 1) are close to zero. Assuming that these higher order derivatives are of order ( and ), the former equation reduce to:
References
 Alard, C. 2000, A&AS, 144, 363 [NASA ADS] [CrossRef] [EDP Sciences]
 Alard, C. 2007, MNRAS, 382, L58 [NASA ADS]
 Alard, C. 2008, MNRAS, 388, 375 [NASA ADS] [CrossRef]
 Blandford, R. D., & Kovner, I. 1988, Phys. Rev. A 38, 4028
 Brewer, B., & Lewis, G. 2006, ApJ, 637, 608 [NASA ADS] [CrossRef]
 Cabanac, R. A., Alard, C., DantelFort, M., et al. 2007, A&A, 461, 813 [NASA ADS] [CrossRef] [EDP Sciences]
 Clowe, D., Bradac, M., Gonzalez, A. H., et al. 2006, ApJ, 648, L109 [NASA ADS] [CrossRef]
 Diego, J. M., Protopapas, P., Sandvik, H. B., et al. 2005, MNRAS, 360, 477 [NASA ADS] [CrossRef]
 Evans, N., & Witt, H. 2003, MNRAS, 345, 1351 [NASA ADS] [CrossRef]
 Grossman, S., & Narayan, R. 1988, ApJ, 324L, 37 [NASA ADS] [CrossRef]
 Kochanek, C., Blandford, R., Lawrence, C., et al. 1989, MNRAS, 238, 43 [NASA ADS]
 Krist, J. 1995, ASPC, 77, 349 [NASA ADS]
 Limousin, M., & the SL2S collaboration 2008 [arXiv:0812.1033]
 Lynds, C., & Petrosian, V. 1986, BA&AS, 18, 1014
 Nelder, J., & Mead, R. 1965, Comp. J., 7, 308
 Paczynski, B. 1987, Nature, 325, 572 [NASA ADS] [CrossRef]
 Peirani, S., Alard, C., Pichon, C., Gavazzi, R., & Aubert, D. 2008, MNRAS, 390, 945 [NASA ADS] [CrossRef]
 Press, W., Teukolsky, S., Vetterling, W., et al. 2007, Numerical Recipes (Cambridge University Press)
 Saha, P., & Read, J. 2009, ApJ, 690, 154 [NASA ADS] [CrossRef]
 Suyu, S., Marshall, P., Hobson, M., et al. 2006, MNRAS, 371, 983 [NASA ADS] [CrossRef]
 Suyu, S., Marshall, P., Blandford, R., et al. 2009, ApJ, 691, 277 [NASA ADS] [CrossRef]
 Trotter, C., Winn, J., & Hewitt, J. 2000, ApJ, 535, 671 [NASA ADS] [CrossRef]
 Van Cittert, P. H., & Physik, Z. 1931, 69, 298
 Vegetti, S., & Koopmans, L. 2009, MNRAS, 392, 945 [NASA ADS] [CrossRef]
 Wallington, S., Kochanek, C., & Narayan, R. 1996, ApJ, 465, 64 [NASA ADS] [CrossRef]
 Warren, S., & Dye, S. 2003, ApJ, 590, 673 [NASA ADS] [CrossRef]
 Wayth, R., Warren, S., Lewis, G., et al. 2005, MNRAS, 360, 1333 [NASA ADS] [CrossRef]
 Wucknitz, O., Biggs, A., & Browne, I. 2004, MNRAS, 349, 14 [NASA ADS] [CrossRef]
All Figures
Figure 1: Color image for SL2SJ021400532 reconstructed from 3 noise filtered HST images in 3 bands. 

Open with DEXTER  
In the text 
Figure 2: Two example of possible shapes for the field . The black dots represents the mean image position. The red curve represents the local behavior of the field near the image, while the black curve is a smooth extrapolation. The upper curve present a solution where the field has local linear behavior near the images. In the lower curve a higher order behavior is presented. In dark areas the field must be larger that the source radius to avoid forming images (see Eq. (8)). Note also that the field must also have zero sum. 

Open with DEXTER  
In the text 
Figure 3: An estimate of the field. The blue line segments represent the local linear guess at the image location, while the light blue curves are the second order interpolations between images. The red line is an adjustment of a 6th order Fourier series to the blue curve. 

Open with DEXTER  
In the text 
Figure 4: Best solution for the fields: red curve, , and black curve . The blue dashes indicates areas where images are present. 

Open with DEXTER  
In the text 
Figure 5: Results from image reconstruction: HST data ( upper right), model ( upper left), data and model subtraction ( lower left), source ( lowerright). Note the reproduction of the source details in the images of the model. An image of the source has been superimposed on the model ( upper left), to illustrate the scaling. 

Open with DEXTER  
In the text 
Figure 6: Correspondence between the features in the model ( right side) and the HST data ( left side). The color code of the relevant source feature is presented in Fig. 7. The green ellipses indicate areas where cosmic cleaning was not satisfactory. There are a very large number of cosmic in the images, and sometimes several of them falls nearly in the same area, making cosmic cleaning very difficult. 

Open with DEXTER  
In the text 
Figure 7: The source features visible in the images (see Fig. 7). The 3 brightest features are outlined by red, yellow, and blue ellipses. These features are used to estimate the accuracy of the reconstruction. The part of the source at the lower left of the image is effectively visible only in one of the images (see Sect. 5 for more explanations). 

Open with DEXTER  
In the text 
Figure 8: Histogram of the difference image between the model and data. The black curve is the histogram of the pixels in the difference image. For each pixel the difference has been weighted by the local Poisson noise expectation. The red curve is the theoretical Gaussian expectation for the normalized residuals. The closeness between the black and red curves indicates that the result of the fit is very close to optimal. The blue curve is the histogram of normalized residuals of the data without subtracting the model. This curve is plotted here for reference; it illustrates the effect of subtracting the model from the data. 

Open with DEXTER  
In the text 
Figure 9: The 3 color diagrams obtained by taking the sequence of raw flux in the 3 HST bands, as well as the sourcecaustic system ( lower right). The upper left presents the color digram of the last image (red), and of the 3 other images. The upper right presents the color diagram for the bright parts of the first 3 images (red curves), while the blue curves represent the faint parts of the source. The lower left plot, compares the color diagram of the first half of the last image (black curve), to the mean diagram of the 3 other images (blue curve), and also compares the color diagram of the second half of the last image (green curve), to the mean color diagram of faint source parts in the 3 other images (red curve). The last figure presents the source and the system of caustics in the source plane. 

Open with DEXTER  
In the text 
Figure 10: The upper plot presents the potential isocontours, the lower plot presents the density isocontours. In the upper plot, the red curves represent the two 3contours of the potential isocontours. The dashed curve is the unit circle. The blue curve represents the local deviation from circularity of the potential. The distance between the blue curve and the unit circle is proportional to the deviation from circularity. The lower plot presents the two 3 isocontours of the density. In each plot the small elliptical contours represent the galaxies in the lens. The small black contours represent the galaxies. 

Open with DEXTER  
In the text 
Figure A.1: Plot of the function ( and y_{0}=0 in Eq. (A.2)). Images are formed when . The nature of the image is related to the behavior of near its minima, the arc corresponds to a large minima, with cancellation of the first and second order derivatives of . For the small image, the behavior of near zero is linear, and corresponds to a much smaller minima of . 

Open with DEXTER  
In the text 
Copyright ESO 2009