# The influence of the matter along the line of sight and in the lens environment on the strong gravitational lensing

###### Abstract

We investigate the influence of the matter along the line of sight and in the lens environment on the image configurations, relative time delays, and the resulting models of strong gravitational lensing. The distribution of matter in space and properties of gravitationally bound haloes are based on the Millennium Simulation. In our numerical experiments we consider isolated lens in a uniform universe model and the same lens surrounded by close neighbours and/or objects close to the line of sight which gives four different descriptions of the light propagation. We compare the results of the lens modeling which neglects effects of the environment and line of sight, when applied to image configurations resulting from approaches partially or fully taking into account these effects. We show that for a source at the redshift the effects are indeed important and may prevent successful fitting of lens models in a substantial part of simulated image configurations, especially when the relative time delays are taken into account. To have good constraints on the models we limit ourselves to configurations of four images. We consider eighty lenses and large number of source positions in each case. The influence of the lens neighbourhood and the line of sight introduces the spread into the fitted values of the deflection angles which translates into the spread in the lens velocity dispersion 4 per cent. Similarly for the lens axis ratio we get the spread of 10 per cent and for the Hubble’s constant of 6 per cent. When averaged over all lenses and all image configurations considered, the median fitted values of the parameters (including the Hubble’s constant) do not differ more than 1 per cent from their values used in simulations.

###### keywords:

gravitational lensing: strong and weak - large-scale structure of the Universe^{†}

^{†}pagerange: LABEL:1–LABEL:10

^{†}

^{†}pubyear: 2014

## 1 Introduction

The strong gravitational lensing has many astrophysical and cosmological applications (see e.g. Kochanek, 2006 and references therein for a review). The qualitative understanding of the majority of multi-image configurations attributed to strong lensing can usually be based on a model of single galaxy-lens in a uniform universe. The quantitative description requires more complicated models, starting with so called external shear (Chang & Refsdal, 1984), invoking a galaxy cluster (Young et al., 1981), or taking into account another galaxy in the lens vicinity (Koopmans et al., 1994, MacLeod et al., 2006). The influence of the mass distribution in the lens vicinity and along the line of sight has been investigated by many authors (Kochanek & Apostolakis, 1988, Keeton et al., 1997, Bar-Kana, 1996, Chen et al., 2003, Wambsganss et al., 2004, Wambsganss et al., 2005, and D’Aloisio & Natarajan, 2011 to cite few).

Photometric survey of several strong lens surroundings (Williams et al., 2006) shows that many of them lie in poor groups of galaxies and that other groups near the line of sight are not uncommon. Spectroscopic observations (Momcheva et al., 2006, Auger et al., 2007) give the distribution of the galaxies along the line of sight and allow more detailed study of their grouping and influence on strong lensing in several cases. The inclusion of the directly observed objects around the lens in modeling greatly improves the quality of fits.

Keeton & Zabludoff (2004) investigate theoretically the problem of the main galaxy close neighbours constructing a poor group of galaxies. They check the image configurations corresponding to various source positions behind the group, different group members playing the role of the main lens, and others playing the role of the environment. They thoroughly analyze the influence of lens environment on the values of the fitted parameters. They show that by neglecting the objects around the lens one introduces bias to the fitted parameter values, which plays the role of a systematic error.

In this paper we continue our investigation of the environmental and line of sight effects which influence the action of strong gravitational lenses using the results of the Millennium Simulation (Springel et al., 2005) from its online database (Lemson & Springel, 2006). We basically follow the approach of Jaroszynski & Kostrzewa-Rutkowska (2012) (hereafter Paper I) including also the time delays in our considerations. We attempt to quantify the influence of matter in the strong lens environment (ENV) and along the line of sight (LOS) on the results of modeling multiple image configurations with measured time delays. We simulate the propagation of light in four different ways. In the most simplified approach we include only the isolated strong lens in a uniform universe model. Other approaches include the lens environment, or the matter along the line of sight, or both. Using each of the approaches we simulate many multiple image configurations, and attempt to fit them with the same kind of simplified model. The rate of failure (i.e. the fraction of unsuccessful fits in each approach) measures the influence of the environment and the line of sight (or each of them separately) on the strong lens. The differences between the fitted values of model parameters and the parameters used in simulations give the estimate of the systematic errors introduced by the environment and the line of sight. Our goal is the comparison of various effects related to light propagation, not the improvement of strong lens modeling.

In Sec. 2 we describe our approaches to light propagation. Sec. 3 presents tools used to compare different models and the results of such comparison. Discussion and conclusions follow in Sec. 4.

## 2 Model of the light propagation

### 2.1 Ray deflections and time delays

The multiplane approach to gravitational lensing (e.g. Schneider & Weiss, 1988; Seitz & Schneider, 1992) using the results of the Millennium Simulation (Springel et al., 2005) and the non-singular isothermal ellipsoids (NSIE) as models for individual halos (Kormann et al., 1994; Kochanek, 2006) is described in Paper I. Here we augment it with the description of relative time delays.

The evolution of the matter distribution is given by the Millennium snapshots which correspond to several discrete epochs with given redshifts . We assume that for the Millennium cube of epoch adequately describes matter distribution. Thus a ray crosses perpendicular layers of matter of defined thickness cut from the Millennium cubes of different epochs. The cubes are randomly shifted and rotated to avoid effects of periodic boundary conditions of the Simulation (Carbone et al., 2008). Since there are several matter layers between the source at and the observer, they can be treated as thin, and may be represented as surface mass distributions projected into their middle planes.

The matter content of each cube is described as a background component representing matter density averaged on cells plus gravitationally bound haloes given by De Lucia & Blaizot (2007) and Bertone et al. (2007). For the background we calculate the gravitational force in 3D and then use its component perpendicular to a ray to obtain the deflection angle. For ray beams with the small opening angles of the major influence of each background cube is an almost constant deflection angle plus its small variation, which we describe as the action of the background convergence and shear , . (These parameters are defined for each layer separately.)

Each projected halo is represented as a difference between two NSIE distributions with the same characteristic deflection angles , axis ratios , and position angles, but different values of core radii , which makes its mass finite:

(1) |

(compare Paper I). The above formula gives the value of characteristic deflection for a halo of given mass and virial radius . (We use , which validates the approximation). We consider axis ratios which are distributed within with maximum probability at , loosely resembling the results of Kimm & Yi (2007). The position angles in the sky are random. Since the background contains the whole mass, including mass of haloes, the latter must be compensated by some negative density distribution. We use the constant surface mass circular disk of radius () for this purpose (for details see Paper I). A compensated halo does not deflect rays outside its radius, so only a finite number of haloes has to be included in calculations.

A ray coming to the observer from the direction in the synthetic sky, crosses the -th layer () at some position given as (Schneider & Weiss, 1988):

(2) |

where is the angular diameter distance as measured by an observer at epoch to the source at epoch , - the angular diameter distance to the same source measured by a present () observer. In calculations we use comoving distances, which in a flat cosmological model simplifies some of the expressions, but we skip these technical details here. We also apply more efficient recurrent formula of Seitz & Schneider (1992), equivalent to the above equation.

Knowing the light path, one can calculate the geometric part of the relative time delay (as compared with the propagation time along a null geodesics in a uniform universe model) using the formula (Schneider et al., 1992):

(3) |

where we consider a ray coming to the observer from the direction (which defines its earlier path, so all are known) and the factors represent time dilatation.

The deflection in each layer can be calculated as a gradient of the deflection potential, which is also a measure of gravitational time delay (Schneider et al., 1992):

(4) |

The potential is defined up to a constant . The cumulative gravitational time delay after crossing all the layers is:

(5) |

where again defines the path and we account for time dilatation. Finally:

(6) |

The above expression contains unknown additive constant. Only the difference in calculated time delays between two rays can have a clear physical meaning.

### 2.2 Light propagation in selected solid angles

Our backward ray shooting covers separate and randomly chosen maps in the synthetic sky in size. To propagate a wide beam of this size we need deflection angles in overlapping regions in all layers. For this purpose we calculate and store the deflection angles for the evenly spaced positions of interest in each layer. Similarly we store the gravitational time delays caused by each layer. With the help of interpolations one can obtain deflections and time delays for any position in any layer.

The deflection angle and the deflection potential are given by analytical formulae and for each grid point Eq. 4 holds exactly. Interpolations and finite differencing used in the numerical approach make this relation approximate. To check the selfconsistency of our methods we have numerically calculated the rotation of the deflection angle on the grid, for different positions, maps, and differencing spacings. The rotation of the deflection field is present, but its level is fairly low: , where the averaging is over many positions within a single map. The result holds for maps at different redshifts for spacings of . For spacings of the rotation is two times more important.

We choose to calculate the strong lens effects in the nonuniform universe model for sources at the redshift . We examine all haloes within the beam up to this redshift (typically few thousands of them) as lens candidates estimating their Einstein ring radii . We choose ten lenses with the largest Einstein rings as most likely to produce multiple images of a randomly positioned source and investigate them in detail.

When investigating multiple image properties we need zoomed maps of smaller parts of the sky. For better resolution we use finer grids with deflections interpolated from previous calculations with the help of bi-cubic spline, so the interpolated deflection derivatives are continuous. We expect the multiple images to lie within few Einstein radii from the halo centre. We check for the presence of other haloes inside the circle of the radius surrounding a dominating lens. If they are present we enlarge the region of interest including zones around all companions. Finally we repeat backward ray shooting inside a square on the sky overlapping the region of interest. The fine grids giving the deflection angles in consecutive layers encompass still larger areas, so one can follow majority of rays deflected off the main region. We keep control of the rays and if some of them leave the mapped region, we neglect the related lens in further considerations.

We store the result of the ray shooting as a vector and scalar arrays:

(7) |

where gives the positions in the source plane of rays apparently coming from the directions on the observer’s sky. Similarly denotes the total (gravitational plus geometric) time delay along the ray coming from the same direction. Superscripts , enumerate the rays. Calculations for each ray are based on Eq.(2) and Eq. (6) respectively.

### 2.3 Numerical experiments: different descriptions of light propagation

To estimate the importance of the influence of the matter along the line of sight (LOS) and/or the matter in the strong lens environment (ENV) on the lensing properties, we repeat the simulations four times. The most realistic approach takes into account both LOS and ENV. When we take into account the action of the strong lens and neglect all other haloes belonging to the same layer, the influence of ENV is lost, but LOS may still be included. In this way we check the effect of LOS alone. Similarly, assuming that rays are deflected only in the strong lens layer ( for other layers) the effect of ENV alone is modeled. Finally using the same strong lens but in a uniform universe (UNI) we have the clear case unaffected by neither LOS nor ENV, which may be used for comparison.

We use the prismatic transformation (Gorenstein et al., 1988) in all layers making the deflection at the middle points of all maps zero and transforming the time delays accordingly. This implies that removing the background and all other haloes in the major lens plane changes the backward ray paths beyond the plane only in the second order. Thus the rays in the LOS+ENV and LOS approaches travel through very similar surroundings and the comparison of lens modeling in these cases is sensible.

### 2.4 Synthetic multiple image configurations

The lens equation giving the dependence is stored in the array for many positions in the image plane and is sufficient for quite accurate interpolations. Using numerical methods we find the distortion matrix and its determinant on the same grid, and then the critical lines and caustics. Finally we find image positions for many source positions using iterative methods. (Compare Paper I). We consider only sources placed inside caustics since we are interested in multiple image configurations. The source positions are evenly spaced inside the region of interest. Once the images for a given source position are found, we calculate the corresponding flux amplifications and time delays, and store all such information for further analysis.

Different descriptions of light propagation produce different critical lines and caustic structures. While the positions on the observer’s sky can be easily compared independently of light propagation model, the same can not be said of source positions. Thus the comparison between image configurations corresponding to the same source position but obtained with different descriptions of light propagation is meaningless. Only the statistical properties of image configurations corresponding to given propagation model can be compared with one another.

## 3 Fitting the simulated strong lenses with simplified models

### 3.1 Simplified fits and their rate of success

Shear | no | no | no | yes | yes | yes |
---|---|---|---|---|---|---|

Delays | no | yes | yes | no | yes | yes |

no | no | yes | no | no | yes | |

model # | 1 | 2 | 3 | 4 | 5 | 6 |

LOS+ENV | 0.44 0.19 | 0.35 0.19 | 0.390.16 | 0.62 0.17 | 0.52 0.16 | 0.56 0.17 |

LOS | 0.64 0.21 | 0.60 0.20 | 0.63 0.19 | 0.71 0.19 | 0.66 0.18 | 0.67 0.19 |

ENV | 0.65 0.18 | 0.58 0.19 | 0.59 0.18 | 0.85 0.08 | 0.76 0.14 | 0.78 0.11 |

UNI | .995 .004 | .994 .004 | .995 .003 | .998 .001 | .996 .002 | .998 .001 |

Note: The table shows the dependence of the rate of acceptability of fits on the method of treating the light propagation (see the text for details), for models neglecting (“no”) or taking into account (“yes”) the external shear, neglecting or taking into account time delays, keeping fixed (“no”) or modeling (“yes”) the value of the Hubble’s constant. The required accuracy is the same for all models (, , , ). Models are numbered for further reference. |

In the case of pure elliptical, non-singular lens one would expect (in general) 1, 3, or 5 images. Our numerical method of finding the images starts from ray-shooting on a finite resolution grid with a small extended source as a target, so some images may be missed or unresolved. Also the external influence may change the number of images. Since configurations with more images probably better constrain the lens, we concentrate on cases with five or four images.

To estimate the influence of the galaxies in the lens vicinity and the matter along the line of sight on the properties of the strong lensing, we attempt to model all our simulated cases of multiple imaging using a single lens model in a uniform Universe. The single lens we are using in modeling is a non-singular finite isothermal ellipsoid used also in the simulations for each of the haloes. When tracing the rays in simulations we interpolate all deflection angles from earlier computed arrays. Single lens modeling uses analytically calculated deflection angles and their derivatives, so it serves also as a test of ray tracing simulations in UNI case.

We also try a more sophisticated lens model using the same non-singular finite isothermal ellipsoids with external shear. In this way the tidal influence of masses close to the rays is at least partially represented.

Models with included time delays information can be used to measure the Hubble’s constant (Refsdal, 1964). We treat the value of the Hubble’s constant as a free parameter in some of our models, to check if it improves the fits and (which is more important) to find the possible influence of LOS and/or ENV on the accuracy of measurements.

The lens is described by six free model parameters (, , defined in Sec. 2, the position angle in the sky, and two components of the lens position in the sky ; the other characteristic radius of the lens, is kept at constant ratio with : and is not free). The source has some (unobserved) position in the sky (another 2 free parameters). Models with the external shear have another 2 free parameters (, ), and the Hubble’s constant may also serve as a free parameter. Thus the model has up to free parameters.

For a image system we have observed positions and independent flux ratios. If all relative time delays are measured, we have another independent observations. Including lens position in the sky we have (without time delays) or (with time delays) independent observations. To have the positive number of the degrees of freedom in all considered cases, we need . In Paper I we considered all four and five image configurations. (The four image configurations, constituting about 5% of all, were cases, where our numerical method failed to find the weakest central image.) Here we neglect the fifth (weakest) image (even if found) since it is usually not observed.

In this work we require that our models reproduce the image positions with the accuracy of , the lens position with , the image flux ratios (measured in stellar magnitudes) with , and time delays (when taken into account) with the accuracy of . The rate of success of modeling depends on the chosen accuracy parameters (compare Paper I). The values we are using correspond to the typical errors of present day observations (e.g. Keeton & Zabludoff 2004), except for the high time delays accuracy, chosen for the reasons explained in Sec. 4.

We fit our models looking for a minimum of . With the exception of the term describing the errors in time delays modeling, we strictly follow the approach of Paper I. For completeness we repeat the formulae supplementing them with the time delay part:

(8) |

The first term, controlling the lens position, used in the model has the obvious form:

(9) |

where the subscript stands for “lens”.

Using the simplified model we can calculate the source positions related to observed image positions . The deformation matrix and magnification matrix can be calculated at each image position:

(10) |

The mismatch between and implies magnified mismatch between modeled and observed image positions, which gives for (Kochanek, 2006):

(11) |

The fitting statistic for flux ratios is given as:

(12) |

where we calculate the flux ratios relative to the first (brightest) image and express them in stellar magnitudes as . (The role of flux ratios is played by the ratios of corresponding lens magnifications.) We then compare the modeled and simulated . For the time delays we use in full analogy:

(13) |

Again we calculate the time delays relative to the first (brightest) image, comparing the modeled and simulated differences of time propagation. is the number of images taken into account.

In all cases the quantities with an extra superscript “” are taken from simulations and mimic the observed values, while the quantities without it result from modeling.

Following the approach of Paper I, we try to fit simulated image configurations with the simplified models. For the starting model we use the parameters of the lens used in the simulations. If the external shear is included in modeling, we use its estimated value (see the next subsection). Standard minimization algorithm gives the best model. We accept it as a valid solution if its is within the range:

(14) |

where is taken from the table of distribution; statistically in 95% of cases the valid models belong to this range of for given DOF.

When simulating image configurations we use sources evenly spaced within caustic region of each lens. Thus every configuration represents some given solid angle of possible source positions . Its value is the same for all configurations related to a given lens , but differs between the lenses. The magnification bias (Avni 1981; Peacock 1982) makes the strongly amplified sources more likely to be included in the observed sample. In general the effect depends on the shape of the source luminosity function; we adopt the simplification of Keeton & Zabludoff (2004) assuming power law form of this function with the power index , so the relative probability of including given image configuration to the sample is proportional to the total magnification :

(15) |

where is a normalization constant.

Using the probability distribution we calculate the rates of acceptance of models of image configurations obtained using different approaches to light propagation. The differences between individual lenses give the estimated errors of the calculated rates. Our lens population is probably too small to be representative and the rates of acceptance can be used to compare the different approaches, but have no absolute meaning. The results are shown in Table. 1. According to this table both ENV and LOS do have an impact on acceptability of fits and the latter seems more important. With the inclusion of the external shear majority of the image configurations can be reproduced successfully. The time delays present an independent constraint on the solution and lower the rate of acceptance of our models. The inclusion of as a free parameter slightly improves the rate of acceptance of fits.

### 3.2 External shear and convergence

We estimate the external shear in two ways. In the first approach the estimate of the external shear is done by removing the main lens and computing influence of its neighbours (ENV) and matter along the paths of the rays (LOS) in the weak lens approximation. After removing the main lens we calculate the deformation matrix including only terms linear in deflection angles derivatives: , where - is the identity matrix and denotes the first order terms:

(16) |

With the usual convention we identify convergence and shear components:

(17) | |||||

(18) | |||||

(19) | |||||

(20) |

Since we are interested in the mean shear acting in the region of the size similar to the main lens Einstein ring, we calculate the deflection angles derivatives in Eq. 16 using finite differencing with spacing . We repeat shear calculations for many locations inside Einstein ring. We obtain the distributions of estimated shear due to LOS and ENV for each lens. We check the results with ten times smaller spacing obtaining very similar results. In Fig. 1 we compare the distribution of estimated shear with the distribution of its fitted values. This approach to shear estimation loosely resembles the cosmic shear measurement, since in both cases the resulting influence on ray paths/image shapes is taken into account.

The method gives also the estimated values of the convergence. We are not using convergence in modeling, but we need it for the discussion of the results of Hubble’s constant fitting. The distributions of the estimated convergence due to ENV, LOS and both are shown in Fig. 2. The distributions are not Gaussian, so we characterize them by their median values () and the parameter such that 68% of results belongs to the range . The numbers are: (ENV), (LOS), and (LOS+ENV) respectively. Our simulations include only eighty lenses and their environments belonging to eight independent small regions of space, so the convergence distributions are not expected to be universal. In the LOS case the skewness and the range of resemble the results of the recent study by Shang & Haiman (2011), who investigate the influence of weak lensing on the spread in the observed signals from standard candles and standard sirens.

In the second approach, following Auger et al. (2007), we estimate the sources of shear taking into account contributions from individual haloes, which are close to the lens on the sky. We use SIS model for each halo. The combined effect is given as:

(21) | |||||

(22) |

We include halos within the circle of radius around the main lens. The subscript enumerates the haloes and denotes the layer it belongs to, is the deflection by a SIS halo with the virial velocity , is the separation, and the position angle of the -th halo relative to the lens.

The fits to image configurations do not reproduce the shear values estimated using the first approach exactly. The fitting gives results which depend on the image configuration used in the particular optimization process and the estimates depend on the particular ray beams used in calculations. Treating fitting and estimating as two methods of measuring the shear value, we see that they are not in conflict: the distributions of fitted and estimated values are similar for combination of all lenses (compare Fig. 1) and the averaged values for a single lens obtained with the two methods are within standard deviation of each other (Fig. 3). The shear values estimated by the second method seem to be less reliable: comparing the influence of the lens neighbours on the sky within radii of , , and we see the dependence of the result on the size of the region taken into account. Our maps are too small, to check still larger regions (even the radius is not applicable to all our lenses), however the dependence seen in the results suggests that the saturation expected in the limit of large radii has not been reached. The analysis of fitted and inferred shear values for a sample of real quad lenses is presented by Wong et al. (2011). According to their study the values of the lens model shears and environment shears are in general different and in three of six investigated cases they are inconsistent at greater than 95% confidence.

### 3.3 Simulated and fitted time delays

We define

(23) |

as a measure of the error in model time delays.

In Fig. 4 we present the probability distribution for the values of the time delays error. Our aim is to find whether successful fitting of image positions and flux ratios with a simplified model automatically improves the time delay errors. The thick solid lines in Fig. 4 show the distributions of initial time delay errors for all configurations and dotted lines - for configurations successfully fitted with models neglecting time delays. These initial values are calculated for models having the parameters of the lens used in simulations and for the expected value of the external shear. The thin solid lines show the distributions of time delay errors after the fit. Comparing the plots we see that a small initial error in time delays gives a better chance of successful fit. Fitting of image positions redistributes the time delay errors. Examining the numerical data we find, that fitting the image positions increases the fraction of models with acceptable time delay errors () by 10 – 20 per cent, but another 15 – 25 per cent still have too large. (Large fraction of these, but not all, can still be improved with fitting procedure using time delays). Our results show that in – of cases reproducing image positions and flux ratios implies also reproducing time delays with good accuracy.

### 3.4 The Hubble’s constant

The value of the Hubble’s constant defines the length unit in the Universe, , and the time unit, . Scaling the distances and sizes of all mass concentrations, but preserving their velocity dispersions, positions, orientations, and shapes, we get the same lens equation and the same image configurations. Since relative time delays scale as well, fits using time delays can be used to measure the scaling factor and so the Hubble’s constant (Refsdal, 1964).

Using our lens models with treated as a free parameter and applying it to all quad image configurations we obtain distributions of estimated Hubble’s constant values corresponding to different descriptions of light propagation, using different lenses and different image configurations. The UNI case shows the typical accuracy of our method. Despite the fact that an initial model has parameter values the same as used in simulations and is acceptable, the fitting program looks for a minimum of which usually lies close to but not exactly at the initial position. This affects all the fitted parameters.

The distributions of fitted Hubble’s constant values are not symmetric, with outliers at positions far from used in simulations, so we characterize them using median values and absolute deviation such, that 68% of results belong to the range. The distributions of are shown in Fig. 5 for models including the external shear (model 6 in Table 1).

The results of the Hubble’s constant fitting depend to some extent on the way the external shear and convergence are treated. In Table 2 we present the median values and absolute deviations of the fitted Hubble’s constant for models neglecting the external shear (column 1) and models taking it into account (column 2). We also check the distributions of the corrected value (see Keeton & Zabludoff, 2004, Gorenstein et al., 1988)

(24) |

where the convergence is estimated using Eq. 17 and for each lens its averaged value is employed. We do not introduce the convergence into our models but only correct the results of fitting according to Eq. 24.

neglected | fitted | corrected | |

LOS+ENV | 70.5 5.3 | 72.6 4.2 | 72.9 3.5 |

LOS | 70.5 6.4 | 73.1 5.0 | 73.4 4.1 |

ENV | 72.3 1.8 | 71.5 2.7 | 71.2 2.7 |

UNI | 72.8 0.8 | 72.8 1.1 | 72.8 1.1 |

Note: The table shows the dependence of the fitted median value and its absolute deviation () of the Hubble’s constant on the method of treating the light propagation in simulations and on the way the external shear and convergence are treated. The cases with zero shear and convergence (“neglected”), the shear treated as a free parameter and no convergence (“fitted”), and the fitted shear with corrected according to Eq. 24 (“corrected”) are presented in consecutive columns. Naming of rows follows Table 1. |

Statistically the fitted values are in agreement with their expected values. In the most realistic LOS+ENV case 68% of the results, after the correction using Eq. 24 have errors below 5%. The median estimated convergences are close to zero; the LOS case with shows the largest departure. For LOS+ENV and LOS cases the inclusion of the external shear improves the median fitted and its accuracy. The correction improves the accuracy in both cases and the median in the first case. For the ENV case the opposite is true. The convergence we are using is estimated with the method of subsection 3.2, which is not perfect (compare shear estimated and fitted values). This method cannot be applied to real lenses since it is not based on observed properties of the image configuration.

In our statistical approach we have neither studied the individual lenses, nor the individual image configurations. Some experiments suggest that the constraints on the value from a single configuration are rather poor. For many image configurations values of from a wide range may give fits of comparable quality.

### 3.5 Lens parameters

Modeling simulated image configurations with simplified models we get the lens parameters. Since we know the “true” parameters of the lens used in rays tracing, we are able to estimate the systematic errors introduced by LOS and ENV effects. The results are presented in the following figures. We investigate the lens axis ratio , and the characteristic deflection angle . With four image configurations the lens virial radius is not sensibly constrained and the same can be said of the lens virial mass defined in Eq. (1). The knowledge of fitted and allows for a good estimate of mass within a given radius.

For every parameter we find the ratio of its fitted (denoted by the superscript “fit”) to original value. The dotted lines in Figs. 6 – 7 denote the results for the models, which do not use the external shear as a parameter, the solid lines - for models using it.

In Fig. 6 we show the distributions of the ratios for models including or not the external shear and time delay data, and for different approaches to light propagation. (See the caption). As can be seen in this figure (and also in Fig. 7), the distributions are rather complicated and depend both on the lens model and the approach to light propagation. Figs. 6 and 7 suggest that the widths of the distributions increase systematically when going through the ENV, LOS, and LOS+ENV cases, but the absolute deviations plotted in Fig. 8 do not fully support this simplified view. The median value is always close to unity (with 1% accuracy).

In Fig. 7 we show the distribution of the ratios for the characteristic lens deflection angle. The median value is always close to unity (1% accuracy). The characteristic deflection angle is related to the velocity dispersion of stars in the lensing galaxy in the same way as in simpler SIS lens models, i.e. . Since (compare Fig. 8), the velocity dispersion derived from the lens modeling has typical relative error .

We also compare the influence of LOS and ENV on the errors in the fitted parameter values showing the absolute deviations of the fitted and ratios in one diagram. As can be seen in Fig. 8, the typical errors in both parameters are roughly proportional to each other.

## 4 Discussion and conclusions

In this paper we have followed the approach of Paper I, investigating the influence of the matter along the line of sight on the properties of strong gravitational lenses. The main difference is the calculation of propagation time in the simulations and inclusion of the relative time delays in the modeling, which may serve the analysis of the Hubble’s constant measurement.

As compared to Paper I, we have simplified the treatment of the background matter distribution (compare Sec. 2.1), effectively replacing a part of the numerical calculation by the use of the fitted analytical formulae. This approach was not employed in Paper I, so the numerical inconsistencies, measured by the rotation to shear ratio (up to ) were more important and could have affected its results, exaggerating the reported influence of the matter along the line of sight on strong lensing.

Our numerical experiments show that matter in the immediate lens vicinity and matter along the line of sight do influence image configurations and values of fitted lens parameters . If the Hubble’s constant is treated as a free parameter, its fitted value also depends on the matter outside the lens. These conclusions are not new (see e.g. Keeton & Zabludoff, 2004) and have been demonstrated to be true in several cases of real lenses (Williams et al., 2006, Momcheva et al., 2006, Auger et al., 2007). According to our results matter along the line of sight is more important than the lens immediate neighbours. This may cause problems in modeling multiple images and time delays of high redshift sources, where direct observations of the matter along the line of sight may be impossible.

We consider a “typical” strong lens and its environment using a synthetic model of the matter distribution based on the cosmological Millennium Simulation (Springel et al., 2005) with gravitationally bound haloes described in De Lucia & Blaizot (2007) and Bertone et al. (2007). We investigate eight randomly chosen small square regions on the synthetic sky (each on a side) and look for the ten strongest lens candidates in each region, assuming that a source is located somewhere in the background at the redshift . We artificially remove the lens neighbours and/or matter inhomogeneities along the line of sight to compare the resulting image configurations. Thus we obtain four different matter distributions with the same strong lens. To assess the importance of matter outside the lens we try to model the image configurations using a single lens in a uniform universe model. The rate of failure of this approach, as well as the spread of the fitted lens parameters serves as a measure of the influence of the matter outside the lens on its properties. Inspecting Table. 1 one can see that (statistically) matter along the line of sight is more important as compared to the lens close neighbours. The same conclusion follows from inspecting Figs. 6 – 7, which show the distributions of fitted lens parameters.

The inclusion of the time delays in modeling lowers the rate of acceptance of fits not more than 10%, despite the high postulated time measurements accuracy (). We have introduced this accuracy requirement because many of our image configurations give time delays of the order of couple of days. For them the time delay modeling and Hubble’s constant fitting with more realistic accuracy would be meaningless. The time delay measurement accuracy we are using in our calculations is probably beyond the present observational possibilities. On the other hand one can choose appropriate lenses with appropriate image configurations for the purpose of estimating . This however goes beyond the scope of this paper, where we examine a randomly chosen lens population as a whole.

Treating the value of the Hubble’s constant as a free parameter and using the simplified lens model with the external shear we get its median fitted value in agreement with the “true” value used in simulations (compare Table 2 and Fig. 5 in Sec. 3.4). The convergence caused by matter outside the lens, when averaged over all our investigated cases, has median value close to zero (Sec. 3.2), so the bias in fitted values should not be pronounced.

The expected convergence along a randomly chosen line of sight is zero since we measure it relative to the uniform universe model. The haloes are compensated, so for each of them the above statement is also true. Since the positions of haloes in space are correlated, the density of matter near any lens should be on average increased. On the other hand, when calculating the convergence we remove the contribution of the lens itself, which is dominant. Thus the external convergence value close to zero on average, as in the case of lens population investigated here, is not implausible. In a more realistic approach the convergence near the main lens should be positive due to the correlations in matter distribution. Simply increasing the radii of compensating disks would not help, since the influence of each individual halo would diminish, but the number of haloes taken into account would increase in inverse proportion. One may introduce some artificial limiting radius independent of masses and take into account only lenses within such limits. We avoid such approach since the influence of high mass, distant haloes would be lost and the compensation would not be complete. (A different method of compensating mass of haloes was presented in Jaroszynski & Kostrzewa, 2010).

Some of the lenses give fitted Hubble’s constant values far from the median. We have checked the accuracy of our method applying it to lenses surrounded by the uniform matter distribution and found that it gives per cent scatter in fitted . The outliers cannot be explained as a result of unusually large convergence. Probably the perturbations to the outliers are nonlinear and cannot be modeled as the external shear and convergence.

Our calculations show that the external shear used in simplified modeling can be roughly estimated by weak lensing measurements along the line of sight (compare Figs. 1, 3 and Eqs. 18, 19). On the other hand taking into account only galaxies visible on the sky close to the lens position (Fig. 3, Eqs. 21, 22) may be insufficient – we have checked such estimates using lens neighbours on the sky in circles of radii up to and found the dependence of the results on the size of the region. It shows the important role played by the large scale structure, which is difficult to include in modeling. More sophisticated methods using spectroscopic observations of galaxies along the line of sight (Momcheva et al., 2006, Auger et al., 2007) may give better estimates of the external shear.

In this paper we have limited our interest to configurations of four images. In five image configurations we neglect the fifth, weakest image, which is not observed in real lenses. (This is another difference in our approach as compared to Paper I.) Neglecting the fifth image we loose the information related to the lens centre. The nonsingular isothermal ellipsoids we are using have small, approximately constant surface density cores. The size of the core has impact on the fifth image brightness, so without the fifth image it cannot be fitted. In simulation we postulate that the core size is proportional to the lens virial radius, which defines the total mass of the lens. Again, without the fifth image the total mass cannot be fitted. Of course the mass inside the Einstein ring can still be estimated using the fitted value of characteristic deflection angle .

## Acknowledgments

We are grateful to the Anonymous Referee, whose critical remarks greatly improved the paper. The Millennium Simulation databases used in this paper and the web application providing on-line access to them were constructed as part of the activities of the German Astrophysical Virtual Observatory. We are grateful to Volker Springel for providing us with the smoothed Millennium density distribution in the early stage of this project. This work has been supported in part by the Polish National Science Centre grant N N203 581540.

## References

- Auger et al. (2007) Auger, M.W., Fassnacht, C.D., Abrahamse, A.L., Lubin, L.M., and Squires, G.K., 2007, AJ, 134, 668
- Avni (1981) Avni, Y., 1981, ApJ, 248, L95
- Bar-Kana (1996) Bar-Kana R., 1996, ApJ, 468, 17
- Bertone et al. (2007) Bertone S., De Lucia G., Thomas P.A., 2007, MNRAS, 379, 1143
- Carbone et al. (2008) Carbone C., Matarrese S., 2008, MNRAS, 388, 1618
- Chang & Refsdal (1984) Chang K., Refsdal S., 1984, A&A, 132, 168
- Chen et al. (2003) Chen J., Kravtsov A.V., Keeton C.R., 2003, ApJ, 592, 24
- D’Aloisio & Natarajan (2011) D’Aloisio A., Natarajan P., 2011, MNRAS, 411, 1628
- De Lucia & Blaizot (2007) De Lucia G., Blaizot J., 2007, MNRAS, 375, 2
- Gorenstein et al. (1988) Gorenstein, M.V., Falco, E.E., Shapiro, I.I., 1988, ApJ, 327, 693
- Jaroszynski & Kostrzewa (2010) Jaroszynski M., Kostrzewa Z., 2010, AcA, 60, 41
- Jaroszynski & Kostrzewa-Rutkowska (2012) Jaroszynski M., Kostrzewa-Rutkowska Z., 2012, MNRAS, 424, 325 (Paper I)
- Keeton et al. (1997) Keeton C. R., Kochanek C. S., Seljak U., 1997, ApJ, 487, 42
- Keeton & Zabludoff (2004) Keeton C. R., Zabludoff, A.I., 2004, ApJ, 612, 660
- Kimm & Yi (2007) Kimm T., Yi S.K., 2007, ApJ, 670, 1048
- Kochanek (2006) Kochanek C.S., 2006, in: Gravitational Lensing: Strong, Weak and Micro, Saas-Fee Advanced Course 33, Springer, Berlin
- Kochanek & Apostolakis (1988) Kochanek C.S., Apostolakis J., 1988, MNRAS, 235, 1073
- Koopmans et al. (1994) Koopmans L.V.E., de Bruyn A.G., Jackson N., 1998, MNRAS,
- Kormann et al. (1994) Kormann R., Schneider P., Bartelmann M., 1994, A&A, 284, 285
- Lemson & Springel (2006) Lemson G. , Springel V., 2006, Astr.Soc.Pac.Conf., 351, 212
- MacLeod et al. (2006) MacLeod C.L., Kochanek C.S., Agol, E., 2009, ApJ, 699, 1598
- Momcheva et al. (2006) Momcheva, I., Williams K., Keeton C., Zabludoff A.,2006, 641, 169
- Peacock (1982) Peacock, J.A., 1982, MNRAS, 199, 987
- Refsdal (1964) Refsdal, S., 1964, MNRAS, 128, 307
- Schneider et al. (1992) Schneider P., Ehlers J., Falco E.E., 1992, “Gravitational Lenses”, Springer-Verlag Berlin Heidelberg New York
- Schneider & Weiss (1988) Schneider P., Weiss A., 1988, ApJ, 327, 526
- Springel et al. (2005) Springel V., White S.D.M., Jenkins A., et al., 2005, Nature, 435, 629
- Seitz & Schneider (1992) Seitz S. Schneider P., 1992, A&A, 265, 1
- Shang & Haiman (2011) Shang, C., Haiman, Z, 2011, MNRAS, 411, 9
- Wambsganss et al. (2004) Wambsganss J., Bode P., Ostriker J.P., 2004, ApJ, 606, L93
- Wambsganss et al. (2005) Wambsganss J., Bode P., Ostriker J.P., 2005, ApJ, 635, L1
- Williams et al. (2006) Williams K.A., Momcheva I., Keeton C.R., Zabludoff A.I., Lehar J., 2006, ApJ, 646, 85
- Wong et al. (2011) Wong, K.C., Keeton C.R., Williams K.A., Momcheva I.G., Zabludoff A.I., 2011, ApJ, 726, 84
- Young et al. (1981) Young P., Gunn J.E., Oke J.B., Westphal J.A., Kristian J., 1981, ApJ, 244, 736