text:
# 1 Introduction

## 1 Introduction

It seems quite natural to expect that the harder a high-energy collision is, the higher is the number of fragments. One of the most tractable and widely explored processes where this phenomenon can be seen is deeply inelastic scattering (DIS) with a possibility to change a hard scale ($`Q^2`$, gauge boson virtuality) and to detect its influence (if any) on the hadronic invariant mass ($`W`$) ”efficiency”. The simplest measure of this efficiency is the multiplicity of secondaries.

In 1969 Yang and his collaborators , basing themselves on the “fragmentation picture“ of violent collisions, made a qualitative prediction: “…for larger values of the momentum transfer $`t`$, the breakup process favors larger multiplicities of hadrons“ (at fixed hadronic mass). Early searches for this effect were inconclusive in both theory and experiment .

In the framework of QCD a quantitative result has been obtained in Refs. : it appears that QCD gluon bremsstrahlung leads to an increase of the hadron multiplicities in DIS, with an increase of $`Q^2`$ at fixed hadronic mass $`W`$, but that this increase is very slow. The distinctive feature of this result is that $`n^{DIS}(W,Q^2)`$ has a finite limit at $`Q^2\mathrm{}`$ and $`W`$ fixed.

Later, another result has been claimed in Ref. , which predicted an infinite and quite rapid growth of $`n^{DIS}(W,Q^2)`$ with $`Q^2`$. In the course of the inference of this result it was supposed that the influence of the (non-perturbative) composite structure of the nucleon is negligible while in  it plays a key role in the slowness of the $`Q^2`$ dependence of the multiplicity at fixed $`W`$. There is even more trivial objection. On general grounds, the infinite growth with $`Q^2`$ at fixed $`W`$ is impossible because of the apriori kinematical bound

$$n^{DIS}(W,Q^2)\frac{W}{m},$$
(1)
where $`m`$ is some effective mass.

Quite recently, a weak dependence on $`Q^2`$ in the framework of the Dual Parton Model was mentioned in .

Experimentally a statistically significant effect of the slow growth of $`n^{DIS}(W,Q^2)`$ was established in $`\nu (\overline{\nu })p`$ interactions  and in $`\mu ^+p`$ interactions . The results of the EMC  have been described in the framework of QCD in Ref. .

However, subsequent measurements at HERA (H1) were interpreted as a practical $`Q^2`$ independence  of $`n^{DIS}(W,Q^2)`$ (for the current hemisphere in the hadronic c.m.s.), while H1  and ZEUS  reported quite fast $`Q^2`$ dependence for the current hemisphere in the Breit frame. It should be noted, however, that this last result concerns different bins in $`W`$ for changing $`Q^2`$ values. Anyway, the situation is controversial and therefore very interesting.

In this paper we give our own interpretation of the HERA data on charged hadron multiplicities in the current fragmentation region; as will be seen in the text below, these are in agreement with Yang’s general hypothesis  and our early QCD results  (see also the review ).

## 2 Hadronic Spectrum and Multiplicity in DIS

According to the factorization for inclusive spectra in DIS, the hadronic spectrum in DIS is represented by two terms:

$$\frac{dn}{dy}^{DIS}(W,Q^2)=_{x_0}^1\frac{dz}{z}w(x,z,Q^2)\frac{d\widehat{n}}{dy}(W_{eff},Q^2)+\frac{dn_0}{dy},$$
(2)
where $`x_0=x+(1x)(m_h/W)\mathrm{exp}(y)`$, $`y`$ is the rapidity of the detected hadron and $`m_h`$ is its mass. In Eq. (2) $`d\widehat{n}/dy`$ defines the hadronic spectrum in partonic subprocess, while the quantity $`dn_0/dy`$ describes the spectrum of the proton remnant. The latter does not contribute to the current fragmentation region of DIS at HERA energies.

Correspondingly, the average hadronic multiplicity in DIS is represented by

$$n^{DIS}(W,Q^2)=_{x_0}^1\frac{dz}{z}w(x,z,Q^2)\widehat{n}(W_{eff},Q^2)+n_0.$$
(3)
For small $`x`$ the weight $`w(x,z,Q^2)`$ in (2), (3) is of the form:

$`w(x,z,Q^2)`$ $`=`$ $`D_g^q({\displaystyle \frac{x}{z}},Q^2,Q_0^2)f_g(z,Q_0^2)`$ (4)
$`\times `$ $`\left({\displaystyle _{x_0}^1}{\displaystyle \frac{dz}{z}}D_g^q({\displaystyle \frac{x}{z}},Q^2,Q_0^2)f_g(z,Q_0^2)\right)^1.`$
As can be seen, the hadronic spectrum in partonic subprocess, $`d\widehat{n}/dy`$, and the hadronic multiplicity $`\widehat{n}`$ depend on the effective energy, which is smaller than $`W`$:

$$W_{eff}^2=\frac{zx}{1x}W^2.$$
(5)
In what follows, we shall work in the c.m.s. of the final hadrons. In terms of rapidity, the current region in the c.m.s. corresponds to

$$Y<y<0$$
(6)
(it is assumed that the proton goes in the positive direction).

In our papers  it has been established that the total hadronic multiplicity in the partonic subprocess of DIS is related to the hadronic multiplicity in $`e^+e^{}`$ annihilation:

$$\widehat{n}(W,Q^2)n^{e^+e^{}}(W)$$
(7)
(up to small NLO corrections, which decrease in $`Q^2`$).

In the partonic subprocess, the rapidity varies in the range

$$\widehat{Y}<y+y_0<\widehat{Y},$$
(8)
where $`\widehat{Y}=\mathrm{ln}(W_{eff}/m_h)`$ and

$$y_0=\frac{1}{2}\mathrm{ln}\left(\frac{1x}{1z}\right).$$
(9)
The quantity $`y_0`$ determines the rapidity of the centre of mass of the partonic subprocess in the centre of mass of the complete process. On integration over $`z`$, the region (8) is ”smeared” into the region

$$Y<y<Y,$$
(10)
with $`Y=\mathrm{ln}(W/m_h)`$.

The average value of the effective energy in (2), (3), available for particle production, appeared to be dependent on both $`W`$ and $`Q^2`$ :

$$W_{eff}^2\kappa (Q^2)W^2.$$
(11)
The efficiency factor $`\kappa (Q^2)`$, which stands in front of $`W^2`$ in (11), is much less than 1 and grows slowly in $`Q^2`$.

From formulas (3), (7) and (11), one can see that the rise of the average hadronic multiplicity in DIS has the same physical nature as in $`e^+e^{}`$ annihilation. For the first time this behaviour has been experimentally established by H1 in 1996 .

However, the QCD growth of $`\widehat{n}`$ is delayed in DIS by the bound-state effects and the slow QCD evolution of the structure function. This is why we predicted that the $`Q^2`$ dependence of $`\widehat{n}`$ at fixed $`W`$ should remain numerically weak at HERA energies .

It follows from (8), (9) that the centre of the spectrum is shifted to the region of positive rapidities and tends to zero at asymptotically high $`Q^2`$ :

$$y_0|{}_{Q^2\mathrm{}}{}^{}\frac{1}{\mathrm{ln}(\mathrm{ln}Q^2))}.$$
(12)
The hadronic spectrum in partonic subprocess in the c.m.s. of DIS has the form

$$\frac{d\widehat{n}^h}{dy}=n^{e^+e^{}}(W_{eff},Q^2)\overline{D}^h(W_{eff},y).$$
(13)
Normalization in the RHS of Eq. (13) is done in agreement with formula (3).

## 3 Hadronic Multiplicities in the Current <br>Fragmentation Region

To calculate the multiplicity of charged hadrons in the current fragmentation region, we have to define expressions of the quark distribution at small $`x`$, of the hadronic spectrum in the partonic subprocess as well as of the multiplicity of charged hadrons in $`e^+e^{}`$ annihilation.

For the quark distribution, we use an analytical expression from Ref. , in the case of soft initial conditions. As was shown in , at small $`x`$ it is in good agreement with the data on the structure function from HERA in the wide ranges of $`Q^2`$. Namely, at high $`Q^2`$, we have

$$D_g^q(z,Q^2)rI_1(t)\mathrm{exp}(d\xi /2),$$
(14)
with $`d=\beta _0+20N_f/27`$. The variable

$$t=2\sqrt{6\xi \mathrm{ln}\left(\frac{1}{z}\right)}$$
(15)
is related to the QCD evolution parameter

$$\xi =\frac{2}{\beta _0}\mathrm{ln}\left(\frac{\alpha (Q_0^2)}{\alpha (Q^2)}\right),$$
(16)
where $`\beta _0=112N_f/3`$ is the $`\beta `$-function in lowest order and

$$r=\frac{t}{2\mathrm{ln}(1/z)}.$$
(17)
The quark and gluon distributions from Ref.  obey the GLAP evolution equations .

The expression of the initial gluon distribution at $`z`$ closed to $`1`$ is chosen to have the following form

$$f_g(z,Q_0^2)|{}_{z1}{}^{}(1z)^{n_g}.$$
(18)
We have omitted constant factors in the RHS of Eqs. (14) and (18) as they do not influence our final results.

The spectrum of hadrons in the partonic process $`\overline{D}^h`$ was calculated by many authors. We use the expression from Refs.  ($`N`$ is a normalization factor):

$$\overline{D}^h(W,\zeta )=\frac{N}{\sigma \sqrt{2\pi }}\mathrm{exp}\left[\frac{1}{8}k\frac{1}{2}s\delta \frac{1}{4}(2+k)\delta ^2+\frac{1}{6}s\delta ^3+\frac{1}{24}k\delta ^4\right]$$
(19)
calculated in the variable

$$\zeta =\mathrm{ln}\left(\frac{W}{E_h}\right).$$
(20)
Here $`E_h`$ is the energy of the detected hadron.

The average value of $`\zeta `$, $`\zeta _0`$, and its dispersion $`\sigma `$ are given by the formulas:

$$\zeta _0=\frac{1}{2}\tau \left(1+\frac{\rho }{24}\sqrt{\frac{48}{\beta _0\tau }}\right)\left(1\frac{\omega }{6\tau }\right),$$
(21)
$$\sigma =\sqrt{\frac{\tau }{3}}\left(\frac{\beta _0\tau }{48}\right)^{1/4}\left(1\frac{\beta _0}{64}\sqrt{\frac{48}{\beta _0\tau }}\right)\left(1+\frac{\omega }{8\tau }\right),$$
(22)
where

$$\tau =\mathrm{ln}\left(\frac{W}{\mathrm{\Lambda }}\right)$$
(23)
and

$$s=\frac{\rho }{16}\sqrt{\frac{3}{\tau }}\left(\frac{48}{\beta _0\tau }\right)^{1/4}\left(1+\frac{\omega }{4\tau }\right),$$
(24)
$$k=\frac{27}{5\tau }\left(\sqrt{\frac{\beta _0\tau }{48}}\frac{\beta _0}{24}\right)\left(1+\frac{5\omega }{12\tau }\right),$$
(25)
$$\delta =\frac{\zeta \zeta _0}{\sigma }.$$
(26)
Here $`\rho =11+2N_f/27`$, $`\omega =1+N_f/27`$.

At low (effective) energies we use the fit of the low-energy data on multiplicity of charged hadrons in $`e^+e^{}`$ annihilation from Ref. :

$$n^{e^+e^{}}=2.67+0.48\mathrm{ln}W^2,$$
(27)
while for high energies ($`W_{eff}>10`$ GeV) we apply the fit from Ref. , which well decsribes $`e^+e^{}`$ data up to LEP energies:

$$n^{e^+e^{}}=1.66+0.866\mathrm{exp}(1.047\sqrt{\mathrm{ln}W^2}).$$
(28)
We have corrected (27) for a fraction of the charged particles from $`K_s^0`$ and $`\mathrm{\Lambda }(\overline{\mathrm{\Lambda }})`$ decays.

Figure 1 represents the result of calculations of charged multiplicity in current hemisphere of the c.m.s. by using formulas (2) and (14) (solid curves) in comparison with the H1 data from Ref. . As can be seen, our QCD predictions are in very good agreement with the data. These are quite compatible with a slow growth of $`n^{DIS}(W,Q^2)`$ in $`Q^2`$ at fixed $`W`$.

Figure 2 demonstrates a rapid rise of $`n^{DIS}(W,Q^2)`$ in the variable $`W`$ for different values of $`Q^2`$, which was predicted many years ago in Refs.  and seen previously in $`e^+e^{}`$ annihilation (the very values of $`Q^2`$ taken from ). Let us note that the H1 data presented in Fig. 2 (see Table 4 in ) do not correspond to some fixed values of $`Q^2`$, in contrast with the experimental points in Fig. 1.

In order to obtain multiplicity of charged hadrons in current region of the Breit frame, the c.m.s. spectrum (13) must be integrated in the region

$$Y<y<y_B,$$
(29)
where

$$y_B=\frac{1}{2}\mathrm{ln}\left(\frac{1+v}{1v}\right)\frac{1}{2}\mathrm{ln}\left(\frac{1}{x}\right).$$
(30)
The quantity

$$v=\sqrt{14x(1x)}$$
(31)
in the RHS of Eq. (30) is the velocity of the Breit frame in the c.m.s. So, $`y_B`$ corresponds to zero rapidity in this frame.

The results of our calculations of the multiplicity of charged hadrons in the current region of the Breit frame are presented in Fig. 3 as a function of $`Q^2`$, in comparison with the H1 data  (solid squares) and ZEUS data  (solid circles).

The theoretical curves in Figs. 1–3 correspond to the following values of parameters:

$$Q_0^2=1\text{ GeV}^2,\mathrm{\Lambda }=0.25\text{ GeV}.$$
(32)
The parameter $`n_g=6.1`$ in Eq. (19) is taken from one of the MRST sets of parton distributions .

It should be noted that the strong $`Q^2`$ dependence seen by H1 and ZEUS in the Breit frame has nothing to do with the $`Q^2`$ dependence of $`n^{DIS}(W,Q^2)`$ in the c.m.s. (see Fig. 1); to a large extent it has a kinematical origin. The point is that an increase of $`Q^2`$ at fixed $`W`$ is equivalent to an increase of $`x`$. As a result, the current region of the Breit frame (29) enlarges. Thus, the rapid growth of hadronic multiplicity in the Breit frame in $`Q^2`$ (at fixed W) reflects a strong increase of the hadronic spectrum towards the central region.

As for the increase of $`n^{DIS}(x,Q^2)`$ in $`Q^2`$ at fixed $`x`$ in the Breit frame, it has been found that it is similar to that in $`e^+e^{}`$ annihilation at high $`Q^2`$, while there is a discrepancy between DIS and $`e^+e^{}`$ data at low $`Q^2`$ . Within our approach, it can be understood as follows. On the one hand, the increase of $`Q^2`$ results in an increase of the height of the spectrum (because $`W`$ grows). On the other hand, the position of the spectrum, $`y_0`$ as defined in (9), tends towards the region of positive rapidities.

These two phenomena go in opposite directions. At $`Q^2`$ high enough, $`y_0`$ varies very slowly with $`Q^2`$ (12). As a result, the rise of hadronic multiplicity in the Breit frame is analogous to that in $`e^+e^{}`$ annihilation. At low $`Q^2`$, $`y_0`$ changes more significantly . This effect partially compensates the growth of the spectrum and there appears a significant difference between DIS and $`e^+e^{}`$ data.

Finally, it is interesting to analyse the case when $`x`$ increases while $`Q^2`$ remains fixed. In this, $`W`$ decreases, which results in a rapid decrease of the spectrum. At the same time, however, the current region in the Breit frame becomes larger in accordance with formulas (29), (30). The effects are of the same order but opposite in sign. The ZEUS data in the Breit frame (see Table 2 in Ref. ) show that there is a slow rise of the hadronic multiplicity in the current hemisphere in the variable $`x`$ at different fixed values of $`Q^2`$.

## 4 Conclusions

In this paper we have presented the results of QCD calculations of the multiplicity of charged hadrons in the current hemisphere of DIS. Both the c.m.s. and the Breit frame are considered. We have shown that the H1 data are in agreement with Yang’s hypothesis and our QCD predictions. Namely, the efficiency of high energy collisions does weakly depend on the hard scale of the process (momentum transfer $`Q^2`$). It means that at fixed energy the efficiency of a particle production in hard processes increases with the shrinking of the interaction region ($`1/Q`$, in DIS).

The observed rise of the hadronic multiplicity with $`Q^2`$ in the current region of the Breit frame has both a dynamical and kinematical origin. This is why a direct comparison of available DIS data in the Breit frame with $`e^+e^{}`$ data is not completely correct.

New data from HERA on the hadronic multiplicity in the current region of the c.m.s. as well as measurements of the total multiplicity in DIS as functions of two variables ($`Q^2`$ and $`W`$/$`x`$) would be very important.

## Acknowledgements

We would like to thank Professors G. Altarelli, S. Catani, A.B. Kaidalov, M. Mangano and G. Veneziano for useful discussions. One of us (A.V.K.) is indebted to Professors G. Altarelli, S. Catani and M. Mangano for their support. We also thank Professor L.M. Shcheglova for stimulating discussions of the HERA data and Professor A.L. Kataev for sending us information on parton distributions.


format: text

text:
# Detection of weak gravitational lensing distortions of distant galaxies by cosmic dark matter at large scales

## Wide-field imaging with control of systematic shape errors

We observed three “blank” (i.e., not containing any known mass concentrations) fields at 23h48m, +0057 J2000; 04h29m, -3618; and 11h38m, -1233 over a period of several years, using the Big Throughput Camera, an array of four large, blue-sensitive CCDs at the Cerro Tololo Inter-American Observatory’s 4-m Blanco telescope. Constructed originally for weak lens observations, this camera covers a 35 arcminute field of view with 0.43 arcsecond pixels. We took multiple 500 second exposures shifted by 5-7 arcminutes and combined them to cover a 43 arcminute field. Before combining, we took several steps to reduce systematic errors arising from the optical system. First, we registered all the images onto a common linear coordinate system, free of the known radial distortion of the telescope optics. We then used the shapes of stars, which are foreground point sources free of the gravitational lensing effect, to correct any additional anisotropies in the point-spread function (the response of the optical system to point sources), such as those due to astigmatism and guiding errors. As described below, our observations covered multiple wavelengths; this enables filtering for certain types of stars and for distant galaxies.

The shape of a star or galaxy can be described by its second central moments, $`I_{xx}\mathrm{\Sigma }Iwx^2`$, $`I_{yy}\mathrm{\Sigma }Iwy^2`$, and $`I_{xy}\mathrm{\Sigma }Iwxy`$, where $`I(x,y)`$ is the intensity distribution above the night sky level, $`w(x,y)`$ is a weight function, the sum is over a contiguous set of pixels defined as belonging to the galaxy, and the coordinate system has been translated so that the first moments vanish. The second moments can be combined to form a size, $`I_{xx}+I_{yy}`$, and two components of a pseudo-vector ellipticity, $`e_1(I_{xx}I_{yy})/(I_{xx}+I_{yy})`$ and $`e_22I_{xy}/(I_{xx}+I_{yy})`$, which vary in the range $`[1,1]`$ (ellipticity in its colloquial sense is the amplitude of this pseudo-vector, $`ϵ\sqrt{e_1^2+e_2^2}`$ with its range $`[0,1]`$). Traditional intensity-weighted moments are calculated with $`w=1`$, but this produces ellipticity measurements with noise properties that are far from optimal or even divergent. In cases of white noise the formal optimal weight for an elliptical source is a noise-free image of that elliptical source. In the absence of such an image, weak lensing measurements are generally made with circular Gaussian weights. We use an elliptical Gaussian as the weight function, which places more weight on the high-signal-to-noise inner parts of the galaxy image, and is nearly optimal for most point-spread functions and for typical exponential galaxy profiles.

The moments of the Gaussian weight ellipse are iterated (from initial values provided by unweighted moments) to match the size and shape of the object, in order to obtain the highest possible signal-to-noise and to insure that the measured ellipticity is not biased toward the shape of the weight function. This “adaptively weighted moments” scheme has been extensively tested on simulated and real data, and has been shown to be unbiased. On simulated data, the algorithm recovers a somewhat higher fraction of the artificially induced shear than does simple intensity weighting. However, our final results do not depend on this particular weighting scheme. Its real benefit lies in rejection of peculiar objects, the vast majority of which are overlapping galaxies seen in projection. If the final centroid of a galaxy differs significantly from the starting centroid, that galaxy is rejected. If the centroid does not shift, the galaxy is accepted (and probably suffers little contamination by its neighbor). Object candidates are also rejected if the centroid or the ellipticity fails to converge; if they are too near the edge of the image; if the size grows too large; or if the moments are negative. About one-third of candidates found by the detection software (which provides the initial unweighted moments, and can “detect” occasional noise peaks) are rejected. For candidates which survive, the measurement error in the final ellipticity is accurately estimated by propagating the Poisson photon noise through the moment equations.

We used foregound stars at many positions in the field of view to measure and correct for systematic ellipticity error. Stars are distinguished from galaxies by their clear separation at the bright end of a size–flux density diagram. We identified roughly 100 such stars on each exposure of each CCD and made a least-squares fit (with 3$`\sigma `$ clip) of a second-order polynomial to the spatial variation of their ellipticity components $`e_1^{}(x,y)`$ and $`e_2^{}(x,y)`$, which would be zero at all points in an ideal observation free of point-spread function anisotropy. Fischer and Tyson have shown that nonzero $`e_1^{}`$ and $`e_2^{}`$ can be cancelled by convolution with a small (three pixel square) flux-conserving kernel with ellipticity components equal and opposite to those of the stars. Simulations as well as weak lensing data on clusters of galaxies show that faint galaxy induced systematics are also removed in this process of circularising stars. We convolved each image with its resulting position-dependent circularising kernel, after which the stellar ellipticities show little variation as a function of position. We then combined the images by averaging with a 3$`\sigma `$ clip, and repeated the point-spread function rounding on the combined image (using roughly 1000 stars and a fourth-order polynomial in this case). Fig. 2 depicts the evolution of one of our worst raw images through this process.

## Catalogues of distant galaxies

We repeated the observing and image processing for each field in three wavelength bands centered on 450 nm, 650 nm, and 850 nm, and for two of the fields we also took 550 nm images. The mean exposure time at each wavelength was 3400 sec. In each field, we used standard software  to identify object positions and fluxes on the 650 nm image (which is the deepest image in each field), yielding roughly 150,000 objects per field. We have confirmed the robustness of the weighted intensity moments in our detected object catalogues by using different (FOCAS  with adaptive circular kernel ) detection and evaluation software. At each object’s position, we evaluated the weighted moments at each wavelength, retaining only the measurements which the iterative weighted moment algorithm did not flag as suspect. Measurements with small sizes ($`I_{xx}<1.0`$ or $`I_{yy}<1.0`$) were also excluded as suspect. The result is a list of multiple independent ellipticity measurements (and corresponding estimated measurement errors $`\sigma _i`$) for each object.

We then computed the best estimate of each galaxy’s ellipticity by averaging the remaining measurements, weighted inversely by their estimated errors. If either of the ellipticity components at any wavelength deviated from this mean by more than $`3\sigma _i`$, that wavelength was eliminated and the process repeated. This step thus eliminates individual galaxy ellipticity measurements at wavelengths at which objects were noisy or blended, and it also reduces the systematic errors because the images at different wavelengths do not share the same residual point-spread function anisotropy. Finally we rejected objects with $`ϵ>0.6`$ as likely to be blends of more than one object, and applied flux density criteria ($`1.8\mu Jy>F_\nu >0.11\mu Jy`$ through the 650 nm filter, 23–26 R magnitude) to select objects likely to be distant galaxies. We use these same selection criteria in calibrating the typical redshift of the background galaxies (below). The final catalogues contain about 45,000 galaxies in each field. A visual inspection of the final catalogues indicates that they are free of spurious objects such as bits of scattered light around bright stars.

These observed ellipticities must be corrected for the overall broadening effect of the point-spread function, which makes elliptical galaxies appear more circular even if the point-spread function itself is perfectly isotropic. To calibrate this effect, we took a deep image with a very small point spread (the Hubble Deep Field South) and convolved it to the resolution of our final images, which is 1.07–1.25 arcsec as measured by the full width at half-maximum at 650 nm. (The resolution on individual exposures, or “seeing” was better, but the stellar size is larger in the final image with systematic shape errors removed from the point-spread function.) While some isolated galaxies became broader and less elliptical as predicted, most merged with their neighbors, producing many more elliptical objects than predicted and preventing the construction of a clear relationship between observed and true ellipticity for individual galaxies.

Instead we calibrated the fraction of cosmic shear recovered, as a function of resolution, from the ensemble of galaxies matching our selection criteria. We induced a known shear into the Hubble Deep Field South, convolved to the desired resolution, applied the same galaxy measurement and selection routines (at 650 nm only), and measured the mean ellipticity of the resulting sample. We averaged over repeated shears in several different directions to assess the measurement errors. The ratio of induced to recovered ellipticity was $`4.5\pm 0.5`$, with no clear trend as a function of resolution. The lack of such a trend would be quite surprising for isolated galaxies, but the coalescence of galaxy images appears to be the dominant effect. In the fairly small range of 1.07–1.25 arcsec resolution, this effect does not change the recovery factor by more than the measurement error of 0.5, so we adopt 4.5 as an overall ellipticity recovery factor.

## Ellipticity correlations of distant galaxies

Miralda-Escudé  has defined two physically revealing ellipticity correlation functions. In this approach, the ellipticity components of a galaxy $`i`$ are calculated not with respect to the arbitrary $`x`$ and $`y`$ axes of the image, but with respect to the line joining it to another galaxy $`j`$ (Fig. 3). Averaging over all galaxies $`i`$ and $`j`$ separated by angle $`\theta `$ on the sky, the correlations $`\xi _1(\theta )e_{1i}e_{1j}`$ and $`\xi _2(\theta )e_{2i}e_{2j}`$ have a unique signature in the presence of gravitational lensing, explained in detail in Fig. 3. We recently reported the detection of a cosmic shear signal in the quadrature sum of these correlation functions.

Fig. 4 shows the ellipticity correlations for each of the three fields in the angular separation range 2-36 arcmin (top panels). The plotted errors indicate 68% confidence intervals determined from 200 bootstrap-resampled realizations of the final galaxy catalogue in each field. Note that the measurements in different angular separation bins are not statistically independent, but $`\xi _1`$ and $`\xi _2`$ are independent from each other, as are the three fields. At $`\theta =6.1`$ arcmin, the confidence that $`\xi _1>0`$ is 97%, 99.5%, and 99.5% for the three fields in the order shown in Fig. 4. Similarly, the confidence that $`\xi _2>0`$ at the same angular scale is 87.5%, $`>99.5\%`$, and 97% respectively. Some “cosmic variance”, or real systematic differences among fields of this size, is expected, but the statistical errors are too large to examine this effect. We plot the average over the three fields in the lower panels of Fig. 4, with $`1\sigma `$ errors in the mean derived from the variance among the fields (black points and errors). The signature of gravitational lensing by large-scale structure is evident: $`\xi _1`$ declines as the angular scale increases, but is positive at all scales, while $`\xi _2`$ matches $`\xi _1`$ at small scales but drops below zero at large scales. This result is robust: Similar, but lower signal-to-noise, profiles are obtained if we use unweighted moments or moments from the 650 nm images only.

We performed several tests for systematic errors. The effects of residual point-spread function anisotropy are demonstrated by plotting the correlation functions of the stars (blue in Fig. 4). These are far closer to zero than are the galaxy correlations. Only $`\xi _1`$ in the innermost bin has an apparently significant stellar correlation. To test the effect that this might have on the galaxy correlations, we computed the star-galaxy correlations $`e_{1,star}e_{1,gal}`$ and $`e_{2,star}e_{2,gal}`$ (green in Fig. 4). The star-galaxy correlations are extremely close to zero in this bin. There are also tests involving the galaxy sample alone. The cross correlation $`\xi _3\frac{1}{2}(e_{1i}e_{2j}+e_{2i}e_{1j})`$ should vanish in the absence of systematic errors (red in Fig. 4). The result is reassuringly close to zero. The plotted errors for $`\xi _3`$ can also be taken as an indicator of the statistical error associated with the number and distribution of galaxies included in the catalogues (but reduced due to the averaging of two functions in $`\xi _3`$). This estimate of statistical error agrees roughly with that shown for $`\xi _1`$ and $`\xi _2`$. Finally, the weak lensing signature disappears if we randomise the galaxy positions.

Apart from these null tests, there are also affirmative tests. One test is to take similar data centered on a cluster of galaxies of known mass. A 650-nm image of massive cluster at redshift 0.45, taken with the same camera and processed in the same way, exhibits correlation functions ($`\xi _1`$ and $`\xi _2`$) of the expected angular dependence and of larger amplitude than in any of the three blank fields, despite likely contamination of the galaxy sample by cluster members. $`\xi _3`$ also vanishes in this field. Another test involves inverting the background galaxy ellipticity distribution to yield a map of projected mass in each “blank” field. We find occasional mass concentrations which can often be identified with likely foreground clusters, but no linear features or pileups at the edges of the image which might indicate problems in the background galaxy catalogues. Furthermore, when a mass map is made using only those galaxies likely to be behind a serendipitous cluster (based on colour information), the lensing signal from that cluster increases markedly. This corroborates the idea that the correlation functions are accumulating over many sources and many overdensities spread throughout the line of sight. All these tests indicate that we have indeed measured cosmic shear in our “blank” fields and that contamination from surviving systematic error is low. We now turn to comparisons with theoretical predictions of this effect.

## Comparison with theoretical predictions

Ellipticity correlations increase strongly with background galaxy redshift, so we must first constrain the source redshift distribution $`N(z)`$. Very little is known about the redshift distribution of galaxies as faint as those used here, so we assume a simple model $`N(z)z^2exp(z/z_0)`$, and adjust $`z_0`$ to match weak gravitational lensing observations of a high-redshift galaxy cluster of known velocity dispersion (MS1054 at $`z=0.83`$). We observed this cluster with the same camera and telescope and reduced the data in the same way as for the blank fields, and compared the faint galaxy ellipticities (tangential to the cluster center) to that expected for a range of $`z_0`$. We found that $`z_0=0.5`$ was the best match.

This model $`N(z)`$ was used as input to a cold dark matter simulation code by W. Hu and J. Miralda-Escudé (see ref ), which computes the shear power spectrum and correlation function for any given cosmology, using the prescription of Hamilton et al. and Peacock & Dodds to calculate the mass power spectrum in the non-linear regime when the growth of gravitationally collapsed dark matter structures modified the mass spectrum. Results were obtained for three cosmological models and are plotted along with our seeing-corrected measurements on a logarithmic scale in Figure 5. Two current models were normalised to the microwave background fluctuations (COBE) at large angle and to local galaxy cluster abundance (assuming mass traces light) at small angle: an open universe with $`\mathrm{\Omega }_{\mathrm{matter}}=0.45`$ (orange in Figure 5), and a flat universe dominated by a cosmological constant $`\mathrm{\Lambda }=0.67`$ (green, solid line). The agreement between the data and the two viable cosmological models is impressive for a first measurement. For comparison purposes we also show the old standard cold dark matter flat cosmology (blue), which is only COBE normalised. (A full listing of the parameters used in these models is shown in Table 1.) To illustrate the effect of varying $`N(z)`$, we also plot the $`\mathrm{\Lambda }`$-dominated cosmology with $`z_0=0.3`$ (green, dotted line). Since our model $`N(z)`$ peaks at $`z=2z_0`$, this lowers the typical redshift from 1.0 to 0.6. Decreasing $`z_0`$ decreases the amplitude of the correlations, but has little effect on their shapes. The uncertainty in $`N(z)`$ implies a factor of several uncertainty in the amplitude of the correlation, and is by far the dominant calibration error.

Despite this uncertainty, COBE-normalised standard cold dark matter is ruled out by the measured values of $`\xi _1`$. While this is not surprising, it is the first cosmological constraint from wide-field weak lensing, and it agrees with several other methods which disfavor this model. The other two models are consistent with the data at the $`3\sigma `$ level. The indication of a low $`\mathrm{\Omega }_{\mathrm{matter}}`$ universe here is in agreement with a remarkable array of independent methods, including type Ia supernovae, cosmic microwave background anisotropies, cluster baryon fraction together with cluster mass (lensing) and primeval deuterium, and the age of the oldest stars coupled with the Hubble constant. However, the shape of $`\xi _2`$ is not a good fit to either of these two model cosmologies, which are based on a single power-law mass spectrum. If confirmed by further data, this would suggest the need for a more complicated mass spectrum.

This technique can further distinguish between open and $`\mathrm{\Lambda }`$-dominated universes if extended to the somewhat larger angular scales where those cosmologies predict $`\xi _2`$ will cross zero as shown in Figure 5. A survey of many $`2^{}\times 2^{}`$ fields now underway will rule out one or more of these cosmologies at the $`8\sigma `$ level at 10 arcmin angles ($`3\sigma `$ level for a differential measure of the slope of the power spectrum). Separating the background galaxies into discrete redshift bins based on multi-colour photometry will enable measurement of the ellipticity correlation (or equivalently the dark matter power spectrum) as a function of cosmic time; wide-field weak lensing surveys deep enough to identify galaxies at $`z2`$ and measure their shapes will constrain several cosmological parameters. Ultimately, the combination of all the power spectrum probes (lensing, cosmic microwave background, galaxy distributions, and peculiar velocities) will tightly constrain theories of the origins of fluctuations in the early universe and their growth into galaxies and large-scale structure.

### Acknowledgements

We gratefully acknowledge help from Wayne Hu and Jordi Miralda-Escudé on theoretical predictions of several cosmological models. We thank Steven Gentile for his artwork, and the staff of CTIO for their help with the BTC project and for their upgrading and maintenance of the delivered image quality of the Blanco telescope. Cerro Tololo Inter-American Observatory is a division of National Optical Astronomy Observatory (NOAO), which is operated by the Association of Universities for Research in Astronomy, Inc., under Cooperative Agreement with the National Science Foundation. BTC construction was partially funded by the NSF.

| Model (Fig 5 colour) | $`\mathrm{\Omega }_b`$ | $`\mathrm{\Omega }_{\mathrm{matter}}\mathrm{\Omega }_b`$ | $`\mathrm{\Omega }_\mathrm{\Lambda }`$ | $`H_0`$ | n | $`\sigma _8`$ | normalization |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Standard cold dark matter (blue) | 0.05 | 0.95 | 0 | 50 | 1.0 | 1.17 | COBE only |
| $`\mathrm{\Lambda }`$-dominated, flat (green) | 0.039 | 0.291 | 0.67 | 70 | 0.94 | 0.84 | COBE+clusters |
| Open universe (orange) | 0.045 | 0.405 | 0 | 65 | 1.01 | 0.71 | COBE+clusters |

These cosmological models were chosen in order to put our ellipticity correlation measurements in context. The old standard cold dark matter model in which the universe is nearly closed by cold dark matter, is also disfavored in other observations. Its rms mass contrast is normalized to that found 300,000 years after the Big Bang via the cosmic microwave background radiation fluctuations observed with the COsmic Background Explorer satellite ($`COBE`$). The other two models agree with a wide variety of observations, but only the cosmological constant ($`\mathrm{\Lambda }`$-dominated, flat) cosmology also agrees with the recent evidence from supernova studies for accelerated expansion. $`\mathrm{\Omega }_b`$ is the fraction of critical density in to ordinary (baryonic) matter. $`\mathrm{\Omega }_{\mathrm{matter}}`$ is the fraction in all matter (mostly dark matter), and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ is the fraction in dark energy (the cosmological constant). $`H_0`$ is the Hubble constant in units of km s-1 Mpc-1. The power spectrum $`P(k)`$ of mass density fluctuations is often plotted in terms of inverse size: the wave number k is inversely proportional to length. The parameter n is the slope of the primeval density power spectrum as a function of k: $`P(k)k^n`$. For a scale-free power spectrum of density fluctuations, n = 1. The parameter $`\sigma _8`$ is the current rms mass contrast in a random sphere of radius 8($`100/H_0`$) Mpc compared with that for numbers of galaxies. The choice of n=1 and $`COBE`$ normalization for standard cold dark matter results in too much mass fluctuation on galaxy cluster scales. By adjusting the slope n and current rms mass contrast $`\sigma _8`$, models can be forced to fit the rms mass contrast now on galaxy cluster scales as well as the $`COBE`$ normalization.


text:
# Systematic uncertainties in gravitational lensing models: a semi-analytical study of PG1115+080

## 1 Introduction

Gravitationally lensed systems are powerful probes of galactic potentials and the scale of the universe. The advantage over the traditional stellar dynamical method is that from the image positions we can measure the shape and the mass of the dark halo well beyond the half-light-radius of a faraway lens galaxy, whether it is virialized or not. The time delay between two images is a measurement of the difference in the length of the two bent light paths, and scales with the distances to the lens and source. So the Hubble constant can be constrained once the redshift of the lens and the source and the time delay are measured. This way of getting $`H_0`$ has the advantage that the underlying physics (general relativity) can be rigorously modeled. The limitation is that there is often a sequence of lens models that can fit the image positions.

Presently about twenty strongly lensed systems are known, half of them being quadruple-imaged quasars and half being double-imaged quasars (e.g., Keeton & Kochanek 1996). We will concentrate on quadruple-imaged systems. They are better constrained than double-imaged systems, since the lens model needs to fit more image positions and also the ratios of time delays between any two pairs of images. Quadruple systems typically involve a quasar source well-aligned with the center of the lens potential well. Presently only two such systems (PG1115+080 and B1608+656) have accurately determined image positions and time delays.

A significant amount of numerical computation is usually required to invert image positions to intrinsic parameters of the potential (cf. Schneider, Ehlers & Falco 1992). The degeneracy of the resulting potential is often not fully explored because of the need to cover a large parameter space, particularly for flattened potentials. Previous authors have often restricted their studies to isothermal or power-law spherical models (Evans & Wilkinson 1998 and references therein), and elliptical models (Witt & Mao 1997 and references therein) and other simple models (Kassiola & Kovner 1993, 1995) with or without external shear. The fully general non-parametric method, e.g., the pixelated lens method of Williams & Saha (2000), is very powerful in demonstrating the complete range of the degeneracy in the lens models, but it involves significant amount of numerical computation and does not provide a clear insight to the relations between the characteristic parameters of the lens. For these reasons it is still desirable to find analytical, yet general potentials, which allow a quick exploration of the model parameter space. For example, it is of interest to generalize the analytical work of Witt & Mao (1997) to non-elliptical lenses. It would also be interesting to find analytical expressions for the time delay in these general lenses, which could help us to understand how the radial profile and lens shape affect the predictions on the Hubble constant.

Here we study a broad class of analytical models with non-axisymmetric, non-elliptical shape and semi-power-law radial profile (§2), and show how to calculate the lens shape and radial profile parameters directly from the image positions (§3). We apply the models to PG1115+080 (§4) and show that the images can be fit perfectly by a large range of lens models (§5). We summarize our results in §6 and conclude with the implications on the Hubble constant.

## 2 Lens equation in a general class of models

### 2.1 Decoupling of angular dependence and radial profile

Any two-dimensional lens potential can be cast in the following form,

$$\psi (x,y)=\psi (\omega ),\omega =\omega (r,\theta ),$$
(1)
where $`\psi `$ has the dimension of square arcsec, $`(x,y)`$ defines a rectangular coordinate system (in units of arcsecs to the West and North of the lens galaxy center) and $`(r,\theta )`$ the corresponding polar coordinate system with

$$(x,y)=(r\mathrm{sin}\theta ,r\mathrm{cos}\theta ),$$
(2)
where $`\theta `$ is the position angle, counterclockwise from North. Unless otherwise specified we shall follow the notations of Schneider et al. (1992). Here $`\omega `$ is defined to have the dimension of radius, so that constant $`\omega `$ curves correspond to equal-potential contours, and define the shape and flattening of the potential. The radial profile of the potential is a smooth function $`\psi (\omega )`$ of the radius $`\omega `$. We can also define

$$m(\omega )\omega \frac{d\psi }{d\omega }$$
(3)
as the mass (in units of square arcsec) enclosed inside the radius $`\omega `$ (in units of arcsec). For example, we have $`m=\alpha \psi `$ for a power-law model $`\psi \omega ^\alpha `$ with slope $`\alpha `$. For a source at redshift $`z_s`$ and angular distance $`D_{os}(z_s)`$, and a lens at redshift $`z_l`$ and distance $`D_{ol}(z_l)`$, the physical mass $`M(\omega )`$ is related to $`m(\omega )`$ by the scaling

$$M(\omega )\pi \mathrm{\Sigma }_cm(\omega ),$$
(4)
where

$$\mathrm{\Sigma }_c\frac{c^2}{4\pi G}\frac{D_{ol}D_{os}}{D_{ls}}\left(\frac{1\mathrm{radian}}{206265\mathrm{arcsec}}\right)^239\frac{D_{ol}D_{os}}{D_{ls}}$$
(5)
is the critical density in units of $`M_{}/\mathrm{arcsec}^2`$, and $`D_{ls}(z_l,z_s)`$ is the relative distance of the lens and source and all distances are in units of parsec.

A light ray from a source at $`(r_s,\theta _s)`$, being deflected to a direction $`(r,\theta )`$ by the lensing galaxy with potential $`\psi (\omega )`$ located at the origin, will experience a time delay $`\mathrm{\Delta }t`$ given by

$$\mathrm{\Delta }t(r,\theta )=h^1\tau _{100}(z_l,z_s)\left[\frac{1}{2}r^2rr_s\mathrm{cos}(\theta \theta _s)+\frac{1}{2}r_s^2\psi (\omega )\right],$$
(6)
where $`h`$ is the Hubble constant $`H_0`$ rescaled to $`100`$ km/s/Mpc. A characteristic value for the time delay is

$$\tau _{100}(z_l,z_s)\frac{D_{ol}D_{os}}{D_{ls}}\frac{1+z_l}{c},\text{for }H_0=100\text{ km/s/Mpc}.$$
(7)
According to Fermat’s principle, the images lie at the minimum of $`\mathrm{\Delta }t`$, so the lens equation is given by

$`{\displaystyle \frac{}{r}}\mathrm{\Delta }t`$ $`=`$ $`rr_s\mathrm{cos}(\theta \theta _s)m(\omega )_r\mathrm{ln}\omega =0`$ (8)
$`{\displaystyle \frac{}{r\theta }}\mathrm{\Delta }t`$ $`=`$ $`r_s\mathrm{sin}(\theta \theta _s){\displaystyle \frac{m(\omega )}{r}}_\theta \mathrm{ln}\omega =0.`$ (9)
Interestingly, the radial part $`m(\omega )`$ can be eliminated by simply combining the two equations,

$$\frac{r_s\mathrm{sin}(\theta \theta _s)}{rr_s\mathrm{cos}(\theta \theta _s)}=\frac{_\theta \mathrm{ln}\omega }{r_r\mathrm{ln}\omega },$$
(10)
similar as in Witt & Mao (1997). This implies that there is a relation linking the image positions to the shape of the potential directly, independent of the radial profile.

### 2.2 Property of the image positions: the semi-hyperbolic curve

For simplicity we shall concentrate on self-similar models with

$$\omega (r,\theta )=rf(\theta ),$$
(11)
where $`f(\theta )`$ defines the shape of the equal potential contours. In principle the shape function can be bi-symmetric or lopsided as long as the corresponding surface density is positive everywhere in the lens plane. To be specific, we will restrict our discussions to bi-symmetric potentials with an angular part

$$f(\theta )=\left|1\delta \mathrm{cos}2\theta ^{}\right|^\beta ,$$
(12)
which is a three-parameter $`(\beta ,\delta ,\theta _p)`$ function of the angle $`\theta `$, where $`\beta `$ is a constant, the parameter $`\delta `$ is a flattening indicator, and $`\theta _p`$ is the position angle of a principal axis of the potential. The angle $`\theta ^{}`$ is the azimuthal angle $`\theta `$ except for a rotation with

$$\theta ^{}\theta \theta _p,(x^{},y^{})=(r\mathrm{cos}\theta ^{},r\mathrm{sin}\theta ^{})$$
(13)
so that $`(x^{},y^{})`$ defines a rectangular coordinate system with the axes coinciding with the principal axes of the lens. For example, the source at

$$(x_s,y_s)=(r_s\mathrm{sin}\theta _s,r_s\mathrm{cos}\theta _s)$$
(14)
in the original rectangular coordinate system would be at

$$(x_s^{},y_s^{})=(r_s\mathrm{cos}\theta _s^{},r_s\mathrm{sin}\theta _s^{}),\theta _s^{}\theta _s\theta _p,$$
(15)
in the rotated rectangular coordinate system.

With these we can compute the right-hand part of eq. 10,

$$\frac{_\theta \mathrm{ln}\omega }{r_r\mathrm{ln}\omega }=\frac{d\mathrm{ln}f(\theta )}{d\theta }=\frac{1}{v\mathrm{cot}\theta ^{}+u\mathrm{tan}\theta ^{}},$$
(16)
where

$$u\frac{1+\delta }{2\mathrm{\Delta }},v\frac{1\delta }{2\mathrm{\Delta }},\mathrm{\Delta }=2\beta \delta ,$$
(17)
are shape indicators just like $`\beta `$ and $`\delta `$. Substituting in eq. 10, and rewriting the image radius $`r`$ as a function of $`\theta ^{}`$, we find the images fall on a family of curves $`r=r(\theta ^{})`$ defined by

$$r(\theta ^{})=r_s\left[\mathrm{cos}(\theta ^{}\theta _s^{})+\mathrm{sin}(\theta ^{}\theta _s^{})\left(v\mathrm{cot}\theta ^{}+u\mathrm{tan}\theta ^{}\right)\right].$$
(18)
An alternative expression for these curves can be obtained by expanding the sinusoidal terms in eq. 18 so that

$$r=(X_a\mathrm{cos}\theta ^{}+Y_a\mathrm{sin}\theta ^{})+\left(\frac{X_b}{\mathrm{cos}\theta ^{}}\frac{Y_b}{\mathrm{sin}\theta ^{}}\right),$$
(19)
which is now linear in four new parameters $`(X_a,Y_a,X_b,Y_b)`$; these parameters are related to the lens shape parameters $`(v,u)`$ and the source position $`(x_s^{},y_s^{})`$ by

$$X_a\left(1+vu\right)x_s^{},Y_a\left(1u+v\right)y_s^{},X_bux_s^{},Y_bvy_s^{}.$$
(20)
These curves (cf. eq. 18) have the nice property that they go through all image positions independent of the radial profile of the lensing galaxy. An example is the semi-hyperbolic curve in Fig. 1. The curve is determined by the source parameters $`(r_s,\theta _s^{})`$, the lens shape parameters $`(v,u)`$ and the lens position angle $`\theta _p`$. The radial profile can take any general physical profile, isothermal or power-law.

The boxiness parameter $`\beta `$ is such that the shape function $`f(\theta )`$ reduces to the usual elliptical form when $`\beta =\frac{1}{2}`$ (Witt & Mao 1997). In the case that $`\psi (\omega )`$ is linear in $`\omega `$, the models reduce to the simple models of Kassiola & Kovner (1995) when $`\beta =1`$. Interestingly, elliptical models with $`\beta =\frac{1}{2}`$ have $`uv=\frac{\delta }{\mathrm{\Delta }}=\frac{1}{2\beta }=1`$ (cf. eq. 17), hence $`X_a=Y_a=0`$, and eq. 19 reduces to

$$1=\frac{X_b}{x^{}}\frac{Y_b}{y^{}}$$
(21)
after applying $`x^{}=r\mathrm{cos}\theta ^{},y^{}=r\mathrm{sin}\theta ^{}`$. This equation prescribes a hyperbolic curve, which is consistent with Witt (1996) and Witt & Mao’s (1997) finding that all four image positions and the source position lie on a certain hyperbolic curve. A hyperbolic curve has a maximum of 5 free parameters, thus they cannot fit image positions of a general quadruple system; four image positions yield a minimal of eight constraints. Experience with fitting several quadruple systems (G2237+0305, CLASS1608+656, HST12531-2914) shows that models with $`|\beta |1/2`$ often give unphysical mass-radius relations $`m(\omega )`$. We find that $`\beta =\frac{1}{8}`$ gives a fair approximation to realistic models. These models have non-elliptical contours, and often yield physical density distributions. Fig. 2 shows that they also cover a sufficiently wide range of axis ratios for the potential and the density so that we can explore the shape of the lens galaxy in fitting the image positions. The expressions for the axis ratios of a power-law lens are given in Appendix A.

## 3 Results

### 3.1 Lens shape directly from fitting image positions

Our models can be used to fit image positions and derive lens shape parameters ($`\delta ,\beta `$) and the source position $`(r_s,\theta _s)`$, free from assumptions of the lens radial profile, but subject to the assumption that the lens’ angular profile obeyes eq. 12. The four unknowns can be derived from the four observed image positions, i.e., the *eight* observables $`(r_i,\theta _i)`$ with $`i=1,2,3,4`$. The position angle $`\theta _p`$ of the lens principal axis is treated as a free variable.

The procedure is simple. First substitute the four observed image positions in eq. 19 to obtain the following four linear equations

$$\frac{\mathrm{cos}\theta _i^{}}{r_i}X_a+\frac{\mathrm{sin}\theta _i^{}}{r_i}Y_a+\frac{1}{r_i\mathrm{cos}\theta _i^{}}X_b\frac{1}{r_i\mathrm{sin}\theta _i^{}}Y_b=1,i=1,2,3,4,$$
(22)
of the four unknown parameters $`(X_a,Y_a,X_b,Y_b)`$. After solving these, either analytically or numerically, these parameters are substituted in eqs. 2015, and 17 to yield the source position and the lens shape in terms of $`(r_s,\theta _s,\delta ,\beta )`$. In fact, we can recast eq. 20 to a set of four simple linear equations of four new unknowns $`(u,v,1/x_s^{},1/y_s^{})`$ by moving the terms $`x_s^{}`$ and $`y_s^{}`$ to the left hand side of the equations, i.e.,

$$X_a/x_s^{}v+u=1,Y_a/y_s^{}v+u=1,X_b/x_s^{}u=0,Y_b/y_s^{}v=0.$$
(23)
The lens shape parameters $`(\delta ,\beta )`$ can then be computed from

$$\delta =\frac{uv}{u+v},\beta =\frac{1}{2(uv)},$$
(24)
and the source position $`(r_s,\theta _s)`$ from

$$r_s=\sqrt{x_{s}^{}{}_{}{}^{2}+y_{s}^{}{}_{}{}^{2}},\theta _s=\mathrm{arctan}(y_s^{},x_s^{})+\theta _p.$$
(25)
Thus we have effectively reduced the problem of fitting image positions to successively solving linear equations, which is a straightforward task.

### 3.2 Non-parametric radial profile and power-law slope

The radial part of the potential can also be extracted from eq. 8 without parameterization. At the positions of the images we have

$$m(\omega )=r^2rr_s\mathrm{cos}(\theta \theta _s)=x^2+y^2xx_syy_s,\omega =rf(\theta ),$$
(26)
where we have used $`_r\mathrm{ln}\omega =r^1`$. Thus we have obtained a mass-radius relation directly from the observed image positions, assuming that the source position $`(x_s,y_s)`$ and the flattening and position angle of the potential $`(\delta ,\theta _p)`$ have been determined by fitting the curve (cf. eq. 18). Interestingly, in the limit that the source is at the center of a circular lens, we have $`r_s0`$, $`\beta 0`$ and $`mr^2`$.

It is useful to characterize the radial profile of a lens galaxy, which is generally not a power-law, by an effective power-law slope $`\alpha (\omega )`$, which varies with the radius $`\omega `$ except for scale-free models. There are several ways of estimating the characteristic power-law slope. Taking any two images $`i`$ and $`j`$, we can form a characteristic power-law slope $`\alpha _{ij}`$ from the mass $`m_i`$ and $`m_j`$ at the image radii $`\omega _i`$ and $`\omega _j`$ with

$$\alpha _{ij}\frac{\mathrm{log}m_i\mathrm{log}m_j}{\mathrm{log}\omega _i\mathrm{log}\omega _j}=\frac{2\mathrm{log}\frac{r_i}{r_j}+\mathrm{log}\frac{1r_sr_i^1\mathrm{cos}(\theta _i\theta _s)}{1r_sr_j^1\mathrm{cos}(\theta _i\theta _s)}}{\mathrm{log}\frac{r_i}{r_j}+\beta \mathrm{log}\frac{1\delta \mathrm{cos}(\theta _i\theta _s)}{1\delta \mathrm{cos}(\theta _j\theta _s)}},$$
(27)
where we have used the mass-radius eq. 26.

Alternatively we can estimate the power-law slope from the observed time delay. First we rewrite the time delay eq. 6, so that the time delay, $`t_{ij}`$, between any two images $`i`$ and $`j`$ is given by

$$\frac{ht_{ij}}{\tau _{100}}=\frac{1}{2}(x_i^2+y_i^2)\frac{1}{2}(x_j^2+y_j^2)(x_ix_j)x_s(y_iy_j)y_s(\psi _i\psi _j),$$
(28)
where $`\psi _i\psi _j`$ is the difference in the lens potential between the two images. Rewriting eq. 28 we form a new estimator $`\alpha _{ij}^t`$ for the power-law slope with

$$\alpha _{ij}^t\frac{m_im_j}{\psi _i\psi _j}=\frac{(x_i^2+y_i^2)(x_j^2+y_j^2)(x_ix_j)x_s(y_iy_j)y_s}{\frac{1}{2}(x_i^2+y_i^2)\frac{1}{2}(x_j^2+y_j^2)(x_ix_j)x_s(y_iy_j)y_s\frac{ht_{ij}}{\tau _{100}}},$$
(29)
where we have used eq. 26.

Thus we have two direct estimators $`\alpha _{ij}`$ and $`\alpha _{ij}^t`$ of the radial profile, computed from the observed images and time delays. In the limit of scale-free power-law models $`\alpha _{ij}=\alpha _{ij}^t=\alpha (\omega )=cst`$. So the deviation from scale-freeness can be estimated by taking the differences such as $`\alpha _{12}\alpha _{34}`$, $`\alpha _{12}^t\alpha _{34}^t`$, or $`\alpha _{14}\alpha _{14}^t`$.

### 3.3 Hubble constant and time delay ratios

To apply the above estimates of the power-law slopes, we have assumed that we know the rescaled Hubble constant $`h`$ from independent observations. Alternatively the Hubble constant $`H_0=100h`$ can be estimated from the time delay $`t_{ij}`$ between two images. Letting $`\alpha _{ij}^t=\alpha _{ij}`$, we get

$$h=\frac{\tau _{100}}{t_{ij}}\left\{\left(\frac{1}{2}\frac{1}{\alpha _{ij}}\right)\left[(x_i^2+y_i^2)(x_j^2+y_j^2)\right]\left(1\frac{1}{\alpha _{ij}}\right)\left[(x_ix_j)x_s+(y_iy_j)y_s\right]\right\}.$$
(30)
The time delay ratio can also be predicted with

$$\frac{t_{ij}}{t_{kl}}=\frac{(\frac{1}{2}\frac{1}{\alpha _{ij}})[(x_i^2+y_i^2)(x_j^2+y_j^2)](1\frac{1}{\alpha _{ij}})[(x_ix_j)x_s+(y_iy_j)y_s],}{\left(\frac{1}{2}\frac{1}{\alpha _{kl}}\right)\left[(x_k^2+y_k^2)(x_l^2+y_l^2)\right]\left(1\frac{1}{\alpha _{kl}}\right)\left[(x_kx_l)x_s+(y_ky_l)y_s\right]},$$
(31)
which is obtained by rewriting eq. 30.

## 4 Application: the surface density and time delay models for PG1115+080

### 4.1 Data

As a simple application, we model the image positions and time delays of the well-studied quadruple system PG1115+080. This system has been extensively studied ever since the first models by Young and collaborators (1981), and has received closer attention after Schechter et al.’s (1997) measurements of its time delay. All models, except those of Saha & Williams (1997), adopt elliptical/circular shapes and a few common radial profiles, with the models of Keeton & Kochanek (1997) and Impey et al. (1998) being the most comprehensive. Only one year after the discovery of the legendary double-image radio-loud quasar Q0957+561, this system was identified as a multiple-imaged system by Weymann et al. (1980) in their survey of nearby bright QSOs. It is now known to consist of four images with the names $`A_1`$, $`A_2`$, $`B`$ and $`C`$ (with flux ratios about 4 : 2.5 : 0.7 : 1, cf. the HST observations of Kristian et al. 1993) of a radio-quiet QSO at redshift $`z_s=1.722`$. The images $`A_1`$ and $`A_2`$ are within $`0.5^{\prime \prime }`$; see the inset of Fig. 1. Interestingly, the lens galaxy is also one of the bright members of a galaxy group ($`N10`$) at redshift $`z_l=0.310`$, first mapped by Young et al. (1981). The center of the group is to the south-west of the lens, roughly at $`r_g=(20^{\prime \prime }\pm 2^{\prime \prime })`$ and $`\theta _g=(117^o\pm 3^o)`$. The lens galaxy has been resolved by both HST and the 8.2-m Subaru telescope in $`0.3^{\prime \prime }`$ seeing. It appears to be an early type galaxy with a de Vaucouleurs profile and a half-light radius of $`0.55^{\prime \prime }`$. There is no sign of differential dust-extinction in the lens galaxy. While NICMOS observations by Impey et al. (1998) show no flattening for the lens, ground infrared images by Iwamuro et al. (2000) found it to be an E1 galaxy elongated towards $`\theta _p65^o`$. Both observations reveal a $`1^{\prime \prime }`$ infrared Einstein ring connecting the four images, which is thought to be the infrared image of the QSO host galaxy. PG1115+080 is also one of the two quadruple systems where the time delay between images has been measured, the other one being the radio-loud quasar B1608+656 from the CLASS survey (cf. Fassnacht et al. 1999). Although two different sets of values are quoted in the literature (Schechter et al. 1997, Barkana 1997), the leading image is the furtherest image (the image C), and the innermost image (image B) arrives last. The time delay ratio $`r_{ABC}=t_{AC}/t_{BA}`$, for the delay between image $`C`$ and $`A_1+A_2`$ vs. image $`B`$ and $`A_1+A_2`$, provides an extra discriminator of the models; the images $`A_1`$ and $`A_2`$ are within $`0.5^{\prime \prime }`$ of each other, and the small relative delay is undetected. Schechter et al. first reported $`r_{ABC}=0.7\pm 0.3`$ from their photometric monitoring program in 1995-1996. Later analysis by Barkana (1997) found $`r_{ABC}=1.13\pm 0.18`$, after taking into account correlations of errors amongst the time delays. The delay between images B and C, $`t_{BC}=25.0\pm 1.7`$ days.

Here we illustrate the application of our models to the most recent data from Impey et al. (1998) of PG1115+080. We do not attempt to model B1608+656 because the lens appears to be a merging pair of galaxies, and the morphology is too complex for our model. We denote with the index $`i=1,2,3,4`$ the four images $`A_1`$, $`A_2`$, $`B`$ and $`C`$. All results are quoted for a standard flat universe without a cosmological constant $`(\mathrm{\Omega },\mathrm{\Lambda })=(1,0)`$. Table 1 gives the relevant quantities to calibrate our results to other universes. The angular distance from redshift $`z_1`$ to $`z_2`$ in a universe of a matter and vacuum density $`(\mathrm{\Omega },\mathrm{\Lambda })`$ times the closure density is generally given by

$$D(z_1,z_2)=\frac{c}{H_0(1+z_2)\sqrt{\mathrm{\Omega }_c}}\mathrm{sinh}\left\{\sqrt{\mathrm{\Omega }_c}_{z_1}^{z_2}𝑑z\left[\mathrm{\Lambda }+\mathrm{\Omega }_k(1+z)^2+(1+z)^3\mathrm{\Omega }\right]^{\frac{1}{2}}\right\},\mathrm{\Omega }_k=1\mathrm{\Lambda }\mathrm{\Omega }.$$
(32)
The predicted Hubble constant should be reduced from the value achieved with the standard $`(\mathrm{\Omega },\mathrm{\Lambda })=(1,0)`$ universe by 3% in the presently favored $`\mathrm{\Lambda }`$-dominated universe with $`(\mathrm{\Omega },\mathrm{\Lambda })=(0.3,0.7)`$ from surveys of distant supernovae.

### 4.2 The source position, the lens shape and mass inside images

First we solve for the lens shape and source position from the linear equations 19 and 23. The solutions for PG1115+080 are given in Table 2, sorted according to the value of $`\beta `$ or $`\theta _p`$. The resulting potential model with $`\beta =\frac{1}{8}`$ or $`\theta _p=60.5^{}`$ has a flattening of between E0 and E1, and interestingly the lens principal axis points towards the location of the galaxy group in the lens plane. Models with other values of $`\beta `$ or $`\theta _p`$ will be discussed in section 5.

To proceed with determining the lens mass at each image position, we substitute the now known flattening parameters $`(\delta ,\theta _p)`$ and the source positions $`(r_s,\theta _s)`$ in eq. 26, to predict four independent data points $`(\omega _i,m_i)`$ with $`i=1,2,3,4`$ in the radius vs. the enclosed mass plane. Fig. 3 shows the predicted lens mass enclosed at the four image positions. Note that the mass rises faster than the light, implying a growing dark mass component at large radius; the light distribution is modeled as an observed de Vaucouleurs $`r^{\frac{1}{4}}`$-law with a half-light radius of $`0.55^{\prime \prime }`$.

### 4.3 Piecewise power-law model

So far we did not enforce any strict parameterization of the radial profile. We only restrict the profile to be of the form of eq. 11; in practice, this is an insignificant restriction. In the following sections we show several ways of modeling the radial profile assuming an isolated lens. None of the models is completely satisfactory.

First we use a minimal model, which assumes that the mass-radius relation is a piecewise power-law, that is, we connect a straight line through two images in the $`\mathrm{log}m`$ vs. $`\mathrm{log}\omega `$ plane. The piecewise values for the power-law slope and axis ratios are given in Table 2.

The Hubble constant $`H_0=100h`$ can be estimated by normalizing the time delay to the observed value $`t_{AC}`$ (Schechter et al. 1997),

$$h=\frac{\tau _{100}}{2t_{AC}}\left[P_{14}+\left(1\frac{1}{\alpha _{14}}\right)S_{14}\right],$$
(33)
where

$$P_{ij}=\left[(x_i^2+y_i^2)(x_j^2+y_j^2)\right]$$
(34)
depends on the image positions $`(x_i,y_i)`$ alone and

$$S_{ij}=2P_{ij}2\left[(x_ix_j)x_s+(y_iy_j)y_s\right]$$
(35)
depends on the source position $`(x_s,y_s)`$ as well. The model yields a Hubble constant $`H_030`$ km/s/Mpc, much lower than determined by other authors (e.g., Impey et al. 1998). The low $`H_0`$ is a result of the high value for the power-law slope $`\alpha _{AC}=1.6`$. Other values are given in Table 2 for various image pairs and observed time delay. The time delay ratio can also be predicted with (cf. eq. 31)

$$r_{ABC}\frac{t_{14}}{t_{32}}=\frac{P_{14}+\left(1\frac{1}{\alpha _{14}}\right)S_{14}}{P_{32}+\left(1\frac{1}{\alpha _{32}}\right)S_{32}}.$$
(36)
Substituting the slopes $`\alpha _{AC}`$ and $`\alpha _{BA}`$ to eq. 31 we find the ratio between images $`r_{ABC}=t_{AC}/t_{BA}0.65`$.

Note that the time delay predictions here are robust and independent of details of the density, since the discontinuity in the density is completely smoothed out in the lensing potential.

### 4.4 Single power-law model

The piecewise-power-law model above necessarily implies a discontinuous density profile. This could be cured if we enforce a single power-law model, that is, we fit a straight line to the four points in the $`\mathrm{log}m`$-$`\mathrm{log}\omega `$ plane. This would give us a power-law slope $`\alpha =1.38`$. Fig. 1 shows our model surface density contours. Similar to the non-parametric models of Saha & Williams (1997) we find a peanut-shaped lens.

We can estimate the goodness of the fit by recomputing the image positions from a given mass model. The general procedure of simulating images of our theoretical lens model is as follows: First combine eq. 26 with the power-law radial profile to get

$$r^2rr_s\mathrm{cos}(\theta \theta _s)=m(\omega )=b_0r_0\left(\frac{\omega }{r_0}\right)^\alpha ,\omega =r\left|1\delta \mathrm{cos}(2\theta 2\theta _p)\right|^\beta .$$
(37)
Then upon substitution of eq. 18 to eliminate $`r`$, we obtain a one-dimensional non-linear equation for the image position angle $`\theta `$, which can be solved easily numerically. As a comparison, one would be dealing with a minimum-finding or a root-finding numerical problem in a two-dimensional plane $`(r,\theta )`$ in the general case without the separation of the angular vs. radial part. For our model of PG1115+080, the image solutions are shown in Fig. 1 as dashed circles, together with the input observed image positions (solid circles). The predicted four images are off by about 60 milli arcsecs, a residual which is inconsistent with the $`3`$ milli arcsecs accuracy of Impey et al. positions, and is marginally consistent with earlier data by Kristian et al. with a quoted error of 50 milli arcsec for the lens galaxy and 5 milli arcsec for the images. The images $`A_1`$ and $`A_2`$ are at two sides of the critical curve (the line of infinite amplification in the source plane), hence are highly amplified with opposite parity. The time delay ratio can be estimated with eq. 31, assuming a constant power-law slope $`\alpha =\alpha _{BC}=1.38`$. This yields $`r_{ABC}1.3`$, close to the Barkana (1997) value of $`1.13\pm 0.18`$. The large difference here is due to the large residual in terms of fitting the images $`A_1`$ and $`A_2`$ with a straight power-law.

We also compute the amplification patterns by taking double derivatives of the time delay surface. The circles in Fig. 1 show the observed (solid circles) and predicted (dashed circles) fluxes of each image and the source (half-closed circle), with the area of each circle in proportion to the flux. The $`B`$ to $`C`$ ratio is well-reproduced and the predictions for images $`A_1`$ and $`A_2`$ are also consistent with observations at about the 0.3 magnitude level.

Finally, if we put a host galaxy around the QSO, the model predicts that the image of the host galaxy will be stretched into an arc. We see a nearly closed ring. This agrees very well with the diffuse ring that Impey et al. discovered in their NICMOS images. The ring maps back to the source plane as a disk with an area of $`0.03`$ square arcsec.

### 4.5 Double power-law model

Alternatively we can fit a smooth, five-parameter lens model with a lens potential

$$\psi (\omega )=c_0\left(\frac{\omega }{a_0}\right)^{\alpha _{in}}\left[1+\left(\frac{\omega }{a_0}\right)^n\right]^{\frac{\alpha _{out}\alpha _{in}}{n}},\omega =r\left|1\delta \mathrm{cos}(2\theta 2\theta _p)\right|^\beta .$$
(38)
The corresponding mass profile is given in Appendix B. This model assumes the lens potential (as well as the lens mass) increases like a double-power-law with an inner slope $`\alpha _{in}`$ and outer slope $`\alpha _{out}`$. The transition is defined by the normalization constant $`c_0`$, the radius $`r=a_0`$ and the sharpness parameter $`n`$; bigger $`n`$ corresponds to sharper transition. When fitting the mass model to the four points in Fig. 3, it turns out that the inner slope $`\alpha _{in}`$ is fully unconstrained. The other four free parameters $`(c_0,a_0,\alpha _{out},n)`$ are determined from fitting the four data points; the procedure is explained in Appendix B. The values of $`(c_0,a_0,\alpha _{out})(0.6,1.1,1.63)`$, nearly independent of the value for the inner slope $`\alpha _{in}`$. Note $`\alpha _{out}\alpha _{AC}`$, consistent with the fact that the mass at $`A_1`$, $`A_2`$ and $`C`$ follows nearly a power-law (cf. Fig. 3). The value of $`n`$ increases from $`15`$ to $`39`$ for $`\alpha _{in}=02`$. That $`n2`$ means that the density changes sharply at the transition. All fits have zero residual and predict nearly identical mass profiles at radii between images $`B`$ and $`C`$. They differ only in the mass profile inside image $`B`$, where we have no direct constraint on the dark matter profile. For all purposes it is sufficient to set $`n=20`$ so that $`\alpha _{in}2`$, in which case the model has a finite core at small radii. From the potential model we can compute the dimensionless surface density contours (cf. Fig. 4). Near the position of the images, the density has a flattening of E2-E3, flatter than the potential, as expected. The potential model can also be substituted in eq. 6 to predict the time delay contours, shown in Fig. 4. The observed images $`C`$, $`A_2`$, $`A_1`$, $`B`$ (the solid circles) are exactly the valley, peak, saddle, valley points of the delay contour, where the theoretical images should lie according to Fermat’s principle. The images $`A_1`$ and $`A_2`$ have nearly the same arrival time. Substituting the lens potential in eq. 28 we can predict the time delay ratio for the double-power-law model. The result, $`r_{ABC}0.65`$, is nearly the same as the piece-wise power-law model. These predicted ratios are in good agreement with the earlier measurement of $`0.7\pm 0.3`$ of Schechter et al. (1997).

We can recompute the image positions from the double-power-law model by solving

$$r^2rr_s\mathrm{cos}(\theta \theta _s)=m(\omega )=c_0\left(\frac{\omega }{a_0}\right)^{\alpha _{in}}\left[1+\left(\frac{\omega }{a_0}\right)^n\right]^{\frac{\alpha _{out}\alpha _{in}}{n}},\omega =r\left|1\delta \mathrm{cos}(2\theta 2\theta _p)\right|^\beta ,$$
(39)
and eq. 18.

Unlike the single-power-law model it predicts five images (cf. Fig. 4). Four predicted images fall exactly on the observed positions. But the model predicts an extra image near image B. This can be identified with the usually highly demagnified fifth image that is sometimes observed in lensed systems. The extra image the model predicts is, however, too bright to be consistent with observations. This is a direct consequence of the fact that $`\kappa `$ at its radius is nearly constant and $`1`$. The amplification can be calculated from $`\kappa `$ as $`(1\kappa )^225`$.

All images lie on the curve defined by eq. 18. The extra image arises because the model density profile (cf. Fig. 5) is non-monotonic near $`1^{\prime \prime }`$ radius, the transition radius $`a_0`$ of the double-power-law potential model. A wiggle in density happens when the transition parameter $`n>2`$, i.e., a sharp transition of the potential.

## 5 Discussion

### 5.1 Power-law slope of PG1115+080 vs. $`H_0`$ and time delay ratio

The previous sections have illustrated the procedure for fitting known quadruple image systems with our models, and have also revealed some puzzling problems with PG1115+080, for example, the low and non-unique $`H_0`$ and the peculiar and non-unique lens density. Fig. 6 shows the values for $`\alpha `$ predicted from either the mass-radius relation (cf. eq. 27) or from the time delay (cf. eq. 29). Two comments are in order: (1) There is a large spread with the power-law slope $`\alpha `$ depending on how it is predicted. To fit the observed time delay within $`2\sigma `$ with a reasonable $`H_0`$ (between 50 and 100 km/s/Mpc), we obtain only a loose constraint $`0.3<\alpha <1.6`$. (2) The rise of the power-law slope $`\alpha `$ with radius in the piecewise-power-law model is somewhat unphysical; realistic models typically have a monotonically steeper density profile (smaller $`\alpha `$) with increasing radius. These results are consistent with the finding of a wiggle near the transition radius ($`1^{\prime \prime }`$) in the surface density of the smooth double-power-law fits. Nevertheless the density is positive everywhere.

A detailed treatment of these problems should include the effect of the neighboring group. Such models are presented elsewhere, since they are beyond the interest of this paper on methods for quadruple systems in general. Nevertheless some further study of the model parameter space is clearly needed. Here it is sufficient to comment on the range of the predicted model parameters, in particular, the power-law slope and the value of $`H_0`$ when we vary the parameter $`\beta `$ or $`\theta _p`$.

Fig. 7 shows the parameter space of the isolated lens model by varying the orientation of the lens principal axis $`\theta _p`$ (cf. eq. 12 and 13). It turns out a one-parameter family of models fits the observations; we effectively vary $`\beta `$. In fact, $`\theta _p61^o+4\beta `$, and $`\mathrm{\Delta }=2\beta \delta \frac{\beta }{2}`$ to a fair accuracy. We find that models with $`\beta <1/3`$ will generally result in negative density models, models with $`\beta 0`$ have a radially-increasing density, and models with $`|\beta |0`$ also produce a nearly constant axisymmetric density, hence implying an unphysically small Hubble constant. Only in the range $`\frac{1}{4}\beta <0`$ do we find physical models with a reasonable flattening. The source positions for different models are shown in Fig. 8. Table 2 compares the predicted parameters of the $`\beta =\frac{1}{4}`$ or $`\theta _p=60^{}`$ model with the $`\beta =\frac{1}{8}`$ or $`\theta _p=60.5^{}`$ model. Fig. 3 shows the predicted radial profile for both models. Both models would have a significant residual if we fit the image positions with a straight power-law, or strange wiggles and extra images if we fit a double-power-law. We find this is generally the case with all isolated lens models without the shear from the neighbouring group.

Fig. 9 shows that there is significant degeneracy with the value for $`H_0`$, depending on the adopted model and the adopted time delay. The most reliable estimate is from the delay between the innermost image B and the outermost image C using

$$h=\frac{\tau _{100}}{2t_{BC}}\left[P_{14}+\left(1\frac{1}{\alpha _{14}}\right)S_{14}P_{32}+\left(1\frac{1}{\alpha _{32}}\right)S_{32}\right],$$
(40)
where we neglected the time delay between the images $`A_1`$ and $`A_2`$. The $`H_0`$ thus predicted scales approximately as $`60(2\alpha )`$ km/s/Mpc (cf. Fig. 9), where $`\alpha =\alpha _{BC}`$ is the average power-law slope. The factor $`2\alpha `$ is the average power-law slope of the model surface density, and is not well-constrained at the radius of the images. Models with $`\beta \frac{1}{4}`$ or $`\theta _p=60^{}`$ predict a steeper density profile, and closer to isothermal than models with $`\beta \frac{1}{8}`$ or $`\theta _p=60.5^{}`$, hence yield a more plausible $`H_0`$ around $`60`$ km/s/Mpc. Models with very shallow profiles ($`\alpha >1`$) are unfavored by the consensus value for $`H_0`$.

In all our isolated lens models, the delay ratio is closer to the Schechter value, while the Barkana value seems to be the widely accepted one. This explains why we see in Fig. 9 a tighter prediction of $`H_0`$ using $`t_{AC}`$ and $`t_{BA}`$ of Schechter values than the Barkana values. Interestingly in the limit that the isothermal model is applicable, $`\alpha =1`$, the time delay ratio can be estimated from the image positions directly (i.e. without computing the source position) with the following equation:

$$r_{ABC}=\frac{P_{14}}{P_{32}}=\frac{r_1^2r_4^2}{r_3^2r_2^2}0.65,\text{for isothermal lens},$$
(41)
where $`r_i`$ is the observed radius of the image $`i`$ from the lens center. The ratio increases when we introduce the shear from the nearby group.

### 5.2 Rotation curve and mass-to-light ratio

Fig. 10 compares the predicted velocity dispersion for the lens with the observed value. The dispersion is predicted with the following formula,

$$\sigma _i^2=V^2\frac{m_i}{2r_i},V^2\frac{D_{os}}{D_{ls}}\frac{(3\times 10^5\mathrm{km}/\mathrm{s})^2}{2\pi \mathrm{radian}}\frac{1\mathrm{radian}}{206265\mathrm{arcsec}},$$
(42)
where the predictions are made from each image $`i=1,2,3,4`$ (strictly speaking the formula is valid only for a singular isothermal lens). Applying this for each value of the lens position angle $`\theta _p`$ we get the full range of possible values for the dispersion. The difference of the dispersion among the four images shows the deviation from an isothermal model. On average the predicted lens galaxy velocity dispersion, 200-300 km/s, suggests that the lens is close to a massive $`L_{}`$ galaxy. The observed dispersion $`\sigma =(281\pm 50)`$ km/s ($`95\%`$ confidence, shaded region) from the Keck spectrum of Tonry (1998). This value is comparable to the theoretical models, but on the high side. A similar problem has been noted in previous models (Schechter et al. 1997). Tonry noted that the discrepancy might be due to a steep radial fall-off of the velocity dispersion or a simple over-estimation of the observed dispersion; the spectrum of the faint lens was made from subtracting two Keck spectra taken with a $`1^{\prime \prime }`$-slit, one passing the lens and the image B, one passing the images $`A_1`$ and $`A_2`$.

Fig. 11 compares the increase of the lens mass from the innermost image $`B`$ to the outermost image $`C`$ vs the increase of the lens light. The mass ratio is predicted using

$$\frac{M_C}{M_B}=\frac{m_4}{m_3}=\left(\frac{r_4}{r_3}\right)^{\alpha _{BC}}.$$
(43)
where the power-law slope $`\alpha _{BC}`$ depends the principal axis $`\theta _p`$. The light ratio

$$\frac{L_C}{L_B}=\frac{_0^{r_4}𝑑r2\pi r\mathrm{exp}\left[7.67\left(\frac{r}{r_e}\right)^{\frac{1}{4}}\right]}{_0^{r_3}𝑑r2\pi r\mathrm{exp}\left[7.67\left(\frac{r}{r_e}\right)^{\frac{1}{4}}\right]},r_e=0.55^{\prime \prime },$$
(44)
where the observed de Vaucouleurs law is used for the projected light. The fact that $`\frac{M_C}{M_B}>\frac{L_C}{L_B}>1`$ implies that the mass grows faster than the light at these radii for all these models, consistent with a large amount of dark matter between 0.8 and 1.5 arcsec.

## 6 Conclusion

Here we itemize our main results.

(i) We found a very general class of lens models that allow for non-elliptical and non-scale-free lenses. These models include previous isothermal elliptical models as special case. We can derive the radial mass distribution of the lens in a non-parametric way. We can study the deviation from the usually assumed straight power-law profile. We also give simple formulas for computing the time delay ratios and for estimating $`H_0`$.

(ii) The models are very easy to compute and can be used for efficient exploration of a large parameter space because the lens equations can be reduced to a set of linear equations.

(iii) We have applied the models to PG1115+080, and have explored a large parameter space of isolated lens models. All models using piece-wise-power-law or double-power-law fit the positions of the images exactly. The models can also reproduce the flux ratios between the images and the stellar velocity dispersion of the lens approximately.

(iv) Our models are consistent with a dark halo up to a radius of three times the half-light-radius of the lens. The enclosed mass increases much faster than the enclosed light as we move radially from the innermost image to the outermost image.

(v) We reconfirm earlier results by (e.g. Schechter et al. 1997) that the principal axis of the lens potential points, within a few degrees, to the external group, and is consistent with the observed value (Iwamuro et al. 2000).

(vi) Our models do not yield a unique prediction of $`H_0`$, because it is sensitive to the lens power-law slope $`\alpha `$ (cf. Fig. 9), which appears to have a large spread, depending on how it is predicted. The power-law slope is also sensitive to small uncertainties of the position angle of the lens $`\theta _p`$. A change of $`\theta _p`$ by a few degrees from the observed value can change $`\alpha `$ from $`0`$ to $`2`$ and $`H_0`$ from 100 km/s/Mpc to 0.

(vii) Our models predict consistently a low time delay ratio $`r_{ABC}`$, which fits only the Schechter value. This and the peculiar oscillation in the predicted density profiles, we believe, are due to the neglect of the external group. This will be dealt with in detail in a follow-up paper.

In summary, our semi-analytical lens models allow us to explore a large range of physical lens density distributions. We find that there is a large systematic uncertainty in the lens models from fitting image positions and time delays of PG1115+080 and isolated lens models are always unsatisfactory for this quadruple system.

The authors thank Tim de Zeeuw for a careful reading of the manuscript and the referee for a constructive report. HSZ thanks Paul Schechter for discussions and encouragement.

## Appendix A Power-law models

Within a certain range of radii, realistic lens mass distributions can be approximated by the following bi-symmetric power-law of the radius $`\omega `$

$$m(\omega )\omega \frac{d\psi }{d\omega }=b_0r_0\left(\frac{\omega }{r_0}\right)^\alpha ,\omega =rf(\theta )=r\left|1\delta \mathrm{cos}(2\theta 2\theta _p)\right|^\beta ,r_0=1^{\prime \prime },0<\alpha <2,$$
(45)
where $`b_0`$ is the characteristic deflection strength at the characteristic length scale $`r_0`$ (one arcsec) of the lens system, and $`\theta _p`$ defines the principle axis of the lens. The corresponding lens potential is related to the enclosed mass $`m(\omega )`$ by

$$\psi (\omega )=\mathrm{const}+\frac{m(\omega )}{\alpha }.$$
(46)
The axis ratio (or its inverse) of the bi-symmetric potential

$$q_\psi \left(\frac{1\delta }{1+\delta }\right)^\beta 1\mathrm{\Delta },\mathrm{\Delta }=2\beta \delta ,$$
(47)
and the axis ratio of the density

$$q_\kappa \left(\frac{1\delta }{1+\delta }\right)^{\frac{1\alpha \beta }{2\alpha }}\left[\frac{1(\frac{4\beta }{\alpha }1)\delta }{1+(\frac{4\beta }{\alpha }1)\delta }\right]^{\frac{1}{2\alpha }}1\mathrm{\Delta }(1+2\alpha ^1),$$
(48)
where the approximations are valid if the flattening is small. The surface density $`\kappa `$ of the model is given by

$$\kappa =\frac{1}{2}br^{\alpha 2}g^{\alpha \beta }(\alpha g^2+(\beta ^2\alpha \beta )g^2+\beta g^{\prime \prime }g),$$
(49)
where $`g1\delta \mathrm{cos}(2\theta 2\theta _p),g^{}2\delta \mathrm{sin}(2\theta 2\theta _p),g^{\prime \prime }4\delta ^2\mathrm{cos}(2\theta 2\theta _p)`$.

## Appendix B Double-power-law models

For the double power-law lens potential

$$\psi (\omega )=c_0\left(\frac{\omega }{a_0}\right)^{\alpha _{in}}\left[1+\left(\frac{\omega }{a_0}\right)^n\right]^{\frac{\alpha _{out}\alpha _{in}}{n}},$$
(50)
the corresponding mass profile is given by

$$m(\omega )=\omega \frac{d\psi (\omega )}{d\omega }=\alpha (\omega )\psi (\omega ),$$
(51)
where

$$\alpha (\omega )=\frac{d\mathrm{log}\psi (\omega )}{d\mathrm{log}\omega }=\frac{\alpha _{in}+\alpha _{out}\left(\frac{\omega }{a_0}\right)^n}{1+\left(\frac{\omega }{a_0}\right)^n}.$$
(52)
Generally speaking, a double-power-law fit (cf. eq. 38) involves searching for solutions of the following equations in a five-parameter space $`(n,\xi _0,\alpha _{out},\gamma _0,a_0)`$:

$$\mathrm{log}m_i=\left[\xi _0+(1+\gamma _0)\alpha _{out}\mathrm{log}\omega _i\right]\left(1+\frac{\gamma _0\alpha _{out}}{n}\right)\mathrm{log}\left[1+\left(\frac{\omega }{a_0}\right)^n\right]+\mathrm{log}\left[1+\gamma _0+\left(\frac{\omega }{a_0}\right)^n\right],i=1,2,3,4,$$
(53)
where

$$\gamma _0\frac{\alpha _{in}}{\alpha _{out}}1,\xi _0\mathrm{log}c_0+\mathrm{log}\alpha _{out}(1+\gamma _0)\alpha _{out}\mathrm{log}a_0,$$
(54)
and $`(\omega _i,m_i)`$ for $`i=1,2,3,4`$ are the four data points in the mass-radius diagram (cf. Fig. 3). Note that the equations 53 are linear in $`(\xi _0,\alpha _{out})`$, so these two variables can be eliminated by Gaussian substitution. If we fix $`n`$, say $`n=20`$, then we are left with two equations for two variables $`(\gamma _0,a_0)`$, which can be solved with a two-dimensional iterative root finding routine. The solution converges rapidly from an initial guess of $`\gamma _0=0`$ and $`a_0=1^{\prime \prime }`$.