text:
# Degeneracies and scaling relations in general power-law models for gravitational lenses

## 1 Introduction

Gravitational lenses offer a unique possibility to determine cosmological distances and hence the Hubble constant (Refsdal, 1964). The method avoids using any form of distance ladder and is almost independent of a deeper understanding of astrophysical processes. This has the advantage that possible statistical and systematic uncertainties can be controlled much better than with other methods. Besides the measurement errors of redshifts, positions and time delays, the most important source of uncertainties in using gravitational lenses to determine $`H_0`$ is the mass model of the lens. Constraints provided by observational data are never sufficient to fix the mass distribution uniquely (Saha & Williams, 1997; Williams & Saha, 2000). Simple parametric models, based on the knowledge about typical galaxies in the local universe or desired mathematical properties, are normally used to overcome this difficulty. This approach, while leading to consistent results, has the disadvantage of hiding the underlying uncertainties and making it difficult to quantify them.

In this paper, we consider a family of mass models with separable radial and angular dependence of the potential. External perturbations are included as an external shear. In this way, the different parts can be studied independently. For the radial dependence, we choose a simple power-law and are especially interested in the influence the radial slope has on the results. Several authors (Wambsganß & Paczyński, 1994; Witt et al., 1995; Wucknitz & Refsdal, 2001) have found approximate scaling relations of $`H_02\beta `$ for spherical models with external shear. The same dependence was found before by Chang & Refsdal (1976) for doubles in spherical models without shear (see also Refsdal & Surdej, 1994). The work of Williams & Saha (2000) showed, that the radial slope also has important effects in non-parametric models. Due to the very general nature of these models, no exact scaling law could be found for them.

We follow an intermediate approach by using an arbitrary function only for the angular part. The models include elliptical mass distributions and potentials as well as other shapes. It has been observed in several lens systems, that the external shear required to fit the data with simple elliptical models is much higher than expected. This might indicate a kind of asymmetry of the galaxy that cannot be accounted for by elliptical models. The general angular part of the potential we use can describe this ‘internal shear’ more accurately than simple parametric approaches.

We extend and generalize the work of Witt et al. (2000) by including the three independent time delays in quadruple systems as constraints for the models. We do not use magnification ratios for several reasons. One is the observational problem of reliably measuring correct flux ratios. Fluxes are influenced by microlensing and extinction. These effects can be very strong in the optical and in some cases they are still significant for radio wavelengths. The other reasons are related to our formalism and will be discussed later.

With this approach, we can derive explicit solutions for the Hubble constant as a function of the observables, the power-law exponent $`\beta `$, and the external shear. We use the results to find a simple but rigorous scaling law describing the dependence of $`H_0\mathrm{\Delta }t`$ on $`\beta `$ in lenses with quadruple images. This scaling law of $`H_0\mathrm{\Delta }t(2\beta )/\beta `$ is also valid for more special parametric power-law models within the allowed range of $`\beta `$ and is therefore inherent to these models and not an artefact of the general angular dependence we discuss.

The qualitative effect of this scaling law is easily understood when comparing it with the $`2\beta `$ scaling discussed before and with the mass-sheet degeneracy, which leads to $`H_0\mathrm{\Delta }t1\kappa `$. In all cases, a shallower density profile (larger $`\beta `$ or $`\kappa `$) leads to smaller values of $`H_0\mathrm{\Delta }t`$. The flat part of the density distribution ($`\kappa `$ in the mass-sheet degeneracy) amplifies the deflection angles but leaves the time delays unaffected. To fit the given geometry, the total deflection angles have to be constant, therefore the time delays ($`H_0\mathrm{\Delta }t`$) decrease.

We also discuss the effect of the external shear on the time delays and the measured Hubble constant. If the shear is changed, the internal asymmetries of the mass distribution have to be adjusted to compensate for the changes. The total effect of these two contributions will be described by the new concept of a ‘critical shear’ $`𝜸_\mathrm{c}`$. The measured Hubble constant is a linear function of the external shear and becomes zero for $`𝜸=𝜸_\mathrm{c}`$. We will present a simple interpretation of the critical shear in terms of the image geometry.

The main goal of this work is to investigate the uncertainties in measurements of $`H_0`$ or cosmological parameters from time delays in gravitational lenses. The combination of distances that can be determined not only depends on the redshifts and $`H_0`$ but also on the other cosmological parameters ($`\mathrm{\Omega }`$ and $`\lambda `$ in standard cosmology) and the smoothness of the mean density distribution (‘lumpiness’ parameter $`\alpha `$). The results will also help in minimizing possible errors, either by selecting lenses with the least uncertainties or by using constraints of the model parameters which are lensing-independent.

Finally, we apply the formalism to several known lens systems, with and without measured time delays, to compare our analytical results with those from numerical model fitting using parametric models. For 2237+0305, we compare with own numerical models and find a very good agreement even though the time delays are not used as constraints in this case.

Our results can be used directly to determine $`H_0`$ from time delays without explicit modelling but should not be used as a substitute for it. They are rather meant to provide an explanation for the degeneracies and scaling relations that have been found for several families of lens models. Nevertheless, we show that the direct application to real systems is possible.

In the appendix, we discuss possible Einstein rings in our family of models. Interestingly, infinitely small sources can be mapped as elliptical rings for arbitrary values of the external shear. The ellipticity of these rings is directly determined by the shear.

## 2 The lens model

For the lens equation and deflecting potential $`\psi `$, we use the following notation:

$`𝒛_𝐬`$ $`=𝒛𝜶(𝒛)`$ (1)
$`𝜶(𝒛)`$ $`=\mathbf{}\psi (𝒛)`$ (2)
$`𝒛`$ $`=r(\mathrm{cos}\theta ,\mathrm{sin}\theta )`$ (3)
We use a power-law approach with a general separable angular part for the main lensing galaxy.

$$\psi ^{(\mathrm{g})}(𝒛)=r^\beta F(\theta )$$
(4)
This family of models includes both elliptical mass distributions and elliptical potentials with arbitrary radial mass index $`\beta `$. In the following, we assume all image positions as known. For the chosen model, this implies also the knowledge of the position of the galaxy centre itself, because coordinates are relative to this centre. Later we will see, however, that some of the equations are translation invariant, leaving the results unchanged when shifting the galaxy. The limiting cases of $`\beta 0`$ and $`\beta =2`$ correspond to a point mass and a density constant in the radial direction, respectively. Between these two, the isothermal models are another special case with $`\beta =1`$. The normalized surface mass density of the general model is

$`\sigma (𝒛)`$ $`={\displaystyle \frac{1}{2}}\mathbf{}^2\psi `$ (5)
$`={\displaystyle \frac{r^{\beta 2}}{2}}\left(\beta ^2F(\theta )+F^{\prime \prime }(\theta )\right)\text{.}`$ (6)
Here and in the rest of this paper, primes indicate derivatives with respect to $`\theta `$. A very simple relation holds for the radial derivative of the potential which we will need in the time delay equations later:

$$𝒛\mathbf{}\psi ^{(\mathrm{g})}=\beta \psi ^{(\mathrm{g})}$$
(7)
To account for nearby field galaxies or the contribution of a cluster, we include an external shear plus an additional constant mass density or convergence $`\kappa `$ into our models. This $`\kappa `$ can not be determined from observations as a consequence of the so called mass-sheet degeneracy, first discussed by Falco et al. (1985) and Gorenstein et al. (1988). In the following, we therefore always use a fixed $`\kappa `$. We parametrize the shear as follows:

$`\psi ^{(\gamma )}(𝒛)`$ $`={\displaystyle \frac{1}{2}}𝒛^\text{t}𝒛`$ (8)
The shear matrix includes both the convergence and the shear itself.
$`=\left(\begin{array}{cc}\kappa \gamma _x& \gamma _y\\ \gamma _y& \kappa +\gamma _x\end{array}\right)`$ (9)
$`𝜸`$ $`=\gamma (\mathrm{cos}2\theta _\gamma ,\mathrm{sin}2\theta _\gamma )`$ (10)
With this definition, $`\theta _\gamma `$ points in the direction of the external perturbing mass or equivalently in the opposite direction. Note that this part of the potential is a special case of the power-law with $`\beta =2`$.

## 3 Time delays

The light travel time $`t`$ for a certain image at $`𝒛`$ can be calculated for arbitrary lens models by using the equation

$$t=\frac{1}{c}D_{\mathrm{eff}}\left(\frac{1}{2}|𝜶(𝒛)|^2\psi (𝒛)\right)+\text{const}\text{with}$$
(11)
$$D_{\mathrm{eff}}=\frac{D_\mathrm{d}D_\mathrm{s}}{D_{\mathrm{ds}}}(1+z_\mathrm{d})\text{.}$$
(12)
Here $`D_\mathrm{d}`$ and $`D_\mathrm{s}`$ are angular size distances of the deflector and the source from the observer, while $`D_{\mathrm{ds}}`$ is the distance of the source measured from the deflector. We separate this equation as usual into one part depending only on cosmology (without $`H_0`$), the contribution from the Hubble constant, and a lens dependent part. This is done by using reduced distances $`d_x`$ (with $`D_x=d_xc/H_0`$) that only depend on the cosmological parameters (and the redshifts) but not on $`H_0`$ itself. We now use a scaled version of the Hubble constant,

$$h=\frac{H_0}{d_{\mathrm{eff}}}\text{,}$$
(13)
which includes the cosmological factors and will later be determined by the lensing effect. This $`h`$ is directly proportional to the Hubble constant itself when the cosmology is fixed. For two images $`i`$ and $`j`$, we get a time delay of $`\mathrm{\Delta }t_{ij}`$ with

$`h\mathrm{\Delta }t_{ij}`$ $`=h(t_it_j)`$ (14)
$`={\displaystyle \frac{1}{2}}\left(|𝜶_i|^2|𝜶_j|^2\right)(\psi _i\psi _j)\text{.}`$ (15)
By using the lens equation (1), we can transform this expression into a linear functional of the deflecting potential and its derivatives. Mixed terms like $`𝒛_i𝜶_j`$ for $`ij`$ can be eliminated, so that the resulting time delay can again be written as the difference of the light travel times $`t_i`$ themselves.

$$ht_i=\frac{1}{2}r_i^2+𝒛_i\mathbf{}\psi _i\psi _iC$$
(16)
Although it is not immediately noticeable, this equation still is invariant under translations of coordinates. To prove this, one has to apply the lens equation.

For the special potential we discuss here, the relation (7), which also holds for $`\psi ^{(\gamma )}`$ with $`\beta =2`$, can be used to eliminate the derivatives of the potential and obtain a simple expression only depending on both parts of the potential at the image position:

$$ht_i=\frac{1}{2}r_i^2(1\beta )\psi _i^{(\mathrm{g})}+\psi _i^{(\gamma )}C$$
(17)
We retain using the notion of light travel times for each image instead of time delays between the images to keep the equations simple. The $`t_i`$ are defined except for an overall additive constant that gets absorbed into the constant $`C`$.

The last equation was already presented in Witt et al. (2000) for the special cases of $`\beta =1`$ including external shear and for the general shearless power-law model.

## 4 Counting constraints and parameters

Before solving the equations, we want to discuss how many parameters of the model can at most be determined from observations of lens systems with $`n`$ images, especially $`n=4`$. See Table 1 and 2 for a list of parameters and constraining equations. As the time delays do not change when adding the same constant to all $`t_i`$, we have to include the constant $`C`$ as parameter and thus have $`n`$ constraints with one parameter for this. We might as well have used only the (uniquely defined) $`n1`$ time delays without $`C`$. Both possibilities result in $`n1`$ more constraints than parameters just for the time delays.

Even for systems with 4 images of one source, it will be impossible to determine all parameters. We therefore decide to fix $`\beta `$ in the following calculations so that the results can be used to study the dependence of $`h`$ on the radial mass slope. We will see that, with fixed $`\beta `$, all equations stay linear.

Another critical parameter is the external shear, which seems to be higher than expected in many detailed lens models. Fixing $`\gamma `$ for the moment, we will be able to investigate the influence of any uncertainties in the external shear. In the special case of $`\beta =1`$, it will even be possible to determine $`𝜸`$ from the constraints, because a number of other parameters do not contribute then.

We do not include fluxes or magnifications, because they would provide only $`n1`$ constraints (the independent flux ratios) and at the same time add $`n`$ more parameters (the second derivatives of the angular part of the potential at the image positions $`F_i^{\prime \prime }`$). The effect, that models get less constrained, when more observations are included in the analysis, is of course unknown in normal parametric models, where the number of parameters is fixed. Our models have an infinite number of parameters, of which only a finite subset is needed to compare with observations. The number of relevant parameters can change when we include more constraints.

Another reason (besides the difficulties in determining fluxes already discussed) is that, in contrast to deflection angles and time delays, magnification ratios are non-linear functionals of the lensing potential, complicating the equations considerably.

## 5 Lens equations

To exploit the information given by the image positions, we have to insert the power-law potential with shear into the lens equation (1). The derivative of the galaxy part of the potential can be obtained most simply by rotating its polar form to Cartesian coordinates using the transformation matrix

$$\left(\begin{array}{c}_x\\ _y\end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}\theta & \mathrm{sin}\theta /r\\ \mathrm{sin}\theta & \mathrm{cos}\theta /r\end{array}\right)\left(\begin{array}{c}_r\\ _\theta \end{array}\right)\text{.}$$
(18)
Written in a form to take into account the role of $`𝜸`$, $`F_i`$ and $`F_i^{}`$ as unknowns of the equations, this leads to the following equation.

$$𝒛_𝐬=(1\kappa )𝒛_ir_i^{\beta 2}\left(\begin{array}{cc}x_i& y_i\\ y_i& x_i\end{array}\right)\left(\begin{array}{c}\beta F_i\\ F_i^{}\end{array}\right)+\left(\begin{array}{cc}x_i& y_i\\ y_i& x_i\end{array}\right)\left(\begin{array}{c}\gamma _x\\ \gamma _y\end{array}\right)$$
(19)
It is easily seen, that this set of $`2n`$ equations can be used to determine $`F_i`$ and $`F_i^{}`$, assuming $`𝜸`$ as known:

$`\beta F_i`$ $`=r_i^\beta \left((1\kappa )r_i^2x_ix_\mathrm{s}y_iy_\mathrm{s}+\gamma _x(x_i^2y_i^2)+2\gamma _yx_iy_i\right)`$ (20)
$`F_i^{}`$ $`=r_i^\beta \left(y_ix_\mathrm{s}x_iy_\mathrm{s}2\gamma _xx_iy_i+\gamma _y(x_i^2y_i^2)\right)`$ (21)
## 6 The general set of linear equations

We now use (20) to express the galaxy potential in the time delay equations (17). We decide to use the light travel times themselves rather than the $`H_0`$-independent time delay ratios as constraints. In this way we can keep the equations linear and the analysis much simpler. As a result, the following set of equations is obtained for $`i=1\mathrm{}n`$:

$$\begin{array}{cc}\hfill (1\kappa )r_i^2& +(x_i^2y_i^2)\gamma _x+\mathrm{\hspace{0.33em}2}x_iy_i\gamma _y\hfill \\ & +\frac{2\beta }{2\beta }(t_ih+C)=\mathrm{\hspace{0.33em}2}\frac{1\beta }{2\beta }𝒛_i𝒛_𝐬\hfill \end{array}$$
(22)
The most interesting fact, besides the linearity, is the simple way the mass index $`\beta `$ appears in the equations. By combining information from the time delay and lens equations the way we did, it was possible to remove the terms with $`\beta `$-dependent exponents. Now $`\beta `$ only contributes in the scaling factors of the unknown parameters. Having solved the system for one value of $`\beta `$, we immediately find the general solution by scaling $`h`$, $`𝒛_𝐬`$ and $`C`$ with the appropriate factors.

## 7 The isothermal model

In the case $`\beta =1`$, the equations (22) obviously degenerate with respect to the source position $`𝒛_𝐬`$. It is then impossible to constrain the latter, but the same information can now be used to determine the external shear. In this case, the inclusion of the lens equations to determine the parameters $`F_i`$ of the galaxy potential was not really necessary, because $`\psi ^{(\mathrm{g})}`$ does not contribute to the time delay equations (17). We can now directly invert the latter or (22) to obtain solutions for $`h_{\mathrm{iso}}`$ and $`𝜸_{\mathrm{iso}}`$.

The transition from exactly isothermal to nearly isothermal models deserves some attention. Even for models with $`\beta `$ differing only slightly from unity, $`𝜸`$ is a free parameter, while it is fixed by the observational data for $`\beta =1`$. An incorrect external shear for almost isothermal models is usually compensated for by diverging source positions and $`F_i`$ and $`F_i^{}`$, leading to unrealistic models. It therefore seems appropriate to use the correct isothermal shear even for models that do not exactly obey $`\beta =1`$. However, one has to take into account possible measurement uncertainties that introduce errors into $`𝜸_{\mathrm{iso}}`$. Especially the time delays, which all contribute to the solution, can introduce significant uncertainties.

## 8 Solution for the general model

In the general case with $`\beta 1`$, (22) can be solved directly to determine $`h`$ for a given shear $`𝜸`$. Even without writing the explicit solution, we see, that

$$h\frac{2\beta }{\beta }(1+g_x\gamma _x+g_y\gamma _y)$$
(23)
for some constants $`g_x`$ and $`g_y`$. The Hubble constant scales with $`(2\beta )/\beta `$ and is a linear (but not proportional) function of the shear. Notably, the isothermal model plays no special role in this equation, despite the fact that, strictly speaking, $`𝜸`$ can not be chosen freely in this case. As the isothermal shear is usually only weakly constrained due to the limited measurement accuracy, we may, however, use different values of $`𝜸`$ even for $`\beta =1`$. This is always true when the constraints we used here are not all available. Considering this, one might well apply equation (23) regardless of the value of $`\beta `$.

## 9 The ‘Critical shear’

As the Hubble constant is linear in $`𝜸`$, there has to be a one-dimensional family of values of the shear with vanishing $`h`$. The shear with the smallest absolute value from this family will now be called the ‘critical shear’ $`𝜸_\mathrm{c}`$. If we denote the shearless value with $`h_0`$, we can write the effect of external shear as

$$\frac{h}{h_0}=1\frac{𝜸𝜸_\mathrm{c}}{\gamma _\mathrm{c}^2}\text{.}$$
(24)
If $`𝜸`$ and $`𝜸_\mathrm{c}`$ point in the same direction, this scaling factor can be written as $`1\gamma /\gamma _\mathrm{c}`$ which is analogous to the scaling of $`1\sigma /\sigma _\mathrm{c}`$ for an additional mass sheet $`\sigma `$ with critical density $`\sigma _\mathrm{c}`$. The critical shear does not dependent on $`\beta `$ or the time delays $`t_i`$ and can be calculated from the image positions alone.

If an upper bound of $`\gamma _{\mathrm{max}}`$ can be assumed, this translates to a range of

$$h_{\left(\begin{array}{c}\mathrm{max}\\ \mathrm{min}\end{array}\right)}=\left(1\pm \frac{\gamma _{\mathrm{max}}}{\gamma _\mathrm{c}}\right)h_0\text{.}$$
(25)
The critical shear is thus a measure for the amount of external shear that can change $`h`$ significantly. The larger it is, the smaller the influence of uncertainties in the shear on the determined Hubble constant. Shear in the direction of $`𝜸_\mathrm{c}`$ contributes maximally, in a perpendicular direction the influence vanishes. With the ‘direction of shear’, we mean in this context the orientation $`2\theta _\gamma `$ of the vector $`𝜸`$. This is not the direction of $`\theta _\gamma `$ towards the perturbing mass, cf. eq. (10). External field galaxies located perpendicular to the direction of maximal effect change $`h`$ by the same amount but with opposite sign. There are, however, four directions where external masses do not contribute. See section 13 for actual numbers of the critical shear in real observed lens systems.

### Geometrical interpretation

For constant light travel times or $`h=0`$, equation (22) describes an ellipse whose axes $`a`$ and $`b`$ are related to the external shear by

$$\frac{\gamma }{1\kappa }=\frac{a^2b^2}{a^2+b^2}\text{.}$$
(26)
For $`\gamma >1\kappa `$, this becomes a hyperbola. The position angle of the minor axis is the same as that of the perturbing mass responsible for the shear ($`\theta _\gamma `$). This means, that each ellipse/hyperbola passing through all images corresponds to a value of $`𝜸`$ with $`h/h_0=0`$. According to equation (24), for these values

$$𝜸𝜸_\mathrm{c}=\gamma _\mathrm{c}^2$$
(27)
holds. The values of $`𝜸`$ for all these conic sections span the complete subspace for which $`h/h_0=0`$ (see Fig. 1). The critical shear is the smallest of these values and can therefore be calculated with equation (26) for the ‘roundest’ ellipse passing through all images. In Figure 2, we show several such ellipses for some lens systems.

We conclude, that for certain values of $`𝜸`$ (see Figure 1), e.g. for $`𝜸=𝜸_\mathrm{c}`$, all time delays vanish. Each of the fitting ellipses can act as an Einstein ring for the correct value of $`𝜸`$. The light travel time is the same for all parts of such rings as Fermat’s theorem requires. In our consideration, the potential is only constrained at the image positions and arbitrary in other directions. The Einstein ring may therefore break up and form a number of discrete images with still vanishing time delays (see appendix).

## 10 Shifting the lensing galaxy

Surprisingly, most of the quantities we determined do not depend on the position of the lens centre. A shift of the centre is equivalent to adding a constant displacement to the vectors $`𝒛_i`$ and $`𝒛_𝐬`$. If we look at the general set of equations (22), we see that such a shift adds terms linear in $`𝒛_i`$ and a constant to the equations. The constant term can be absorbed in $`C`$ and (for $`\beta 1`$) the linear terms in $`𝒛_𝐬`$. As $`C`$ and $`𝒛_𝐬`$ are of no interest, the equations do not change when applying this shift.

This means, that $`h`$ and the critical shear are translation invariant and can be determined even if the lens position is not known. This is only true for the family of models we analyse here. Simple parametric models usually only fit the data for a specific position of the lens centre.

## 11 Spherical models for nearly Einstein ring systems

For spherical models, the equations are overdetermined. There are nevertheless systems, which can be fit accurately with this kind of model. It can be shown, that $`h_{\mathrm{sph}}/h_0`$ becomes zero for $`\beta 0`$ to cancel the vanishing denominator in (23) and ensure a finite $`h_{\mathrm{sph}}`$. It can also be shown, that $`𝜸`$ is parallel to $`𝜸_\mathrm{c}`$ for arbitrary $`\beta `$. Taken together, this means that for point mass models, we always have $`𝜸=𝜸_\mathrm{c}`$. The geometrical interpretation can be used to determine the direction of external shear for spherical models without calculations. It is given by the orientation of the minor axis of the roundest ellipse passing through all images. For the systems discussed below (and for point mass lenses), even the absolute value of $`\gamma `$ can be determined from this ellipse.

For systems, where the $`n`$ images are all located very close to the Einstein ring at $`r_0`$, we can recover another scaling relation. In this case, the power-law can be interpreted as a local approximation to any radial mass profile, like softened power-law spheres or other models. We assume, that one fitting reference model is known. It is then possible to find a family of other models which are also consistent with the observations.

This was done numerically by Wambsganß & Paczyński (1994) and led to a scaling of $`h_{\mathrm{sph}}2\beta `$. Wucknitz & Refsdal (2001) presented a simple interpretation of this fact in terms of the well known mass-sheet degeneracy (Gorenstein et al., 1988). If we multiply the lens equation $`𝒛_𝐬=𝒛𝜶`$ with $`1\kappa `$, we get another lens equation with the source position and lensing potential (or deflection angle) scaled with the same factor, plus an additional constant convergence $`\kappa `$. This means, that lens models given by the deflection angles $`𝜶`$ and $`\stackrel{~}{𝜶}=(1\kappa )𝜶+\kappa 𝒛`$ are equivalent, but source positions, time delays, etc. are scaled by $`1\kappa `$ in the latter model.

As reference model, we choose an exponent $`\beta _0=1`$. This model is now transformed as just described and then approximated locally near $`r_0`$ by a modified power law with exponent $`\beta `$. Figure 3 illustrates this idea. For the new exponent, we find

$$2\beta =1\kappa \text{.}$$
(28)
As time delays, source position and external shear all scale with $`1\kappa `$, we find the general scaling laws

$$𝒛_𝐬2\beta ,h_{\mathrm{sph}}2\beta ,𝜸2\beta $$
(29)
for an arbitrary mass index $`\beta `$. The time delay ratios do not change.

With $`𝜸=𝜸_\mathrm{c}`$ for point mass systems, the shear for arbitrary $`\beta `$ becomes $`𝜸=𝜸_\mathrm{c}(2\beta )/2`$, leading to a ratio of the Hubble constants for the spherical and shearless case of

$$\frac{h_{\mathrm{sph}}}{h_0}=\frac{\beta }{2}\text{.}$$
(30)
The two models differ significantly for a realistic range of $`\beta `$.

## 12 Influence of the radial mass index $`𝜷`$

One of the interesting properties of equation (23) is the very simple dependence of the results on $`\beta `$. The determined value of the Hubble constant $`h`$ simply scales with the factor $`(2\beta )/\beta `$. The most alarming fact is, that this factor does not depend on the geometry of the lens, the time delay ratios or the amount of external shear. When using the models described here to determine $`H_0`$ from lens time delays, the error due to the assumption of an incorrect $`\beta `$ will be exactly the same for all lens systems as long as the real $`\beta `$ is more or less equal for all lensing galaxies.

Witt et al. (2000) numerically found a scaling of $`(2\beta )/\beta `$ in the case of a power-law model without external shear for orthogonal image pairs ($`\theta _i\theta _j90^{}`$), while these computations lead to $`2\beta `$ for opposed images. In contrast to their work, we have used all time delays as constraints so that they cannot scale differently. The common scale factor of all time delays shows as a scaling of $`h`$ in our calculations.

The reader might feel as uneasy about the seemingly diverging time delays or $`h`$ in the limit $`\beta 0`$ as the author did. This limit is equivalent to point mass models and one should not observe diverging time delays for this kind of lens. The point causing trouble here is the fixed external shear in our considerations. To fit the data with a point mass model, the shear has to be equal to the critical shear $`𝜸_\mathrm{c}`$. Taking this into account, the result gets multiplied by a vanishing $`h/h_0\beta `$ which cancels the $`1/\beta `$ factor. If we now change the shear by a small amount, the potential and thus $`h`$ will immediately diverge in the limit of small $`\beta `$. This has no direct physical implications, because the mass models will become extremely unphysical to compensate for the shear effect. Small relative differences in the $`F_i`$, that will be introduced by an incorrect shear, will lead to enormous asymmetries in the mass distribution. This is related to the fact, that realistic compact mass distributions can provide only almost spherically symmetric potentials. Any multipole moments would radially decrease more rapidly than the monopole term and must be very strong to have any effect. We will discuss this in more detail for a special case in section 13.1.

For spherical systems with images near the Einstein ring, we confirmed the approximate scaling of $`h_{\mathrm{sph}}2\beta `$. This seems to be incompatible with the general scaling law of $`h_0(2\beta )/\beta `$ at first sight. With the factor of $`h_{\mathrm{sph}}/h_0=\beta /2`$ from equation (30), by which the shear changes the result in spherical models, these two results are, however, in perfect agreement.

## 13 Application to special cases

To illustrate the results, we have presented here, and to test their relevance for real lenses, we want to apply the formalism to systems with a special Einstein cross like symmetry and to some of the known real systems that are either useful to actually determine $`H_0`$ (1115+080, 0911+0551, 1608+656) or are interesting because they are very well studied systems like 2237+0305. The calculations will show, that both scaling relations ($`h_{\mathrm{sph}}2\beta `$ and $`h_0(2\beta )/\beta `$) are relevant for the determination of $`H_0`$ from time delays. The detailed numerical models for 2237+0305 will furthermore show, that the scaling also applies if the time delays themselves are not used to constrain the models. All time delays scale almost exactly as predicted by our analytical work. We will also see that parametric models may fit only for a limited range of $`\beta `$. The scaling relations are then valid only within this range.

Besides the effects of the exponent $`\beta `$, the possible strong effects of any external shear will also be confirmed by the numerical models.

### 13.1 Symmetric Einstein cross like systems

A rather special example of systems shall be discussed explicitly in this section. We consider a lens with time delays $`\mathrm{\Delta }t_{12}=\mathrm{\Delta }t_{34}=0`$ and the following image positions:

$`x_{1,2}`$ $`=\pm (1\epsilon )r_0`$ $`y_{1,2}`$ $`=0`$ (31)
$`x_{3,4}`$ $`=0`$ $`y_{3,4}`$ $`=\pm (1+\epsilon )r_0`$ (32)
From these data, we immediately conclude $`𝒛_𝐬=0`$. The $`y`$-component of the shear does not contribute at all and cannot be constrained. The equation determining the Hubble constant now reads

$$h\mathrm{\Delta }t_{13}=(1\kappa )r_0^2\frac{2\beta }{\beta }2\epsilon \left(1\frac{\gamma _x}{\gamma _\mathrm{c}}\right)\text{,}$$
(33)
with the critical shear

$$\gamma _\mathrm{c}=(1\kappa )\frac{2\epsilon }{1+\epsilon ^2}\text{.}$$
(34)
The special symmetry makes it possible to choose any value of external shear even for isothermal models. Furthermore, it is possible to exactly reproduce the data with spherically symmetric models plus external shear. In this case, the shear is uniquely defined:

$`\gamma `$ $`=(1\kappa )(2\beta )\epsilon +𝒪(\epsilon ^2)`$ (35)
The time delay equation now becomes somewhat more complicated than in the non-spherical case. To first order in $`\epsilon `$, it reads

$`h_{\mathrm{sph}}\mathrm{\Delta }t_{13}`$ $`=(1\kappa )r_0^2(2\beta )\epsilon +𝒪(\epsilon ^2)\text{.}`$ (36)
Comparing this with (33), we recover the factor of $`\beta /2`$ between the spherical and shearless case. Both models fit the data exactly and, in this special case, the models are even compatible with highly popular elliptical mass distributions. That means that without any independent information about the external shear (or equivalently the ellipticity of the galaxy itself), we have a factor of two uncertainty even when only considering these two simplest models for $`\beta =1`$. The real situation may be even much worse, when we consider models with internal and external shear. In this case, any small unknown contribution of external shear of the order of $`\gamma _c`$ (which for very symmetric systems becomes arbitrarily small) will change the result significantly.

Witt et al. (1995) discussed exactly the same type of systems with spherical models plus shear. With $`\gamma `$ fixed, they also derived a scaling law of $`(2\beta )/\beta `$ (see their equation 8<sup>1</sup><sup>1</sup>1The exponent of the first term in equation (8) in Witt et al. (1995) is incorrect, it should be the same as that in the second term. (S. Mao, private communication)). When the shear is constrained by the lens equations, the scaling changes to the $`2\beta `$ form.

We finally want to discuss the consequences of diverging time delays in the fixed shear case for $`\beta 0`$ due to equations (23) and (33). For simplicity, we assume $`𝜸=0`$, $`\kappa =0`$ and $`\epsilon 1`$, but the argument is generally also true for other values. We write the potential as a multipole expansion<sup>2</sup><sup>2</sup>2We use an expansion for the principal axes of the lens system. The $`\mathrm{sin}2k\theta `$ terms might be included as well but they would not change the density on these axes, which is what we are interested in.

$$\beta F(\theta )=r_0^{2\beta }\left(1\underset{k=1}{\overset{\mathrm{}}{}}a_k\mathrm{cos}2k\theta \right)\text{.}$$
(37)
To be compatible with equation (20), the coefficients have to meet the condition

$$\underset{k=1}{\overset{\mathrm{}}{}}a_k=(2\beta )\epsilon \text{.}$$
(38)
We notice, that densities (6) can become negative near the axes. To minimize the angular density contrast, we have to keep only the monopole and quadrupole terms and set all higher coefficients to 0. The potential is then equivalent to a density of

$$\sigma =\frac{\beta }{2}\left(\frac{r}{r_0}\right)^{\beta 2}\left(1+\epsilon (2\beta )\left(\frac{4}{\beta ^2}1\right)\mathrm{cos}2\theta \right)\text{,}$$
(39)
which is everywhere positive only for sufficiently high values of $`\beta `$. When using realistic mass models, we can therefore expect a lower bound for $`\beta `$ to achieve acceptable fits. This applies not only to this special symmetric lens system but is true in general. Numerical models presented in the next section will confirm this result for $`\gamma =0`$ (see Figure 4).

### 13.2 The Einstein cross Q2237+0305

This lens is not usually taken into consideration when thinking about determination of $`H_0`$, because the time delays are expected to be very small and can therefore not be determined easily. Here we show that even if all three time delays were known exactly, constraints for the Hubble constant would still be very weak.

The degeneracy caused by the unknown mass index $`\beta `$ was already discussed by Wambsganß & Paczyński (1994) for spherical models plus external shear. The authors found the $`2\beta `$ scaling using numerical models. We now want to investigate how strong the assumption of a spherical main galaxy really influences the results.

No time delays are available for 2237+0305 and they may never be determined. We can nevertheless calculate the critical shear defined before and compare it with the a typical value one gets for spherical models. Positions including error bars used for this were taken from Crane et al. (1991) to make results comparable with Wambsganß & Paczyński (1994).

The critical shear as calculated from these positions<sup>3</sup><sup>3</sup>3All coordinates in this paper: $`x`$ to east and $`y`$ to north. is $`𝜸_\mathrm{c}=(0.0915\pm 0.044,0.0958\pm 0.045)`$ or $`|𝜸_\mathrm{c}|=0.13`$ where the errors are $`1\sigma `$ bounds from Monte Carlo simulations. Numerical modelling results in a shear of $`0.0696`$ almost exactly parallel to $`𝜸_\mathrm{c}`$ for isothermal spherically symmetric potentials. We therefore expect the time delays (or $`h`$ if we take $`\mathrm{\Delta }t`$ as known) of the spherical model with shear to be a factor $`0.4745`$ smaller than in the shearless case. For the moderately small shear of $`0.07`$, this is a huge effect. This factor is in good agreement with the expected value of $`1/2`$ for idealized systems.

To compare results in the general case, we performed numerical model fitting with an elliptical potential approach plus external shear.

$`F(\theta )`$ $`\left(\left({\displaystyle \frac{\mathrm{cos}\theta }{1+\epsilon }}\right)^2+\left({\displaystyle \frac{\mathrm{sin}\theta }{1\epsilon }}\right)^2\right)^{\frac{\beta }{2}}`$ (40)
Elliptical potentials are known to be unphysical for large ellipticities. Although $`\epsilon `$ is small in our case, we may expect unrealistic solutions for small values of $`\beta `$, because the limit of acceptable ellipticities vanishes for $`\beta 0`$ (cf. last section). In fact the fitted $`\epsilon `$ also increases with decreasing $`\beta `$.

We decided not to use invented time delays (calculated for a reference model) to fit the models. In this way, we can check the validity of our results even for cases where multiple time delays are not used as constraints. Plots of the residuals, ellipticity, shear and time delay between component A and B are shown in Figure 4.

The non-vanishing residuals at $`\beta 2`$ might be worrying at first, because a spherical model without external shear can fit any image configuration for $`\beta =2`$. The result would be a sheet of constant density equal to the critical density. There would be no isolated images but an area of constant (and very high) surface brightness. One might thus naively think, that residuals should be very small near $`\beta =2`$. This is not the case. It is true, that the deviations of the projected images will become arbitrary small in the source plane. On the other hand, however, the magnifications diverge in the limit, causing the deviations in the lens plane to stay finite.

Three families of models have to be discussed. First, we fixed the ellipticities at 0 to compare with the results from Wambsganß & Paczyński (1994). The residuals are almost constant for spherical models. The shear and the time delays scale very accurately with $`2\beta `$ as in the idealized considerations.

More interesting in our context is the behaviour in the shearless case where we expect to find a $`(2\beta )/\beta `$ scaling of the time delays. The $`\chi ^2`$ gets unacceptably large for $`\beta 1`$. In the other cases ($`\beta 0.6`$), the agreement with the theoretical predictions also shown in Fig. 4 is very good. For isothermal models, the ratio of time delays calculated for the two models (spherical$`/`$shearless) is 10–20 per cent larger than predicted by the critical shear. This is still a very good agreement considering that the time delays were not used as constraints for the numerical models.

Real lensing galaxies usually are elliptical and also embedded in an external shear field. We therefore also let both $`\epsilon `$ and $`\gamma `$ vary freely in order to minimize $`\chi ^2`$. Counting the formal number of constraints and parameters, we expect a minimum of $`\chi ^2=0`$. The fact, that $`\chi ^2`$ does not vanish is a confirmation of the degeneracy involving shear and ellipticity already discussed by Witt & Mao (1997). The effective number of parameters is therefore smaller then the formal one. Because we notice the ellipticity changing only slightly with $`\beta `$, a scaling like that in the spherical case is expected. This can indeed be seen in the plot, where the time delay scales with $`2\beta `$. Contrary to the spherical case, $`𝜸`$ does not scale proportional to $`2\beta `$ but is additionally shifted by a constant offset. This is due to the fact that part of the shear has been transformed to ellipticity.

We conclude, that both scaling laws for $`h`$ or the time delays can be relevant, depending on the family of models used. In the case of 2237+0305, however, the influence of external shear is stronger than the effect of $`\beta `$ for any realistic values of the latter.

### 13.3 PG 1115+080

Models and their degeneracies for this quadruple system have been studied extensively (e.g. Courbin et al., 1997; Impey et al., 1998; Keeton & Kochanek, 1997; Saha & Williams, 1997; Schechter et al., 1997; Zhao & Pronk, 2001), leading to a variety of more or less realistic mass distributions and a range of values for the Hubble constant between about 40 and $`80\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$.

All these authors agree on the importance of taking the effect of a nearby galaxy group into account. Our formalism includes this group as an external shear of the order $`\gamma 0.1`$. This shear is however not well constrained. Keeton et al. (1997) showed, that the residuals do not change much for ranges of $`\gamma 0.06`$$`0.2`$. With our general family of models, the effect of unknown external shear can be quantified by using the critical shear and equation (24). To calculate $`\gamma _\mathrm{c}`$, we only need image positions relative to an arbitrary reference centre. The uncertainty in the galaxy position, which is usually much higher than in the image positions, does not affect the result. Using the HST observations from Impey et al. (1998) with their claimed accuracy of $`0\stackrel{}{.}002`$ as a basis for Monte Carlo simulations, we obtain a critical shear of $`𝜸_\mathrm{c}=(0.142\pm 0.002,0.167\pm 0.003)`$. Although the ground based positions from Courbin et al. (1997) are not compatible with the HST results within the formal error bars, the critical shear from these data is about the same, $`𝜸_\mathrm{c}=(0.143\pm 0.003,0.160\pm 0.005)`$. This means, that any shear of the order 0.1 can change the results for $`H_0`$ significantly. As the uncertainties in $`\gamma `$ are of this order of magnitude, large effects on $`h`$ result even for fixed $`\beta `$.

The formalism we present in this article was not developed to directly determine $`H_0`$ from observations but to study the model degeneracies and scaling laws. One might nevertheless try to use the prescription from section 7 to obtain an estimate for the Hubble constant and the external shear for an isothermal model. The errors of the observational data of course have to be taken into account.

Time delays derived from the same light curves have been published by Schechter et al. (1997) and Barkana (1997). The A component was not resolved in these observations, and time delays were only determined relative to the sum of A<sub>1</sub> and A<sub>2</sub>. This is justified by the small time delay between the two, which is expected to be the order of hours. For our Monte Carlo simulations, we assumed a time delay between the A images of $`(0\pm 1)`$ days.

Results for $`H_0`$ in the isothermal case differ depending on which set of time delays and positions is used. With the redshifts $`z_\mathrm{s}=1.722`$ and $`z_\mathrm{d}=0.310`$ (Tonry, 1998) and an Einstein-de Sitter universe, we obtain values of $`H_0`$ between 47 and $`58\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ with errors between 12 and 30 per cent ($`1\sigma `$). The external shear is only weakly constrained, but seems to be of the order 0.1. We have to stress, that this result includes all possible isothermal models with arbitrary angular dependence and is thus much more general than elliptical mass distributions. As long as different determinations of the positions and time delays are not consistent with each other within their error bars, any results for $`H_0`$ have to be interpreted with care of course.

### 13.4 RX J0911.4+0551

This quad, initially discovered as a triple (Bade et al., 1997), has unique geometrical properties and is a strong candidate for time delay determination, although no result has been published yet. Rapid variability has been detected in the X-ray regime, providing the possibility of a determination of all three time delays with unprecedented accuracy (Chartas et al., 2001). Our first models (elliptical potential plus shear, $`\chi ^21`$) were presented in Burud et al. (1998). The external shear in the best-fitting model is $`\gamma =0.32`$ and points almost exactly in the direction of a nearby cluster of galaxies.

The redshift and velocity dispersion of this cluster was measured by Kneib et al. (2000) to $`z=0.769`$ and $`\sigma _v=(836\begin{array}{c}+180\\ 200\end{array})\mathrm{km}\mathrm{s}^1`$. From this, they derive an absolute shear of $`\stackrel{~}{\gamma }=0.11\begin{array}{c}+0.05\\ 0.04\end{array}`$ using a SIS model for the cluster.

To compare with lens models, we have to use the reduced shear $`\gamma =\stackrel{~}{\gamma }/(1\kappa )`$, because the convergence $`\kappa `$ caused by the cluster was not taken into account explicitly. The mass-sheet degeneracy simply scales all parameters in (22) with $`1/(1\kappa )`$. For SIS models, $`\stackrel{~}{\gamma }=\kappa `$ holds and we obtain $`\gamma =0.12\begin{array}{c}+0.05\\ 0.04\end{array}`$. This measurement is not in good agreement with the model from Burud et al. (1998). The reduced shear differs by about 0.2, the direction of the two being almost identical. A possible explanation for this discrepancy is the presence of a second galaxy close to the main lens that might change the potential considerably. It is also possible that the internal asymmetry of the main galaxy itself can not be described as an elliptical mass distribution.

To estimate the uncertainty in time delays or the Hubble constant derived from them, we use again the critical shear which is $`𝜸_\mathrm{c}=(0.553\pm 0.013,0.101\pm 0.005)`$ or $`|𝜸|=0.56`$. Even with this very large critical shear, the uncertainty in the real $`\gamma `$ has significant effects, because it is large as well.

### 13.5 B1608+656

This system is the first and up to now only quad for which all three independent time delays have been measured (Fassnacht et al., 1999). This offers the unique possibility to apply our method to a system providing the complete set of constraints. HST images show a main lensing galaxy but also a weaker second galaxy between the four images (Jackson et al., 1998). We nevertheless apply the method to B1608+656, if only to see if the effect of a secondary lens can be detected in this way, e.g. by pretending there is a very large external shear. Data for the positions of the images and the main lensing galaxy were taken from Koopmans & Fassnacht (1999). The formal accuracy of the image positions is extremely high, of the order 2–12 $`\mu `$arcsec. For the Monte Carlo simulations, we used 1 mas scatter in each coordinate to account for possible shifts by local density fluctuations, caused for example by globular clusters (Mao & Schneider, 1998).

We can use the general equation (22) to determine the Hubble constant for shearless models. With the redshifts of $`z_\mathrm{d}=0.6304`$ and $`z_\mathrm{s}=1.394`$ and standard Einstein-de Sitter cosmology, we obtain a value of $`H_0=(37\pm 5)\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ for $`\beta =1`$. The critical shear is $`𝜸_\mathrm{c}=(\mathrm{\hspace{0.17em}0.072}\pm 0.001,\mathrm{\hspace{0.17em}0.069}\pm 0.001)`$. For isothermal models, equation (22) predicts a shear of about $`𝜸=(0.32\pm 0.02,0.11\pm 0.01)`$, depending on which HST image is used to determine the galaxy position. The result for the Hubble constant is $`H_0=(130\pm 15)\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$.

The enormous differences in both $`H_0`$ predictions is a consequence of the large external shear $`\gamma \gamma _\mathrm{c}`$. No external shear at all is needed to fit the data when including the influence of the second galaxy in the field. The models in Koopmans & Fassnacht (1999) even predict velocity dispersions for both galaxies which are of the same order of magnitude. This is surprising, since the secondary galaxy is much weaker in all bands in the optical images.

## 14 Discussion

We used a very general semi-parametric lens model approach to study the changes of time delays and the determined Hubble constant with the assumed radial density slope, quantified by the radial mass index $`\beta `$. By using only linear constraints, it was possible to keep all the fundamental equations linear. This made the study of the $`\beta `$ dependence easy. For fixed external shear in quadruple lenses, this resulted in the simple scaling law $`H_0(2\beta )/\beta `$, independent of the lens geometry, the time delay ratios or the external shear. This means, that a systematic error in the assumed $`\beta `$ will have exactly the same effect on all lenses and will not show as scatter in the results. The good agreement between $`H_0`$ measurements from different lenses (Koopmans & Fassnacht, 1999) should therefore not be taken as evidence for an accurate determination of $`H_0`$. It merely shows, that all lensing galaxies seem to have more or less the same $`\beta `$.

In nearly isothermal models, a systematic error of only 10 per cent in $`\beta `$ will result in an error of about 20 per cent in the deduced Hubble constant. To compare the results from lensing not only with each other but also with results derived from other methods, this possible source of error has to be taken into account.

Furthermore, it is important not only to be aware of this effect, but to try and obtain better constraints on the radial mass profile. Possible ways to do this include modelling lenses with multiply imaged extended sources to obtain constraints for the lensing potential at a wide range of distances to the lens centre. Lenses with multiply imaged point sources can constrain the models only at a small number of radii. This should be done not only for systems with measured time delays, but for as many applicable lenses as possible to acquire reliable measurements of $`\beta `$ for a representative sample of lenses.

Another possibility is the detailed study of the dynamics of lensing galaxies or of other (usually low redshift) galaxies to learn more about the range of realistic galaxy mass profiles and use the results of this analysis to produce more realistic models for the lenses.

We also quantified the effect of external shear by introducing the concept of a ‘critical shear’ $`𝜸_\mathrm{c}`$. The effect of $`\gamma `$ on $`H_0`$ is linear and strongest in the direction of $`𝜸_\mathrm{c}`$. For a fixed direction, its amount is proportional to $`\gamma /\gamma _\mathrm{c}`$. The shear has to be exactly critical to fit the observations in point mass lenses. The value of $`𝜸_\mathrm{c}`$ can be found in a geometrical way. It is given by the ellipticity of the roundest ellipse passing through all images. For $`𝜸=𝜸_\mathrm{c}`$ and non-vanishing $`\beta `$, the time delays $`H_0\mathrm{\Delta }t`$ become zero. This is also true for a whole family of models which are represented by the less symmetric ellipses fitting the images. In the appendix, we discuss the relation to possible Einstein rings which are exactly given by these ellipses.

The effect of shear is also the clue in understanding the compatibility of the general scaling law with the simpler one of $`h_{\mathrm{sph}}2\beta `$ for spherical models. In the latter models, the shear is constrained by the observational data and changes $`h`$ by a factor of $`\beta /2`$ when compared to shearless models. In cases where spherical models are able to fit the data, the allowed range of $`\gamma `$ results in an uncertainty of $`h`$ always covering both models. This may in certain cases only apply for a limited range of $`\beta `$.

Interestingly, the value of the critical shear and the Hubble constant in the general model (with fixed shear for $`H_0`$) do not depend on the position of the lensing galaxy. This may be of use for systems where this position cannot be determined accurately.

When using more general models or less constraints than in our calculations, the scaling laws still apply. They are only valid for a subset of the possible models then and one would have to expect even larger uncertainties when using the whole set. This also applies for lenses with less than four images. More special models on the other hand, like parametrized elliptical power-law models, may be able to constrain the range of possible results much better. Nevertheless, the scaling laws still apply for the range of models that are compatible with the constraints. Even in these cases, our results may be used to determine $`H_0`$ without explicit modelling.

## Acknowledgments

It is a pleasure to thank Ester Piedipalumbo, Shude Mao and Sjur Refsdal for interesting discussions on the subject. The very constructive comments from the referee P. Saha greatly helped in understanding the geometrical properties of the critical shear.

## Appendix A Einstein rings and high image multiplicities

An interesting property of the general power-law models we used in the main part of this paper is the possibility to produce Einstein rings from point sources for arbitrary values of the external shear. For an Einstein ring parametrized by $`r=r(\theta )`$, all points on this ring must have the same light travel time to meet Fermat’s theorem. For this appendix, we set $`\kappa `$ and $`\gamma _y`$ to 0 for simplicity. For an arbitrary direction of the shear, we just have to replace $`\theta `$ by $`\theta \theta _\gamma `$. We also assume, that $`\gamma <1`$.

With $`t=0`$, the general equations (22) describe an ellipse, which is centred on the lens in the special case $`𝒛_𝐬=0`$ or $`\beta =1`$.

$$r=r_0(1+\gamma \mathrm{cos}2\theta )^{1/2}$$
(41)
The minor and major axes are $`r_0/\sqrt{1\pm \gamma }`$, compatible with equation (26). The potential is in this case an elliptical one

$$F(\theta )=\frac{1}{\beta }r_0^{2\beta }(1+\gamma \mathrm{cos}2\theta )^{\beta /2}\text{.}$$
(42)
With respect to the tangential caustic, the effects of ellipticity and shear cancel in these models and the caustic degenerates to a point. This is qualitatively different from elliptical mass distributions where the caustic is deformed and overlaps itself, producing areas of higher multiplicities (see below). With arbitrary $`𝒛_𝐬`$ and $`\beta `$, the centre of the ellipse is shifted to $`𝒛_0`$ with<sup>4</sup><sup>4</sup>4This equation is valid for arbitrary directions of $`𝜸`$.

$$\left(\begin{array}{cc}1+\gamma _x& \gamma _y\\ \gamma _y& 1\gamma _x\end{array}\right)𝒛_0=\frac{1\beta }{2\beta }𝒛_𝐬\text{.}$$
(43)
For $`\beta 1`$, this shift can take any value if $`𝒛_𝐬`$ is varied. To obtain a globally unique function $`F(\theta )`$, the centre of the lens has to be located inside of the ellipse.

Even for lens systems with four images, it is always possible to find an ellipse passing through all of them, which can act as an Einstein ring for the corresponding value of $`𝜸`$ given by the ellipticity. This does not mean, that we always see an Einstein ring for this special value of external shear, as $`F(\theta )`$ is not constrained for angles between the images.

Small deviations from the Einstein ring case can lead to an arbitrary number of images near the former elliptical ring. Special cases of these systems (singular isothermal ellipsoidal mass distributions with shear) with up to eight images were mentioned by Lopez et al. (1998) and Witt & Mao (2000) and discussed in detail by Keeton et al. (2000). Evans & Witt (2001) present results for shearless models with arbitrary $`F(\theta )`$. From (20) and (21) we obtain the following condition for an image of a source at $`𝒛_𝐬=0`$:

$$\frac{F^{}}{F}=\frac{\beta \gamma \mathrm{sin}2\theta }{1+\gamma \mathrm{cos}2\theta }$$
(44)
A global solution for this differential equation is given by (42) which leads to the elliptical ring we discussed before. For a number of discrete images, a more general solution is possible:

$$F(\theta )=f(\theta )(1+\gamma \mathrm{cos}2\theta )^{\beta /2}$$
(45)
At the positions of the images $`\theta _i`$, (44) has to be met, leading to the simple condition

$$f^{}(\theta _i)=0\text{.}$$
(46)
As $`f(\theta )`$ is an arbitrary function, we can easily construct systems with any number of images. The radial coordinates of the images can then be determined to be

$$r(\theta _i)=\left(\beta f(\theta _i)\right)^{1/(2\beta )}(1+\gamma \mathrm{cos}2\theta _i)^{1/2}\text{.}$$
(47)


format: text

text:
# On the Symmetry of Universal Finite-Size Scaling Functions in Anisotropic Systems

(March 6, 2024)

## Abstract

In this work a symmetry of universal finite-size scaling functions under a certain anisotropic scale transformation is postulated. This transformation connects the properties of a finite two-dimensional system at criticality with generalized aspect ratio $`\rho >1`$ to a system with $`\rho <1`$. The symmetry is formulated within a finite-size scaling theory, and expressions for several universal amplitude ratios are derived. The predictions are confirmed within the exactly solvable weakly anisotropic two-dimensional Ising model and are checked within the two-dimensional dipolar in-plane Ising model using Monte Carlo simulations. This model shows a strongly anisotropic phase transition with different correlation length exponents $`\nu _{\mathrm{𝗅𝗅}}\nu _{}`$ parallel and perpendicular to the spin axis.

The theory of universal finite-size scaling (UFSS) functions is a key concept in modern understanding of continuous phase transitions Privman and Fisher (1984); Privman (1990); Henkel and Schollwöck (2001). In particular, it is known that the UFSS functions of a rectangular two-dimensional (2D) system of size $`L_{\mathrm{𝗅𝗅}}\times L_{}`$ depend on the aspect ratio $`L_{\mathrm{𝗅𝗅}}/L_{}`$ Binder (1990). For instance, in *isotropic* systems the scaling function at criticality $`\overline{U}_\mathrm{c}`$ of the Binder cumulant $`U=1\frac{1}{3}m^4/m^2^2`$ Binder (1981), where $`m^n`$ is the $`n`$-th moment of the order parameter, is known to be a universal function $`\overline{U}_\mathrm{c}(L_{\mathrm{𝗅𝗅}}/L_{})`$ for a given boundary condition. This quantity has been investigated by several authors in the isotropic 2D Ising model with periodic boundary conditions Binder and Wang (1989); Kamieniarz and Blöte (1993), while the influence of other boundary conditions on $`\overline{U}_\mathrm{c}(L_{\mathrm{𝗅𝗅}}/L_{})`$ has recently been studied in Refs. Okabe et al. (1999); Kaneda and Okabe (2001).

In *weakly* *anisotropic* systems, where the couplings are anisotropic ($`J_{\mathrm{𝗅𝗅}}J_{}`$ in the 2D Ising case), the correlation length of the infinite system in direction $`\mu =\mathrm{𝗅𝗅},`$ becomes anisotropic and scales like $`\xi _\mu ^{(\mathrm{})}(t)\widehat{\xi }_\mu t^\nu `$ near criticality. ($`t=(TT_\mathrm{c})/T_\mathrm{c}`$ is the reduced temperature and we assume $`t>0`$ without loss of generality.) This leads to a correlation length amplitude ratio $`\widehat{\xi }_{\mathrm{𝗅𝗅}}/\widehat{\xi }_{}`$ different from unity. The UFSS functions then depend on this ratio, i. e. $`\overline{U}_\mathrm{c}=\overline{U}_\mathrm{c}(L_{\mathrm{𝗅𝗅}}/L_{},\widehat{\xi }_{\mathrm{𝗅𝗅}}/\widehat{\xi }_{})`$. However, isotropy can be restored asymptotically by an anisotropic scale transformation, where all lengths are rescaled with the corresponding correlation length amplitudes $`\widehat{\xi }_\mu `$ Landau and Swendsen (1984); Indekeu et al. (1986); Yurishchev (1997). Thus the UFSS functions depend on $`L_{\mathrm{𝗅𝗅}}/L_{}`$ and $`\widehat{\xi }_{\mathrm{𝗅𝗅}}/\widehat{\xi }_{}`$ only through the *reduced* aspect ratio $`(L_{\mathrm{𝗅𝗅}}/\widehat{\xi }_{\mathrm{𝗅𝗅}})/(L_{}/\widehat{\xi }_{})`$.

In *strongly* *anisotropic* systems both the amplitudes $`\widehat{\xi }_\mu `$ as well as the correlation length exponents $`\nu _\mu `$ are different and the correlation length in direction $`\mu `$ scales like

$$\xi _\mu ^{(\mathrm{})}(t)\widehat{\xi }_\mu t^{\nu _\mu }.$$
(1)
Examples for strongly anisotropic phase transitions are Lifshitz points Hornreich et al. (1975) as present in the anisotropic next nearest neighbor Ising (ANNNI) model Selke (1992); Diehl and Shpot (2000); Pleimling and Henkel (2001), or the non-equilibrium phase transition in the driven lattice gas model Katz et al. (1983); Schmittmann and Zia (1995). Furthermore, in dynamical systems one can identify the $`\mathrm{𝗅𝗅}`$–direction with time and the $``$–direction(s) with space Henkel (1999), which in most cases give strongly anisotropic behavior.

Using the same arguments as above we conclude that UFSS functions of strongly anisotropic systems depend on the *generalized* reduced aspect ratio (cf. Binder and Wang (1989))

$$\rho =L_{\mathrm{𝗅𝗅}}L_{}^\theta /r_\xi ,\text{with}r_\xi =\widehat{\xi }_{\mathrm{𝗅𝗅}}\widehat{\xi }_{}^\theta $$
(2)
being the *generalized* correlation length amplitude ratio, and with the anisotropy exponent $`\theta =\nu _{\mathrm{𝗅𝗅}}/\nu _{}`$ Henkel (1999). Up to now no attempts have been made to describe the dependency of UFSS functions like $`\overline{U}_\mathrm{c}(\rho )`$ on the shape $`\rho `$ of strongly anisotropic systems. In particular, it is not known if the anisotropy exponent $`\theta `$ can be determined from $`\overline{U}_\mathrm{c}(\rho )`$. This problem is addressed in this work.

Consider a 2D strongly anisotropic finite system with periodic boundary conditions. When the critical point of the infinite system is approached from temperatures $`t>0`$, the correlation lengths $`\xi _\mu `$ in the different directions $`\mu `$ are limited by the direction in which $`\xi _\mu ^{(\mathrm{})}`$ from Eq. (1) reaches the system boundary first Binder (1990). For a given volume $`N=L_{\mathrm{𝗅𝗅}}L_{}`$ we define an “optimal” shape $`L_{\mathrm{𝗅𝗅}}^{\mathrm{opt}}\times L_{}^{\mathrm{opt}}`$ at which both correlation lengths $`\xi _\mu ^{(\mathrm{})}`$ reaches the system boundary simultaneously, i. e.

$$L_\mu ^{\mathrm{opt}}:=\xi _\mu ^{(\mathrm{})}(t)$$
(3)
for some temperature $`t>0`$ (Fig. 1a). We immediately find using Eqs. (1, 2) that the optimal shape obeys $`\rho _{\mathrm{opt}}1`$ for all $`N`$, giving $`L_{\mathrm{𝗅𝗅}}^{\mathrm{opt}}=r_\xi (L_{}^{\mathrm{opt}})^\theta `$. A system of optimal shape should show the strongest critical fluctuations for a given volume $`N`$ as the critical correlation volume $`\xi _{\mathrm{𝗅𝗅},\mathrm{c}}\xi _{,\mathrm{c}}`$ spans the whole system.

At the optimal aspect ratio $`\rho =1`$ the correlations are limited by both directions $`\mathrm{𝗅𝗅}`$ and $``$ (Fig. 1a). If the system is enlarged by a factor $`b>1`$ in the $`\mathrm{𝗅𝗅}`$–direction (Fig. 1b), the correlation volume may relax into this direction but does not fill the whole system due to the limitation in $``$–direction. A similar situation with exchanged roles occurs if the system is enlarged by a factor $`b>1`$ in the $``$–direction (Fig. 1c). We now *assume* that systems (b) and (c) are similar in the scaling region $`L_\mu ^{\mathrm{opt}}\mathrm{}`$, i. e. that their correlation volumes are asymptotically equal.

Hence we can formulate a *symmetry* *hypothesis*: Consider a system with periodic boundary conditions and optimal aspect ratio $`\rho =1`$ at the critical point. If this system is enlarged by a factor $`b>1`$ in $`\mathrm{𝗅𝗅}`$–direction, it behaves asymptotically the same as if enlarged by *the same factor* $`b`$ in $``$–direction.

To formulate this hypothesis within a finite-size scaling theory, we consider a 2D strongly anisotropic system of size $`L_{\mathrm{𝗅𝗅}}\times L_{}`$ which fulfills the generalized hyperscaling relation $`2\alpha =\nu _{\mathrm{𝗅𝗅}}+\nu _{}`$ Binder and Wang (1989). For our purpose it is sufficient to focus on the critical point. The universal finite-size scaling *ansatz* Privman and Fisher (1984); Binder and Wang (1989); Privman (1990); Binder (1990); Henkel and Schollwöck (2001) for the singular part of the free energy density $`f_\mathrm{c}=F_{\mathrm{s},\mathrm{c}}/(Nk_\mathrm{B}T_\mathrm{c})`$ reads Hucht (2002)

$$f_\mathrm{c}(L_{\mathrm{𝗅𝗅}},L_{})\frac{b_{\mathrm{𝗅𝗅}}b_{}}{N}Y_\mathrm{c}(b_{\mathrm{𝗅𝗅}},b_{})$$
(4)
with the scaling variables $`b_\mu =\lambda ^{\nu _\mu }L_\mu /\widehat{\xi }_\mu `$, where $`\lambda `$ is a free scaling parameter. The scaling function $`Y_\mathrm{c}`$ is universal for a given boundary condition, all non-universal properties are contained in the metric factors $`\widehat{\xi }_\mu `$. These metric factors occur due to the usual requirement that the relevant lengths are $`L_\mu /\xi _\mu ^{(\mathrm{})}(t)`$ near criticality, and cannot be absorbed into $`\lambda `$ in contrast to isotropic systems. For the three systems in Fig. 1 we set $`\lambda =(L_\mu ^{\mathrm{opt}}/\widehat{\xi }_\mu )^{1/\nu _\mu }`$ to get

$`f_\mathrm{c}(L_{\mathrm{𝗅𝗅}}^{\mathrm{opt}},L_{}^{\mathrm{opt}})`$ $``$ $`{\displaystyle \frac{1}{N}}Y_\mathrm{c}(1,1)`$ (5a)
$`f_\mathrm{c}(bL_{\mathrm{𝗅𝗅}}^{\mathrm{opt}},L_{}^{\mathrm{opt}})`$ $``$ $`{\displaystyle \frac{b}{N}}Y_\mathrm{c}(b,1)`$ (5b)
$`f_\mathrm{c}(L_{\mathrm{𝗅𝗅}}^{\mathrm{opt}},bL_{}^{\mathrm{opt}})`$ $``$ $`{\displaystyle \frac{b}{N}}Y_\mathrm{c}(1,b).`$ (5c)
The proposed symmetry hypothesis states that for $`b>1`$ Eqs. (5b) and (5c) are asymptotically equal in the scaling region where $`L_\mu ^{\mathrm{opt}}`$ is large,

$$f_\mathrm{c}(bL_{\mathrm{𝗅𝗅}}^{\mathrm{opt}},L_{}^{\mathrm{opt}})\stackrel{b>1}{}f_\mathrm{c}(L_{\mathrm{𝗅𝗅}}^{\mathrm{opt}},bL_{}^{\mathrm{opt}}).$$
(6)
Hence the scaling function $`Y_\mathrm{c}`$ has the simple symmetry

$$Y_\mathrm{c}(b,1)\stackrel{b>1}{=}Y_\mathrm{c}(1,b).$$
(7)
To rewrite $`Y_\mathrm{c}`$ as function of the generalized aspect ratio $`\rho `$ (Eq. (2)) instead of the quantities $`b_\mu `$, we set $`b_{}=1`$ in system (c) and get, as then $`\lambda =(bL_{}^{\mathrm{opt}}/\widehat{\xi }_{})^{1/\nu _{}}`$,

$$f_\mathrm{c}(L_{\mathrm{𝗅𝗅}}^{\mathrm{opt}},bL_{}^{\mathrm{opt}})\frac{b^\theta }{N}Y_\mathrm{c}(b^\theta ,1).$$
(8)
Eqs. (5c) and (8) are identical and we conclude that $`bY_\mathrm{c}(1,b)=b^\theta Y_\mathrm{c}(b^\theta ,1)`$. At this point it is convenient to define the scaling function $`\overline{Y}_\mathrm{c}(b)=bY_\mathrm{c}(b,1)`$ which fulfills

$$f_\mathrm{c}(L_{\mathrm{𝗅𝗅}},L_{})\frac{1}{N}\overline{Y}_\mathrm{c}(\rho ).$$
(9)
For this scaling function the symmetry reads

$$\overline{Y}_\mathrm{c}(\rho )\stackrel{\rho >1}{=}\overline{Y}_\mathrm{c}(\rho ^\theta ).$$
(10)
We see from Eq. (9) that the critical free energy density $`f_\mathrm{c}`$ is a universal function of the reduced aspect ratio $`\rho =L_{\mathrm{𝗅𝗅}}L_{}^\theta /r_\xi `$ without any non-universal prefactor, and that at criticality *all* system specific properties are contained in the non-universal ratio $`r_\xi `$ from Eq. (2).

Ansatz Eq. (4) can also be made for the inverse spin-spin correlation length at criticality Hucht (2002)

$$\xi _{\mu ,\mathrm{c}}^1(L_{\mathrm{𝗅𝗅}},L_{})\frac{b_\mu }{L_\mu }X_{\mu ,\mathrm{c}}(b_{\mathrm{𝗅𝗅}},b_{}).$$
(11)
The proposed symmetry gives $`X_{\mu ,\mathrm{c}}(b,1)\stackrel{b>1}{=}X_{\overline{\mu },\mathrm{c}}(1,b)`$, where $`\overline{\mu }`$ denotes the direction perpendicular to $`\mu `$. Hence the correlation volumes $`\xi _{\mathrm{𝗅𝗅},\mathrm{c}}\xi _{,\mathrm{c}}`$ of system (b) and (c) in Fig. 1 are indeed equal as assumed above and become $`\xi _{\mathrm{𝗅𝗅},\mathrm{c}}\xi _{,\mathrm{c}}\frac{N}{b}X_{\mathrm{𝗅𝗅},\mathrm{c}}^1(b,1)X_{,\mathrm{c}}^1(b,1)`$.

The correlation length amplitudes $`A_\xi ^\mu `$ in cylindrical geometry ($`b_\mu \mathrm{}`$, $`b_{\overline{\mu }}=1`$), which can be calculated exactly for many isotropic two-dimensional models within the theory of conformal invariance Cardy (1984) generalize to the strongly anisotropic form Henkel and Schollwöck (2001)

$$A_\xi ^\mu =\underset{L_{\overline{\mu }}\mathrm{}}{lim}L_{\overline{\mu }}^{\nu _\mu /\nu _{\overline{\mu }}}\underset{L_\mu \mathrm{}}{lim}\xi _{\mu ,\mathrm{c}}(L_{\mathrm{𝗅𝗅}},L_{}).$$
(12)
Inserting Eq. (11) they become

$$A_\xi ^{\mathrm{𝗅𝗅}}=r_\xi X_{\mathrm{𝗅𝗅},\mathrm{c}}^1(\mathrm{},1),A_\xi ^{}=r_\xi ^{1/\theta }X_{,\mathrm{c}}^1(1,\mathrm{})$$
(13)
which shows that in general $`A_\xi ^\mu `$ ist not universal. The symmetry hypothesis states that both limits of the scaling function $`X_{\mu ,\mathrm{c}}`$ are equal. Denoting this universal limit $`A_\xi :=X_{\mathrm{𝗅𝗅},\mathrm{c}}^1(\mathrm{},1)=X_{,\mathrm{c}}^1(1,\mathrm{})`$ we obtain $`A_\xi ^{\mathrm{𝗅𝗅}}=r_\xi A_\xi `$ and $`A_\xi ^{}=r_\xi ^{1/\theta }A_\xi `$ as well as the amplitude relations

$$A_\xi ^{1+\theta }=A_\xi ^{\mathrm{𝗅𝗅}}(A_\xi ^{})^\theta ,\frac{A_\xi ^{\mathrm{𝗅𝗅}}}{A_\xi ^{}}=r_\xi ^{1+1/\theta }.$$
(14)
These predictions can be checked within the exactly solved weakly anisotropic 2D Ising model with different couplings $`J_{\mathrm{𝗅𝗅}}`$ and $`J_{}`$, where the paramagnetic correlation length reads $`\xi _\mu ^{(\mathrm{})}(t)=(\mathrm{log}\mathrm{coth}(\beta J_\mu )2\beta J_{\overline{\mu }})^1`$ with $`\beta =1/k_\mathrm{B}T`$ Onsager (1944). The amplitude ratio $`r_\xi `$ at the critical point $`\mathrm{sinh}(2\beta _\mathrm{c}J_{\mathrm{𝗅𝗅}})\mathrm{sinh}(2\beta _\mathrm{c}J_{})=1`$ Onsager (1944) becomes $`r_\xi =\mathrm{sinh}(2\beta _\mathrm{c}J_{\mathrm{𝗅𝗅}})`$ Zia and Avron (1982). On the other hand, the inverse correlation length amplitudes in cylinder geometry Eq. (12) has been calculated Nightingale and Blöte (1983) to give $`A_\xi ^\mu =\frac{4}{\pi }\mathrm{sinh}(2\beta _\mathrm{c}J_\mu ),`$ which immediately yields Eqs. (13) if we insert the well known universal value $`A_\xi =4/\pi `$ Luck (1982); Cardy (1984). The left relation of Eqs. (14) has already been derived for several weakly anisotropic models, where it simplifies to $`A_\xi ^2=A_\xi ^{\mathrm{𝗅𝗅}}A_\xi ^{}`$ (Nightingale and Blöte, 1983, Eq. (7)).

To check the symmetry numerically in strongly anisotropic systems, we now focus on the Binder cumulant $`U`$. The scaling *ansatz* at criticality Eq. (4) becomes

$$U_\mathrm{c}(L_{\mathrm{𝗅𝗅}},L_{})\frac{1}{b_{\mathrm{𝗅𝗅}}b_{}}\stackrel{~}{U}_\mathrm{c}(b_{\mathrm{𝗅𝗅}},b_{})=\overline{U}_\mathrm{c}(\rho )$$
(15)
with the scaling function $`\overline{U}_\mathrm{c}(b)=\stackrel{~}{U}_\mathrm{c}(b,1)/b`$, and the calculation is completely analogous to the free energy case. The symmetry hypothesis for the cumulant scaling functions $`\stackrel{~}{U}_\mathrm{c}`$ and $`\overline{U}_\mathrm{c}`$ thus reads (cf. Eqs. (7,10))

$$\stackrel{~}{U}_\mathrm{c}(b,1)\stackrel{b>1}{=}\stackrel{~}{U}_\mathrm{c}(1,b),\overline{U}_\mathrm{c}(\rho )\stackrel{\rho >1}{=}\overline{U}_\mathrm{c}(\rho ^\theta ).$$
(16)
The generalization of the cumulant amplitude $`A_U^\mu `$ Binder (1981); Burkhardt and Derrida (1985) to strongly anisotropic systems is similar to Eq. (12) and gives

$$A_U^\mu =\underset{L_{\overline{\mu }}\mathrm{}}{lim}L_{\overline{\mu }}^{\nu _\mu /\nu _{\overline{\mu }}}\underset{L_\mu \mathrm{}}{lim}L_\mu U_\mathrm{c}(L_{\mathrm{𝗅𝗅}},L_{}).$$
(17)
Inserting the scaling *ansatz* Eq. (15) we now find

$$A_U^{\mathrm{𝗅𝗅}}=r_\xi \stackrel{~}{U}_\mathrm{c}(\mathrm{},1),A_U^{}=r_\xi ^{1/\theta }\stackrel{~}{U}_\mathrm{c}(1,\mathrm{}),$$
(18)
which again are in general not universal. Using the symmetry hypothesis we can define $`A_U:=\stackrel{~}{U}_\mathrm{c}(\mathrm{},1)=\stackrel{~}{U}_\mathrm{c}(1,\mathrm{})`$ and get $`A_U^{\mathrm{𝗅𝗅}}=r_\xi A_U`$, $`A_U^{}=r_\xi ^{1/\theta }A_U`$ as well as the identities (cf. Eqs. (14))

$$A_U^{1+\theta }=A_U^{\mathrm{𝗅𝗅}}(A_U^{})^\theta ,\frac{A_U^{\mathrm{𝗅𝗅}}}{A_U^{}}=r_\xi ^{1+1/\theta }.$$
(19)
The cumulant scaling function $`\overline{U}_\mathrm{c}(\rho )`$ must be extremal at $`\rho =1`$ due to symmetry. Furthermore, as a deviation from the optimal aspect ratio $`\rho =1`$ reduces the cumulant, it has a maximum at this point Binder and Wang (1989). A sketch of $`\overline{U}_\mathrm{c}(\rho )`$ for an assumed anisotropy exponent $`\theta =2`$ is depicted in Fig. 2. For $`\rho >1`$ both $`\overline{U}_\mathrm{c}(\rho )`$ and $`\overline{U}_\mathrm{c}(\rho ^{}=\rho ^\theta )`$ collapse onto a single curve, reflecting the proposed symmetry. It is obvious from Fig. 2 that $`\overline{U}_\mathrm{c}(\rho )`$ (and thus also $`\overline{Y}_\mathrm{c}(\rho )`$ from Eq. (10)) can not be analytic at $`\rho =1`$ in strongly anisotropic systems, as the two branches $`\overline{U}_\mathrm{c}(\rho )`$ and $`\overline{U}_\mathrm{c}(\rho ^{})`$ identical for $`\rho >1`$ fork at $`\rho =1`$ Hucht (2002). On the other hand, $`\overline{Y}_\mathrm{c}(\rho )`$ and $`\overline{U}_\mathrm{c}(\rho )`$ can be analytic at $`\rho =1`$ if the anisotropy exponent $`\theta =1`$, as in the case of the isotropic 2D Ising model (Ferdinand and Fisher, 1969, Eq. 3.37).

To check the symmetry hypothesis in a strongly anisotropic system, I performed Monte Carlo simulations of the two-dimensional dipolar in-plane Ising model Hucht (2002)

$$=\frac{J}{2}\underset{ij}{}\sigma _i\sigma _j+\frac{\omega }{2}\underset{ij}{}\frac{(r_{ij}^{})^22(r_{ij}^{\mathrm{𝗅𝗅}})^2}{|\stackrel{}{r}_{ij}|^5}\sigma _i\sigma _j$$
(20)
with spin variables $`\sigma =\pm 1`$, ferromagnetic nearest neighbor exchange interaction $`J>0`$, and dipole interaction $`\omega >0`$. The distance $`\stackrel{}{r}_{ij}=(r_{ij}^{\mathrm{𝗅𝗅}},r_{ij}^{})`$ between spin $`\sigma _i`$ and $`\sigma _j`$ is decomposed into contributions parallel and perpendicular to the spin axis. In the simulations the Wolff cluster algorithm Wolff (1989) for long range systems proposed by Luijten and Blöte Luijten and Blöte (1995) was used, modified to anisotropic interactions. In contrast to earlier work De’Bell and Geldart (1989); Bulenda et al. (2000) using renormalization group technics it is found that this model shows a strongly anisotropic phase transition. The details of the simulations will be published elsewhere Hucht (2002).

After $`T_\mathrm{c}`$ was determined, systems with constant volume $`N=L_{\mathrm{𝗅𝗅}}L_{}`$ were simulated, which was chosen to have a large number of divisors in order to get many different aspect ratios (e. g. $`N=2^63^35^2=43200`$ has $`84`$ divisors). The resulting critical cumulant $`U_\mathrm{c}(L_{\mathrm{𝗅𝗅}}L_{}^\theta )`$ for two different volumes $`N=4320,43200`$ is depicted in the inset of Fig. 3. As expected, both curves have the same maximum value $`\overline{U}_\mathrm{c}(1)=0.555(5)`$ at criticality. With variation of $`\theta `$ the curves are shifted horizontally and collapse for $`\theta =2.1(3)`$, with maximum at $`r_\xi =0.415(40)`$. To check the proposed symmetry we fold the left branch with $`\rho <1`$ (open symbols) to the right and rescale the $`\rho `$-axis with $`\theta `$. The resulting data collapse for $`\rho >1`$ is shown in Fig. 3. This collapse and the additional condition that both curves must go to zero as $`A_U/\rho `$ allows a precise determination of $`\theta `$ and $`r_\xi `$ as well as of the universal amplitude $`A_U=3.5(2)`$.

In conclusion, I postulate a symmetry of universal finite-size scaling functions under a certain anisotropic scale transformation and generalize the Privman-Fisher equations Privman and Fisher (1984) to strongly anisotropic phase transitions on rectangular lattices at criticality. It turns out that for a given boundary condition the only relevant variable is the generalized reduced aspect ratio $`\rho =L_{\mathrm{𝗅𝗅}}L_{}^\theta /r_\xi `$ and that e. g. the free energy scaling function Eq. (9) obeys the symmetry $`\overline{Y}_\mathrm{c}(\rho )\stackrel{\rho >1}{=}\overline{Y}_\mathrm{c}(\rho ^\theta )`$. At criticality, the free energy density $`f_\mathrm{c}`$, the inverse correlation lengths $`\xi _{\mu ,\mathrm{c}}`$, and the Binder cumulant $`U_\mathrm{c}`$ are universal functions of $`\rho `$, without a non-universal prefactor. All system specific properties are contained in the non-universal correlation length amplitude ratio $`r_\xi `$ (Eq. (2)).

The generalization to higher dimensions is straightforward Hucht (2002), an interesting application would be the precise determination of the exponent $`\theta `$ at the Lifshitz point of the three-dimensional ANNNI model Diehl and Shpot (2000); Pleimling and Henkel (2001). An open question is the validity of the proposed symmetry in non-equilibrium systems with appropriate boundary conditions, which recently have been shown to exhibit Privman-Fisher universality Henkel and Schollwöck (2001).

###### Acknowledgements.

I thank Sven Lübeck and Erik Luijten for valuable discussions and Malte Henkel for a critical reading of the manuscript. This work was supported by the Deutsche Forschungsgemeinschaft through SFB 491.


text:
# The Pairwise Velocity Distribution Function of Galaxies in the LCRS, 2dF, and SDSS Redshift Surveys

## 1 Introduction

If, as generally accepted, clustering evolves through gravitational instability then the properties of galaxy clustering and the distribution of galaxy peculiar velocities contain fundamental information on the global properties of the Universe and provide important constraints on cosmological models. For a given initial mass fluctuation spectrum, clustering evolution and the generation of galaxy peculiar velocities is strongly dependent on the value of the mass density parameter $`\mathrm{\Omega }_m`$. For example, in a low $`\mathrm{\Omega }_m`$ Universe mostly devoid of matter, it is very difficult to generate high peculiar velocities since the Hubble flow will dominate over gravitational instability from early times unless the mass fluctuation amplitude is very high. Consequently, low values of $`\mathrm{\Omega }_m`$ generate relatively smaller peculiar velocities than high values for a given fluctuation amplitude. Once both the present galaxy clustering amplitude and peculiar velocity distribution functions are well-determined, cosmologists will have made significant progress in determining the initial conditions and evolution of the mass density in the Universe.

Unfortunately, neither galaxy positions nor peculiar velocities are presently measurable with significant accuracy to make robust, direct measurements of these two quantities. Rather, most of the information available today is contained in galaxy redshift surveys. As is well-known, the redshift to a galaxy represents the sum of the galaxy’s ‘distance’ and its radial peculiar velocity. Therefore, measurements of the galaxy two-point correlation function are contaminated by these redshift distortions and what is measured is appropriately called the redshift-space correlation function (see also Juszkiewicz et al. 2000).

It is possible, however, to exploit the anisotropies created by the peculiar velocities in the redshift-space correlation function (see Peebles 1980). In this method, the redshift space correlation function is represented in two-dimensions, where the axes correspond to the directions parallel ($`\pi `$) and perpendicular ($`r_p`$) to the line-of-sight. The resultant correlation function $`\xi _z(r_p,\pi )`$ is then anisotropic since the the peculiar velocities distort the correlation function principally along the line-of-sight. Since the galaxy-galaxy correlation function is a two-point estimator, the anisotropies being generated actually reflect the value of the pairwise galaxy peculiar velocity distribution function. Measurement of this distribution function is the focus of this Letter.

In a seminal paper, Davis and Peebles (1983) using redshift information from the first Center for Astrophysics Redshift Survey (CfA1) measured a value for the second moment of the pairwise velocity distribution function (hereafter PVD) of $`\sigma _{12}=340\pm 40`$ $`\mathrm{km}\mathrm{sec}^1`$ and characterized the distribution as an exponential. Subsequent numerical work predicted a much larger value of approximately 1000 $`\mathrm{km}\mathrm{sec}^1`$ for a standard $`\mathrm{\Omega }h=0.5`$ Cold Dark Matter (CDM) model (Davis et al. 1985). However, later work has questioned the accuracy of the original CfA1 result (Mo, Jing, & Börner 1993, Somerville, Davis & Primack 1997).

Subsequent redshift surveys have given somewhat discrepant results. The IRAS 1.2 Jy survey (Fisher et al. 1994) was in good agreement with a value of $`\sigma _{12}=317_{49}^{+40}`$ $`\mathrm{km}\mathrm{sec}^1`$ while Marzke et al. (1995) using the second Center for Astrophysics Redshift Survey (CfA2) together with the Southern Sky Redshift Survey (SSRS) found $`\sigma _{12}=540\pm 180`$ $`\mathrm{km}\mathrm{sec}^1`$. Similar analysis applied to the Las Campanas Redshift Survey, hereafter LCRS, has found $`\sigma _{12}=452\pm 60`$ $`\mathrm{km}\mathrm{sec}^1`$(Lin et al. 1995). Guzzo et al. (1997) find $`\sigma _{12}=345_{65}^{+95}`$ $`\mathrm{km}\mathrm{sec}^1`$ for late-type galaxies using the Pisces-Perseus Redshift Survey and report this to be a fair estimate for field galaxies. Other recent measurements include Small et al. (1999) using the Two Norris Redshift Surveys who reports $`326_{52}^{+67}`$ $`\mathrm{km}\mathrm{sec}^1`$ and Ratcliffe et al. (1998) using the Durham/UKST Galaxy Redshift Survey with an estimate of $`416\pm 36`$ $`\mathrm{km}\mathrm{sec}^1`$.

Numerical work and re-analysis of existing surveys have shown the sensitivity of the measurement of the pairwise velocity dispersion to the idiosyncracies of the data using standard techniques. Mo et al. (1993) found that the estimated dispersion is extremely sensitive to the presence of rich clusters in a sample, $`\sigma _{12}=`$ 300 to 1000 $`\mathrm{km}\mathrm{sec}^1`$ for subsets of the same data. Zurek et al. (1994) using high resolution CDM simulations found large variations in the value of the dispersion on CfA1 size scales. Using real data, Marzke et al. (1995) found that by excluding the rich clusters from their survey, R$`1`$, the measured velocity dispersion dropped to $`\sigma _{12}=295\pm 99`$ $`\mathrm{km}\mathrm{sec}^1`$.

The sensitivity of the standard methods to the presence of rich clusters is primarily due to the fact that they generally only estimate the second moment of the distribution, which is highly sensitive to the existence of hot, rich clusters in the data. Recognizing this problem, several authors have invented new statistics that are less sensitive to rich clusters (see Kepner, Summers & Strauss 1997, Davis, Miller & White 1997, Strauss, Ostriker & Cen 1998, and Baker, Davis & Lin 2000). While new statistics can be a powerful approach, their measures are often somewhat indirect and their results less intuitive.

To circumvent this problem and with the intent to determine the entire peculiar velocity distribution function rather than just the second moment, Landy, Szalay & Broadhurst (1998), hereafter LSB98, developed a method based on a Fourier-space deconvolution of the redshift-space distortions in the correlation function in the $`(r_p,\pi )`$ basis. Using this approach with the LCRS data, a value of $`363\pm 44`$ $`\mathrm{km}\mathrm{sec}^1`$ was obtained. Further, by recovering the Fourier transform of entire distribution function, it was directly shown that the PVD is well-characterized by an exponential.

Other recent measurements of the PVD using the same LCRS data but different techniques have also been reported. Jing, Mo, & Börner (1998) measured the PVD of the LCRS data and obtained $`570\pm 80`$ $`\mathrm{km}\mathrm{sec}^1`$. The discrepancy with the result of LSB98 was ascribed to a failure of LSB98 to account for infall effects (see Jing & Börner 1998). Another analysis of the LCRS data was presented in Jing & Börner (2001), using a method closely related to that of LSB98, with a PVD value of $`510\pm 70`$ $`\mathrm{km}\mathrm{sec}^1`$. The discrepancies between these two latter results and those in this Letter will be presented in the discussion.

Recently, Peacock et al. (2001) have published a measurement of the PVD using a larger sample of the 2dF survey than publicly released as an adjunct to their estimation of the cosmological parameter $`\beta `$. They obtain a value of $`385`$ $`\mathrm{km}\mathrm{sec}^1`$ but do not report standard errors.

## 2 Data

The three data sets used in this analysis are all publicly available. The LCRS data is described in Shectman et al. 1996. The results for that data published here are identical to those reported in LSB98 and the reader is directed there for further information.

The 2dF data is the release of June 2001 (see www.mso.anu.edu.au/2dFGRS/). This data has been cut to $`0.05<z<0.2`$ (150 to 570 $`h^1`$Mpc assuming $`\mathrm{\Omega }_m=0.3,\mathrm{\Omega }_\lambda =0.7`$) and is naturally broken up into two subsets toward the north and south galactic caps. These sets are disjoint on the scales under consideration. The north and south data consisted of 34805 and 42309 galaxies respectively. Random catalogs for this data were generated using the publicly available routines supplied by the 2dF Survey.

The SDSS data is the release of June 2001 (see Stoughton et al. 2001). As with the 2dF data, the redshift limits were set at $`0.05<z<0.2`$ and subsets of the data towards the north and south galactic caps were analyzed independently. This data consisted of 8918 galaxy positions in the North and 7862 in the South. The random catalogs for this data were kindly supplied by Adrian Pope.

## 3 Method

In essence, this approach extracts the Fourier transform of the PVD from the Fourier transform of the galaxy-galaxy correlation function in the $`(r_p,\pi )`$ basis using a deconvolution procedure. By working in Fourier space, the deconvolution is one of simple division that directly returns the Fourier transform of the PVD. Since the method here is identical to that described in LSB98, the reader is directed there for further details. However, given the stability of the reported results, it is worth re-iterating two of the principle advantages of working in Fourier space.

Firstly, by utilizing the Fourier transform of the $`\xi (r_p,\pi )`$ correlation function, the measurement of the PVD is no longer a pair-weighted statistic. This reduces the contamination by rich clusters in the data. More specifically, the presence of rich clusters in the sample super-imposes a number of pairs with a higher dispersion than the thermal dispersion of field galaxies. These pairs predominately effect the tail of the correlation function and produce an overestimation of the PVD, especially of its second moment. By determining the Fourier transform of the PVD rather than the PVD itself, the large number of bins in the tail are all compressed into the innermost few resolution elements in Fourier space. Therefore they are naturally down-weighted in the fitting procedure, which uses equal weights for every cell in $`k`$-space. This method does not exclude the signal from the rich clusters but rather incorporates them in a way which decreases their contribution to the variance in a robust fashion.

Secondly, in the deconvolution procedure, the weighting function is effectively proportional to $`r_p^1`$, which emphasizes the high signal-to-noise core of the correlation function. Since most of the mass of the correlation function lies within the central core, the resulting signal will predominately reflect the value of the PVD within a scale of about 1$`h^1`$Mpc. This both reduces contamination by infall effects, which are expected to be small on these scales, and localizes the measurement to reflect the value of the PVD on small scales.

## 4 Results

The Fourier transform of the PVD for each redshift survey is shown in Figure 1. In all cases, this function is well-characterized by a Lorentzian. Since the Fourier transform of a Lorentzian is an exponential, the PVD itself is well-characterized by an exponential. For comparison, the best-fitting Gaussian distributions are also shown. All values reported correspond to the equivalent exponential value for a given fitted Lorentzian.

The result is very stable across all three redshift surveys. Since the 2dF is a b-band selected survey and the LCRS and SDSS are both r-band selected, it is expected that the 2dF Survey will contain a relatively greater number of field galaxies. Thus, the PVD for the 2dF survey should have a lower value as is found: $`357\pm 17`$ $`\mathrm{km}\mathrm{sec}^1`$ for the SDSS versus $`331\pm 19`$ $`\mathrm{km}\mathrm{sec}^1`$ for the 2dF. However, this difference is not of high significance given the standard errors for the data.

The Lorentzian in the case of the 2dF data shows additional structure around the central peak. This type of structure has also shown up in some of the analyses of subsets of the data. The effect is likely a consequence of the fact that the method necessarily constrains the value of the signal to unity at the origin. Adjusting the fit to exclude the structure in the central peak does not significantly change the measured value of the dispersion.

In order to get a better handle on the uncertainties in the results, the north and south data in the 2dF survey was also divided into two halves and an identical analysis performed. All results are reported in Table 1.

## 5 Discussion

As is evident in Figure 1, the PVD function for all cases is clearly well-characterized by an exponential. Furthermore, Table 1 shows that these measurements are very stable both within and between data sets.

As mentioned in the Introduction and Method sections, one of the major difficulties in the measurement of the PVD has been the contamination by rich clusters in the data. For example, Marke et al. (1995) found that by excluding rich clusters in their data that the signal fell from $`\sigma _{12}=540\pm 180`$ $`\mathrm{km}\mathrm{sec}^1`$ to $`\sigma _{12}=295\pm 99`$ $`\mathrm{km}\mathrm{sec}^1`$. Additionally two analyses of the LCRS data, one using a method based on the correlation function and one presented in this Letter, give $`\sigma _{12}=452\pm 60`$ $`\mathrm{km}\mathrm{sec}^1`$(Lin et al. 1995) and $`363\pm 44`$ $`\mathrm{km}\mathrm{sec}^1`$ (LSB98) respectively. It was also shown by Mo et al. (1993) that the estimated dispersion is extremely sensitive to the presence of rich clusters in a sample, to the order of errors of several hundred $`\mathrm{km}\mathrm{sec}^1`$for subsets of the same data.

These previous findings emphasize the importance of applying identical methods to subsets of data and across surveys. The results presented here reflect a very robust and stable technique. The small variance in the measured signal across and between data sets clearly indicate that the signal is not being dominated by rich clusters in the data.

This conclusion is also supported by the agreement between the results of the two r-band selected surveys, the LCRS and SDSS, both in terms of the value of the measured signal and the small standard errors. The standard errors of the 2DF b-band selected survey are also of a similar magnitude. Since this latter b-band survey is expected to contain relatively fewer late-type galaxies, it would be expected that its standard error would be substantially less than that of the r-band selected surveys whose signal should be more strongly contaminated by clusters in the data.

One criticism of this method (see Jing & Börner 1998) is that it does not adequately take into account infall effects and consequently underestimates the value of the PVD; compare their value of $`570\pm 80`$ $`\mathrm{km}\mathrm{sec}^1`$ for the LCRS data. In a complete characterization (see Peebles 1980), the galaxy peculiar velocity distribution function $`f(v_{12}|r)`$ is usually modeled as an exponential along one dimension with

$$f(v_{12}|r)=\frac{1}{\sqrt{2}\sigma _{12}}\mathrm{exp}\left(\frac{\sqrt{2}|v_{12}\overline{v_{12}}|}{\sigma _{12}}\right).$$
(1)
In our analysis, the infall parameter $`\overline{v_{12}}`$ is not included for several reasons. Firstly, given the excellent fit of the Fourier transform of the PVD with the Lorentzian, there is no evidence that the result is being contaminated by infall effects. Secondly, as was discussed in LSB98, since this method relies on a Fourier transform of the correlation function and most of the mass of the correlation function lies within a central core of about 1 $`h^1`$Mpc, the resulting signal will predominately reflect the value of the PVD within a scale of about 1$`h^1`$Mpc where the infall is expected to be small.

In later work, Jing & Börner (2001) re-analyze the LCRS and obtain a result of $`510\pm 70`$ $`\mathrm{km}\mathrm{sec}^1`$. This method very similar to that of LSB98 although the model includes a large-scale infall parameter in the linear regime characterized by $`\beta `$. In our method by truncating the correlation function at approximately 30$`h^1`$Mpc, the approximate scale of the transition between the linear and non-linear regimes, and windowing the distribution, the signal is limited to the non-linear regime. This is very important considering that the distortions of the redshift correlation function due to the small-scale PVD and the large-scale linear velocity infall are competing effects and so their errors are positively correlated. Additionally, although Jing & Börner (2001) incorporate a model for the infall effects of the $`\beta `$ distortions, they do include it in the fit but rather set it at a value of $`\beta =0.5`$.

The competing effects of these two distortions can be clearly seen in Fig. 4 of Peacock et al. (2001) where a joint maximum likelihood fit of the estimation of $`\sigma _{12}`$ and $`\beta `$ from a larger 2dF data set is presented. The best fit values are $`\sigma _{12}=385`$ $`\mathrm{km}\mathrm{sec}^1`$ and $`\beta =0.43\pm 0.07`$. The LSB98 value of $`363\pm 44`$ $`\mathrm{km}\mathrm{sec}^1`$ is within one standard error of this result and would imply a $`\beta `$ of $`0.4`$, while the measurement of Jing & Börner (2001) $`510\pm 70`$ $`\mathrm{km}\mathrm{sec}^1`$ implies a best-fitting value of $`\beta 0.6`$.

Davis, Miller & White (1997) also present a method for determining the thermal velocity dispersion of galaxies using a galaxy-weighted rather than the more standard pair-weighted measures. An analysis of this result on the LCRS data is reported in Baker, Davis & Lin (2000). Although it is difficult to directly compare their approach to the one presented here, the results are generally consistent.

## 6 Conclusion

Since galaxy distances and peculiar velocities can not be adequately directly measured at this point, every approach to measure the properties of the galaxy peculiar velocity distribution function is somewhat indirect. For example, almost every method attempts to measure the pairwise galaxy peculiar velocity distribution function rather the the single particle function itself. Also, it has been shown that due to the existence of hot clusters in any large redshift survey, that a direct fit to the correlation function itself is problematic.

The method utilized in this Letter has several advantages over other extant techniques. Firstly, it relies on a straightforward Fourier transform of the $`\xi (r_p,\pi )`$ correlation function, a mathematical operation that is well understood. Secondly, in extracting the Fourier transform of the PVD, the need to model the underlying $`\xi (r)`$ is obviated. Thirdly, by restricting the analysis to non-linear scales, the need to model for infall effects is eliminated. Fourthly, the method returns the Fourier transform of the entire PVD function and its functional form can be directly investigated. And lastly, the technique is very stable within and between independent surveys.

Using this method, it has been shown that the galaxy pairwise velocity distribution function can be well-represented by an exponential for the three largest existing redshift surveys. The two r-band selected surveys, the LCRS and the SDSS, give values of with $`363\pm 44`$ $`\mathrm{km}\mathrm{sec}^1`$ and $`357\pm 17`$ $`\mathrm{km}\mathrm{sec}^1`$ respectively. The b-band selected 2dF survey gives slightly lower results with $`331\pm 19`$ $`\mathrm{km}\mathrm{sec}^1`$and this difference is expected given the different selection effects of r-band and b-band selected surveys. Further analysis of subsets of the data has shown that the signal is not sensitive to idiosyncracies in the data sets such as the prevalence of rich clusters.

The author would like to acknowledge support from the Jeffress Memorial Trust and NSF Grant AST 99-00835 and would also like to thank the LCRS, 2dF, and SDSS surveys for making their data available to the astronomical community. The author also acknowledges useful discussions with Alexander S. Szalay and thanks Adrian Pope for random data.