# Pixelated Lenses and 𝐻₀ from Time-delay QSOs

## 1 Introduction

Most ways of measuring the Hubble constant involve a form of distance ladder, which utilizes a number of astrophysical standard candle and standard ruler relations, and is calibrated locally by a geometrical technique such as parallax (e.g., Madore et al. 1999, Madore et al. 1998, Kennicutt 1995). A recent exciting development in this field is to extend the reach of the geometrical rung of the distance ladder by using masers in orbit around galaxy centers to get distances to nearby galaxies thus bypassing Cepheids (Herrnstein et al. 1999). A few methods involve no distance ladder: good examples are (i) inferring the distance of Type II supernovae from their light curves and spectra by modeling their expanding photospheres (Schmidt et al. 1992), and (ii) comparing the $`H_0`$-independent angular extent of galaxy clusters to their $`H_0`$-dependent depth as deduced by the X-ray emission, and the Sunyaev-Zeldovich microwave background decrement due to the cluster (Hughes & Birkinshaw 1998).

But the most ‘one-step’ method of all was proposed by S. Refsdal in 1964, though it has only recently become feasible. The principle of Refsdal’s method is simple. In a system where a high-redshift QSO is split into multiple images by an intervening galaxy lens, the difference in light travel time between different images (observable as time delays if the QSO is variable) is proportional to the scale factor of the universe. The time delay is given by the schematic formula

$`\text{Time delay}`$ $`=`$ $`h^1\times \text{1 month}\times \text{image separation in arcsec}^2`$   
$`\times z_{\mathrm{lens}}\times \text{weak dependence on }z_{\mathrm{lens}}\text{,}z_{\mathrm{QSO}}\text{, and cosmology}`$   
$`\times \text{lens-mass-distribution dependent factor}`$   
where the last two factors are of order unity. To obtain $`H_0`$ using this method one requires three types of input: (i) the observed time delay(s) between QSO images, (ii) knowledge of the cosmology, and (iii) the mass distribution in the lensing galaxy. The first can and has been measured with increasing precision for about eight systems so far. The second is not a serious problem, because the dependence on cosmology is weak and the uncertainty due to it is easy to quantify; in this paper we will refer all results to the Einstein-de Sitter cosmology. The uncertainty in $`H_0`$ is dominated by the third item; the number of usable constraints on the mass distribution in the galaxy are few, while the range of possible distributions is huge. Thus, mass distribution is the major source of uncertainty.

Two different paths can be taken to compensate for our lack of knowledge about the galaxy. One is to assume an exact parametric form for the galaxy mass distribution and fit the observed lensing properties as best as possible; the other is to take the image properties as exact, and try to reconstruct the galaxy mass map as best as possible. Single parametric models fix the last term in (1) and thus cannot account for the uncertainty resulting from it. Blandford & Kundić (1996) advise that even if one finds a parametric galaxy model which is dynamically possible and which reproduces the image properties with acceptably low $`\chi ^2`$, one still has to ‘aggressively explore all other classes of models’ to get the true uncertainty in $`H_0`$. To explore the model space in a systematic fashion one needs to use a representation of the galaxy mass distribution that is general and not restricted to a particular form. One way would be to expand the mass distribution using a set of basis functions, another is to pixelate the galaxy and take each pixel as an independent mass element. We introduced pixelated models in Saha & Williams (1997, hereafter SW97) but at that time did not have any strategy for searching model space. We have now extended that work to explore the model space with the goal of estimating the uncertainty in the derived value of $`H_0`$.

The plan of this paper is as follows. In Section 2 we summarize the observational situation with regard to strongly lensed QSOs. In Section 3 we present the general lensing formalism and point out a few properties of the lensing equations that are useful in interpreting the results of modeling. We also explain the reasons for confining our analysis to PG1115+080 and B1608+656 for now. Sections 4 and 5 describe our method for deriving $`H_0`$ and test it on a synthetic sample via a blind test. Application to the real systems can be found in Section 6. Section 7 discusses our results.

## 2 Observed Time-Delay Lenses

The first piece of input for $`H_0`$ determination is the measurement of time delays between the various QSO images. At the present time, ten multiply-imaged QSOs already have measured time delays or are being monitored: Q0957+561 (Kundić et al. 1997a), PG1115+080 (Schechter et al. 1997, Barkana 1997), B1608+656 (Fassnacht et al. 1999), B0218+257 (Biggs et al. 1999), PKS 1830-211 (Lovell et al. 1998), HE 1104-1805 (Wisotzki et al. 1998), B1030+074, B1600+434 (Burud et al. 1999), J1933+503, and RXJ0911+0551 (Hjorth et al. 1999). In this work we limit ourselves to 4-image lenses with known source and lens redshifts and accurate time delay measurements; PG1115+080 and B1608+656 fit the description.

PG1115 (Weymann et al. 1980) was the second lens to be discovered. The source is a radio-quiet QSO at $`z_s=1.722`$. Accurate positions for the images were measured by Kristian et al. (1993); lightcurves were analyzed by Schechter et al. (1997), and time delays derived by Schechter et al. and Barkana (1997). The main lensing galaxy is an outlying member of a small galaxy group, $`z_l=0.311`$ with an estimated line of sight velocity dispersion of $`270\pm 70`$km s<sup>-1</sup> (Kundić et al. 1997b). A summary of observational results on this system can be found in SW97.

B1608 was discovered in the Cosmic Lens All-Sky Survey (Myers et al. 1995, Myers et al. 1999). The lens is either a perturbed single galaxy or a merging/interacting pair of galaxies superimposed in the plane of the sky. The source and lens redshifts are 1.394 and 0.630 respectively. The time delays were recently reported by Fassnacht et al. (1999) based on VLA observations spanning 7 months. The time delays we use in this work are an earlier determination (Fassnacht, private communication), and are less than 0.5$`\sigma `$ away from the values quoted in Fassnacht et al. (1999); $`\mathrm{\Delta }t_{BA}=28.5`$, $`\mathrm{\Delta }t_{BC}=32`$, and $`\mathrm{\Delta }t_{BD}=77`$.

## 3 Lensing formalism

A photon traveling through a galaxy will take longer to arrive at the observer then an unimpeded photon. Part of the time delay occurs because the path of the ray bundle makes a detour rather than going straight; the time delay is further increased because the photon travels through the gravitational potential well of the galaxy. The total time delay is given by,

$$\tau (\text{`},\text{`}_\mathrm{s})=(1+z_\mathrm{l})\frac{D_{\mathrm{ol}}D_{\mathrm{os}}}{D_{\mathrm{ls}}}\left[\frac{1}{2}(\text{`}\text{`}_\mathrm{s})^2\frac{1}{\pi }d^2\text{`}^{}\kappa (\text{`}^{})\mathrm{ln}|\text{`}\text{`}^{}|\right]$$
(2)   
where \` is the position on the sky, $`\text{`}_\mathrm{s}`$ is the source position, $`D`$’s are the angular diameter distances between the source, the lens and the observer, $`z_\mathrm{l}`$ is the redshift of the lens galaxy, and $`\kappa (\text{`})`$ is the projected mass density in the galaxy in units of $`\mathrm{\Sigma }_{\mathrm{crit}}=(c^2/4\pi G)(D_{\mathrm{os}}/D_{\mathrm{ls}}D_{\mathrm{ol}})`$.

If the lens mass distribution $`\kappa (\text{`})`$ is known then the arrival time surface, Eq. (2) provides us with all the necessary information about the images. Time delay between any two images is just the difference between $`\tau `$ at the relevant locations. According to Fermat’s Principle the images appear at stationary points of the arrival time surface,

$$\frac{\tau }{\text{`}}=0=\text{`}\text{`}_\mathrm{s}\frac{1}{\pi }d^2\text{`}^{}\kappa (\theta ^{})\frac{\text{`}\text{`}^{}}{|\text{`}\text{`}^{}|^2}$$
(3)   
Image distortion and magnification are given by the inverse of the curvature matrix of the arrival time surface

$$\left[\frac{^2\tau }{\theta _i\theta _j}\right]^1$$
(4)   
A few things can be learned by looking at the arrival time and lens equations:

(1) The time ordering of the images can be deduced from the image configuration using the morphological properties of the arrival time surface. The image furthest from the lensing galaxy is always first, and the one nearest the galaxy the last. In four-image QSOs the second image is the one opposite the first. Figure 1 illustrates.

(2) When four images are formed by an isolated galaxy of typical ellipticity the images are located nearly at the same galactocentric distance. This is easy to see by considering the two pieces of the arrival time surface. If the source and the center of the galaxy are not well aligned, i.e., if the ‘bump’ due to the gravitational potential contribution is away from the ‘well’ of the geometrical contribution, then the steepness of the geometrical part allows only two images to form, one roughly on either side of the galaxy center. To get four images, the bump of the gravitational contribution must be centered close to the source location. In such a situation the total arrival time configuration is centrally symmetric and the resulting images are approximately equidistant from the galaxy center.

(3) If the four images of a single source are located at different galactocentric distances the simplest explanation is the presence of external shear. External shear effectively raises the gravitational part of the arrival time surface closest to itself (see Fig. 1 and Eq. 2). The effect is to push the locations of the stationary points away from the source of external shear, hence increasing the radial spread of images. It follows that the direction of external shear can be determined by examining image locations with respect to the galaxy center. PG1115 is a good example; the image closest to the galaxy center, image B, is located between the galaxy-lens and the galaxy group, which is the source of external shear in this case.

(4) Position angles (PA) of images are determined by the ellipticity PA of the galaxy roughly at the radius of the images. When images are spread over a range of radial distances their PA provide information on galaxy ellipticity PA over a range of galactocentric distances. Thus detailed modeling can reveal the twisting of the isodensity contours.

(5) Not all types of information about images are equally useful for modeling purposes. The arrival time surface integrates over $`\kappa (\text{`})`$ twice, making time delays most sensitive to the overall mass distribution in the galaxy, and least dependent on the local small-scale perturbations in the mass distribution. Image positions are determined from the lensing equation which integrates over $`\kappa (\text{`})`$ once. Finally, image magnifications are very dependent on the local behavior of mass, making them the least useful for modeling. This means, unfortunately, that a double like Q0957, though it has well-measured substructure in the images and near-perfect time-delay measurements, provides too few constraints on the lensing mass to usefully estimate $`H_0`$ unless drastic assumptions about the mass distribution are made. In that case, the derived errors will tend to be underestimated as was noted by Bernstein and Fischer (1999) who constructed many types of parametric models for Q0957: ‘The bounds on $`H_0`$ are strongly dependent on our assumptions about a “reasonable” galaxy profile’.

(6) A linear rescaling of the arrival time and lens equations, i.e., multiplying both by a constant factor $`ϵ`$ will not alter the observable properties of images, image separations and relative magnification tensors. Physically the transformation amounts to rescaling the mass density of the lens by $`ϵ`$ and adding a constant mass density sheet. This transformation was first discussed by Gorenstein et al. (1988) with regard to modeling of Q0957, and later became known as the mass sheet degeneracy. Note that a mass sheet extending to infinity is not needed, a mass disk larger than the observed field is enough because an external monopole has no observable effect.

## 4 The method

The first step is to pixelate the lens plane mass distribution of the main lensing galaxy. In practice we use $`0.1^{\prime \prime }`$ pixels, and limit the galaxy to a circular window of radius about twice that of the image-ring. Pixelated versions of Eqs. (2) and (3) are:

$$\tau (\text{`},\text{`}_\mathrm{s})=(1+z_\mathrm{l})\frac{D_{\mathrm{ol}}D_{\mathrm{os}}}{D_{\mathrm{ls}}}\left[\frac{1}{2}|\text{`}|^2\text{`}\text{`}_\mathrm{s}\underset{n}{}\kappa _n\psi _n(\text{`})\right]$$
(5)   
and

$$\text{`}\text{`}_\mathrm{s}\underset{n}{}\kappa _n\stackrel{}{\alpha }_n(\text{`})=0,$$
(6)   
where the summation is over mass pixels and $`\psi _n`$ and $`\alpha _n`$ are integrals over individual pixels and can be evaluated analytically (see Appendix of SW97). A term $`|\text{`}_\mathrm{s}|^2`$ has been omitted from Eq. (3) because a constant additive factor in the arrival time cannot be measured.

Image properties translate into linear constraints in the $`(N+2)`$-dimensional model space, where $`N`$ dimensions represent a pixel each and 2 represent source coordinates. We call these primary constraints. The images can provide us with only a few constraints: in a 4-image system we have $`2\times 4`$ coordinates and 3 time delay ratios: 11 in all. On the other hand, the unknowns are numerous, $`20^2`$ mass pixels plus 2 source coordinates. This results in a plethora of galaxy models each of which reproduces the image properties exactly. Luckily, the bulk of these models can be discarded because they do not look anything like galaxies. In fact, we consider only those models which satisfy the following further (linear) constraints, which we call secondary. These pertain to the main lensing galaxy:

1. mass pixel values, $`\kappa _n`$ must be non-negative;
2. the location of the galaxy center is assumed to be coincident with that of the optical/IR image;
3. the density gradient of the lens must point no more than $`45^{}`$ away from the center of the galaxy;
4. the lens must have inversion symmetry, i.e., look the same if rotated by $`180^{}`$ \[enforced only if the main lensing galaxy appears unperturbed and has no companions close to the QSO images\];
5. logarithmic projected density gradient in the vicinity of the images, $`d\mathrm{ln}\kappa /d\mathrm{ln}\theta =\mathrm{ind}(r)`$ should be no shallower that $`0.5`$. For a power law projected density profile, radial magnification at an image is equal to $`1/\mathrm{ind}(r)`$, therefore a statement that $`\mathrm{ind}(r)<0.5`$ means that images are magnified radially by less than a factor of 2, which is probably a reasonable limit given the appearance of optical Einstein rings seen in some systems, for example, PG1115 and B0218;
6. external shear, i.e., influence of mass other than the main lensing galaxy is restricted to be constant across the image region, and is represented by adding a term $`\frac{1}{2}\gamma _1(\theta _1^2\theta _2^2)+\gamma _2\theta _1\theta _2`$ to the lensing potential, Eq. (2).

All these constraints are non-restrictive and are obeyed by the vast majority of galaxies, thus our analysis explores the widest possible range of galaxy mass distributions.

Obviously, the primary and secondary constraints are not enough to isolate a unique galaxy mass solution. A unique solution can be singled out by further specifying galaxy properties. For example, in SW97 particular galaxy models were found making a trial value of $`H_0`$ as one of the primary constraints, and then and picking the model that followed the observed light distribution as closely as possible given the rigid primary and secondary constraints, see Figures 2–5 of SW97. Here our aim is different.

Any of the infinitely many models remaining after the primary and secondary constraints have been applied could be the real lens, as all of them reproduce the image properties exactly and all look reasonably like a galaxies, therefore any one of the corresponding derived $`H_0`$’s could be the real $`H_0`$. We want to produce an ensemble that samples this model space, and our procedure is as follows.

The allowed models form a simplex in the $`(N+2)`$-dimensional space of mass pixels and source positions, because the constraints are all linear. We start with a random point in the allowed simplex (i.e., an allowed model). Next we choose a random vertex of that simplex, which is easily done by linear programming. Then we consider the line joining the current point with the vertex, and move to a random point on it, taking care to remain inside the simplex. The process is repeated until a sample of 100 models, and hence 100 $`H_0`$ values, is assembled. This procedure is a trivial case of the Metropolis algorithm (see e.g., Binney et al. 1992) for sampling density functions in high-dimensional spaces.

The resulting ensemble of $`H_0`$ values has a straightforward interpretation in terms of Bayesian probabilities. The part of model space allowed by the secondary constraints is the prior (i.e., possibilities allowed before considering tha data). Our prior is uniform, which is to say that we have not incorporated any prior preferences between different models allowed by the secondary constraints. Since the unknowns $`\kappa _n`$ occur linearly in Eqs. 5 and 6, a uniform prior means that any linear interval in $`\kappa _n`$ is a priori as probable as any other interval of equal length. The primary constraints come from data, and the 100 models that satisfy both primary and secondary constraints sample the posterior probability distribution. At the present time there is no clear motivation to use any other but a uniform prior, however, a non-uniform prior, if desired, would modify the method only slightly: one could either keep the same 100 models but weight them according to the prior, or take the prior into account while choosing the models through the Metropolis prescription.

## 5 Blind tests of the method

Before applying the method to real systems we try it on a synthetic situation designed to resemble the real world as close as possible. One of us, “Person A”, picked an $`h`$ value and created a set of four galaxies and the corresponding images of a single background source in each case. Exact values of image positions with respect to the galaxy center and time delays (but not $`h`$, nor information as to whether the galaxy was inversion symmetric or if there was any external shear) were conveyed to the other one of us, “Person B”, who used this information to construct an ensemble of galaxy models and derive $`h`$ distributions for each case separately.

We ran the whole experiment several times to remove bugs, and did not want to fall into the trap of simply publishing the results of the best run. So once we were confident that the experiment worked, we decided that the next four galaxies, whatever the results, would go into the published paper. Figure 2 pictorially illustrates the three stages of our blind test.

Person B applied the reconstruction method to each system twice, once with the assumption of inversion symmetry (i.e., symmetric galaxies, see item 4 in Section 4), and once without. Based on the appearance of the reconstructed mass distribution Person B decided whether inversion symmetry constraint was right in each case. Figures 310 present the results for each of the four galaxies. For galaxies #1, 3 and 4 Person B picked symmetric options, and the asymmetric option for galaxy #2. Panels (a) and (b) of Figures 3, 5, 7, and 9 show the actual projected density distribution and the average of the 100 reconstructed galaxies, for galaxies #1, 2, 3, and 4 respectively. In a map which is an average of many reconstructions, persistent features of individual maps are enhanced while peculiarities are washed out, so the average is a reasonable guess as to what the real galaxy looks like, in a probabilistic sense.

Panels (a) of Figures 4, 6, 8, and 10 plot the slope of density profile, $`\mathrm{ind}(r)`$ vs. derived $`h`$. The ‘real’ value of $`h`$ is 0.025. In all the cases the slope of the density profile, $`\mathrm{ind}(r)`$ in the vicinity of the images correlates with the derived $`h`$ value, though the degree of correlation and its slope is not universal. Qualitatively, the reason for the correlation is easily understood. A relatively flat galaxy density profile, i.e., $`|\mathrm{ind}(r)|`$ is small, translates into a flat gravitational contribution to the arrival time surface, and ‘fills’ the well of the geometrical time delay contribution evenly resulting in small fluctuations in the amplitude of the total arrival time surface. Thus the predicted time delays between images will be small, and to keep the observed time delays fixed the derived $`h`$ has to be small as well.

Panels (b) of Figures 4, 6, 8, and 10 show the derived $`h`$ probability distribution. These distributions look different for all galaxies, because galaxy morphologies are different. Since all four are independent probability distributions based on that galaxy, the overall distribution is just the product of the four, see Figure 11. The solid histogram is the product of the four distributions presented in panels (b) of Figures 4, 6, 8, and 10. The dashed histogram is similar, but results from Person B excluding what appeared to be the best constrained galaxy (# 3), and the dotted histogram represents the case where inversion symmetry was not applied to any of the systems. All three resultant distributions recover $`h`$ fairly well, with the $`90\%`$ of the models contained within $`20\%`$ of the true $`h`$. However the distributions are not the same; the most probable values are different by $`10\%`$. This illustrates how a relatively minor feature in modeling constraints, namely exclusion or inclusion of inversion symmetry, can make a considerable difference in the estimated $`h`$ value when the goal is to achieve precision of $`10\%`$. Based on this observation we conclude that assumed galaxy shape in parametric reconstructions plays a major role in determining the outcome of the $`H_0`$ determination.

How robust are the results to the changes in other modeling assumptions? Changing pixel size by a factor of $`1.5`$, and relaxing mass gradient angle constraint (item 3 in Section 4) does not change our results considerably.

## 6 Application to real systems

### 6.1 PG1115

Figure 12 shows the results of the reconstruction for PG1115. Since the main lensing galaxy has no close companions and its light profile is smooth we have included inversion symmetry as one of the modeling constraints. The average of 100 arrival time surfaces is shown in Figure 12(a); Figure 12(b) shows the corresponding caustics and critical lines. The latter are not as smooth as the former because locations of caustics and critical lines are derived using the gradients the arrival time surface, which are always noisier than the original function. Panels (c) and (d) plot the quantity $`\text{`}\text{`}_\mathrm{s}_n\kappa _n\psi _n(\text{`})`$, and the total arrival time surface, respectively. The plot of the modified gravitational potential, (c), illustrates the effect of external shear which is due to a galaxy group to the lower right of the main galaxy.

Because $`\mathrm{ind}(r)`$ has been measured for the main lensing galaxy in PG1115, the relation between profile slope and derived $`H_0`$, Figure 13(a), can be used to derive an upper limit on $`H_0`$. Impey et al. (1998) fit the galaxy light with a de Vaucouleurs profile of an effective radius $`r_e=0.59^{\prime \prime }`$. At the location of the images, about $`1.3^{\prime \prime }`$ from galaxy center the double logarithmic density slope is $`\mathrm{ind}(r)=2.3`$. Assuming that the mass profile can only be shallower than the light profile, and consulting Figure 13(a) we place an upper limit on $`H_0`$ of 75$`\mathrm{km}\mathrm{sec}^1\mathrm{Mpc}^1`$.

If the true mass density profile slope is isothermal the corresponding $`H_0`$ is 30$`\mathrm{km}\mathrm{sec}^1\mathrm{Mpc}^1`$. A low value of $`H_0`$ was also obtained by parametric models that assumed isothermal models for the galaxy (Schechter et al 1997).

In the blind test, Section 5, we assumed that all time delays are known precisely, which is not currently the case for any of the systems except Q0957. What effect does an error in time delay determination have on the derived $`H_0`$? Figure 13(b) shows two distributions derived using two different $`\mathrm{\Delta }t`$ determinations based on the same lightcurves. There is a $`20\%`$ difference in the most probable value of $`H_0`$ in the two histograms, but overall they are not very different. Both distributions are very broad; $`90\%`$ of the models span the range between 30 and 75$`\mathrm{km}\mathrm{sec}^1\mathrm{Mpc}^1`$.

Figure 14 shows a dense version of the arrival time surface. The regions of the plot where the contours are sparse are the flattest, i.e., most ‘stationary’ regions of the lens plane. This is where one would expect to find images of sources placed close to the main source. For example if the point-like QSO is surrounded by a host galaxy the image of that galaxy will be well delineated by these ‘empty’ regions. In fact, the observed optical ring in the case of PG1115 is well reproduced by the ring in Figure 14.

### 6.2 B1608

The light distribution of the lensing system is rather messy, possibly representing a merging/interacting galaxy pair, therefore inversion symmetry was not used in the following reconstructions. Figures 15 and 16 are the same as Figures 12 and 13, but for B1608. The range 50 to 100$`\mathrm{km}\mathrm{sec}^1\mathrm{Mpc}^1`$ in Figure 16(b) encompasses about 90$`\%`$ of the reconstructions.

### 6.3 Combined p(h) plot

Just like in the case of the blind test we now multiply probability distributions from PG1115 and B1608 to get the combined distribution, Figure 17. $`90\%`$ of all points lie within the range 43–79$`\mathrm{km}\mathrm{sec}^1\mathrm{Mpc}^1`$, while the median of the distribution is 61$`\mathrm{km}\mathrm{sec}^1\mathrm{Mpc}^1`$. Note that the errorbars obtained using our method are substantially larger than what is usually quoted in other studies. We ascribe this increase to the more systematic sampling of the whole image-defined model space unrestricted by the confines of parametric models.

## 7 Discussion and Conclusions

Multiply-imaged QSO systems provide us with an elegant way of measuring $`H_0`$, and a lot of observational and modeling effort has been invested in this enterprise. As the quality of the observational data improves, most of the uncertainty in $`H_0`$ is contributed by the mass distribution in the lens. How to treat this problem is a matter of some debate. Should one use a single, physically motivated mass model, or should one approach the problem with no preconceptions about the galaxy form?

In general, ad hoc restrictions on the allowed mass models translate into a too optimistic and probably biased estimated distribution of $`H_0`$’s. To avoid this trap one has to allow as much freedom for the lens models as possible. On the other hand, to yield a useful estimate of $`H_0`$ one has to restrict the amount of freedom allowed for the models using physically motivated criteria. Ideally one wants to balance these two opposing tendencies and impose just the correct quantity and quality of model constraints. Based on our experience from the present work we conclude that parametric, or any other approach that severely restricts the freedom of the galaxy-lens, has over-constrained their models and thus ended up with unrealistically small errorbars, and biased $`H_0`$’s. As a result different models of the same systems can yield discrepant results. For example, Romanowski & Kochanek (1998) use dynamical methods to model the galaxy in Q0957, and further constrain the galaxy to be similar to nearby ellipticals; they quote $`61_{15}^{+13}`$ at $`2\sigma `$ level. Bernstein & Fischer (1999) analyzed the same system but used a range of astrophysically reasonable parametric models. Their estimate, $`77_{24}^{+29}`$, also at $`2\sigma `$ level, does not agree with that of Romanowski & Kochanek. Our approach is different in that it does not presuppose a galaxy shape, but instead allows us to impose as many or as few constraints as is deemed appropriate.

The most unrestricted models would be constrained solely by what we call the primary constraints, i.e., image observables. By definition these would yield unbiased estimates of $`H_0`$ based on lensing data alone. We chose to go somewhat beyond this and apply what we call secondary constraints, which describe realistic galaxies in most general terms. The derived $`H_0`$ distributions are narrower; the price we pay is a small amount of bias. It can be argued that we are still too generous with our mass models, i.e. other galaxy characteristics can be safely assumed, and hence tighter constraints can be applied to the models without sacrificing the unbiased nature of results. This avenue can be taken in the future work if additional constraints become available.

A potential source of additional modeling constraints are optical rings, lensed images of the QSO galaxy host, which are seen in some cases, for example, in PG1115 and B0218. The orientations and elongations of individual images of QSO host galaxy can be used as linear inequality constraints to narrow down the range of possible galaxy mass distributions.

If two or more sources with known redshifts were lensed by the same foreground galaxy, these could be used to break the mass sheet degeneracy and thus further constrain the galaxy. However in practice cases of two sources at different redshifts lensed by the same galaxy are expected to be very rare because of the small galaxy cross-sections.

Probably the most promising potential constraint is based on the relation between the slope of the projected density profile around the images and the derived $`H_0`$. If the slope can be estimated by means other than lensing, or at least a limit placed on its value as we did in Section 6.1 using the observed slope of the light distribution, then $`H_0`$ can be constrained much better compared to what is currently possible.

With the two systems used in the present work, PG1115+080 and B1608+656, and implementing primary constraints of image properties and secondary constraints describing a few general properties of lensing galaxies, we conclude that $`H_0`$ is between 43 and 79$`\mathrm{km}\mathrm{sec}^1\mathrm{Mpc}^1`$ at $`90\%`$ confidence level, with the best estimate being 61$`\mathrm{km}\mathrm{sec}^1\mathrm{Mpc}^1`$.