text:
# Fitting Gravitational Lenses: Truth or Delusion

## 1 INTRODUCTION

The fitting of models to observational data has always been a major concern in strong gravitational lensing (e.g., Schechter 2000). In most cases, the lens is only composed of 2 (a “doublet”) or 4 (a “quadruplet”) point-like images. For the doublets, the observable constraints are the four relative coordinates of the images with respect to the lensing galaxy, together with the one flux ratio of the first image to the second. For the quadruplets, this becomes eight relative coordinates and three flux ratios. Only for a handful of lens systems are time delays available. It therefore happens very frequently that lens models have more degrees of freedom than the number of constraining observations. In fact, the situation is often even worse than we have just described. Sometimes the centre of the lensing galaxy itself cannot be reliably identified. Sometimes a lens system is complicated by the possible effects of nearby bright galaxies. Very often, the flux ratios are untrustworthy, either because of differential extinction in the optical bands (e.g., Falco et al. 1999) or because of microlensing in the optical and radio (e.g., Irwin et al. 1989; Koopmans & de Bruyn 2000) or because of scintillation, scatter-broadening and free-free absorption in the radio (e.g., Patnaik et al. 1992, Jones et al. 1996, Winn, Rusin & Kochanek 2003). It is therefore already clear that considerable caution must be exercised in interpreting the results of fits to gravitational lens systems.

The approach of modellers has been largely two-pronged. First, general non-parametric methods have been developed (e.g., Saha & Wiliams 1997). These have proved powerful in exploring the range of degeneracies in the lensing mass that can give rise to a particular image configuration. For example, in the case of PG 1115+080, Saha & Williams present three different models obeying the lensing constraints, but for which the lensing mass distribution resembles a face-on spiral, an edge-on disk or a flattened elliptical galaxy respectively. Second, many simple parametric models have been thoroughly explored, such as those based on elliptically stratified potentials (e.g., Witt 1996; Witt & Mao 1997) or densities (e.g., Kassiola & Kovner 1993; Muñoz, Kochanek & Keeton 2001). The advantage of this is that the models may already incorporate some of the known properties of nearby galaxies. However, this approach may also be dangerous because simple ansatze like elliptical potentials or surface mass densities can introduce unexpected properties. These properties may be so severe that, for example, the flux ratios may not be well fitted and wrong conclusions may be drawn. This can be seen most clearly in the case of elliptical potentials, in which there are strong constraints on the flux ratios (e.g., Witt & Mao 2000; Hunter & Evans 2001).

All this attains added significance in the light of recent claims of evidence for substructure from “anomalous flux ratios” in strong lensing. Mao & Schneider (1998) were the first to point out that substructure may be needed to explain the flux ratios in some cases. The instance that they selected, the quadruplet B 1422+231, has three highly magnified bright images and one much fainter image. Three highly magnified images occur generically near a cusp, and a Taylor expansion gives a universal relationship that the sum of the fluxes of the two outer images should equal the flux of the middle image. This is strongly violated in B 1422+231 leaving substructure on top of the smooth model as the believable culprit. This result is supported by the detailed models of B 1422+231 by Bradac et al. (2002) which find that the discrepancy requires substructure on the mass scale $`10^6M_{}`$. Recently, Dalal & Kochanek (2002) looked at seven predominantly radio quadruplets and argued that the flux ratios were anomalous by comparison with those expected for simple isothermal lenses with external shear. They claimed that the anomalous flux ratios implied the existence of $`2`$ per cent of the mass of the lensing galaxy in substructure. Metcalf & Zhao (2002) looked at five optical quadruplets and similarly claimed evidence for anomalous flux ratios on comparison with those expected for simple elliptical power-law potentials with shear. The problem with this procedure is, of course, that anomalous flux ratios may not be the result of substructure at all, but may simply reflect deficiencies in the modelling.

This motivates us to introduce in sections 2 and 3 a new approach for fitting which can incorporate the flux ratios at outset and which can permit an arbitrary azimuthal dependence for the surface density. Importantly, all the free parameters enter linearly into the models and so the lens and flux ratio equations can always be solved by straightforward matrix inversion. The models are only restricted by the fact that the surface mass density must be positive. This algorithm is used in section 4 to assess the evidence for anomalous flux ratios. Readers who are mainly interested in this application may skip section 3 entirely on the first reading. This section is rather mathematical and explicitly derives the linear equations for the parameters to be fitted.

## 2 SCALE-FREE MODELS

### 2.1 Potential and Surface Mass Density

Scale-free galaxy models with arbitrary power-law fall-off are widely used in galactic astronomy and dynamics (e.g., Toomre 1982; Evans 1994; Evans, Carollo & de Zeeuw 2000). The isophotes of a scale-free galaxy have the same shape at every radius and are completely described by a shape function $`G(\theta )`$, which depends only on the position angle $`\theta `$ with respect to the major axis. In such galaxies, the convergence $`\kappa `$ and deflection potential $`\varphi `$ are proportional to a power of $`r`$, namely

$$\kappa =\frac{1}{2}G(\theta )r^{\beta 2},\varphi =r^\beta F(\theta ).$$
(1)
Here, ($`r,\theta `$) are familiar polar coordinates in the lens plane. The power is restricted to lie in the range $`0<\beta <2`$ for galaxy models. For $`\beta <0`$, the mass at the centre of the galaxy does not converge. For $`\beta >2`$, the surface mass density increases with radius, which is unrealistic. The important $`\beta =1`$ case corresponds to an everywhere flat rotation curve.

In gravitational lensing, some simple scale-free potentials were already investigated by Kassiola & Kovner (1993). They were followed by Witt, Mao & Keeton (2000), who pointed out that the time delay is independent of the angular function $`F(\theta )`$ for the special case $`\beta =1`$. Many applications to this general model followed (see Evans & Witt 2001; Zhao & Pronk 2001; Wucknitz 2002; Cardone et al. 2002).

From Poisson’s Law, we have

$$G(\theta )=\beta ^2F(\theta )+F^{\prime \prime }(\theta ).$$
(2)
Given $`F(\theta )`$, it is straightforward to generate $`G(\theta )`$. However, it is also possible to solve this ordinary differential equation for $`F(\theta )`$ in terms of $`G(\theta )`$, as shown in Appendix A. This is useful, as it gives the scale-free potential corresponding to a given set of scale-free equidensity contours.

The critical curves can be given analytically. Using the determinant of the Jacobian, we obtain

$`detJ(r,\theta )`$ $`=`$ $`(\beta 1)r^{2\beta 4}[\beta FF^{\prime \prime }(\beta 1)F^2+\beta ^2F^2]`$ (3)
$`r^{\beta 2}[F^{\prime \prime }+\beta ^2F]+1=0.`$
Setting $`\rho =r^{\beta 2}`$ yields a quadratic equation in $`\rho `$ which can easily be solved. The two roots correspond to the radial and the tangential critical curves respectively. They can be mapped onto the radial and tangential caustics with the lens equation using the components of the deflection angle (35). For the flat rotation curve case ($`\beta =1`$), the equation for the critical curve becomes linear and thus the caustic network is readily established, as shown in Evans & Witt (2001).

### 2.2 Fourier Expansions

Since $`F(\theta )`$ and $`G(\theta )`$ are periodic with period $`2\pi `$, we can expand them as Fourier series (see e.g. Bronshtein & Semendyayev 1998). For convenience, we start with the potential $`F(\theta )`$ and write

$$F(\theta )=\frac{a_0}{2}+\underset{k=1}{\overset{\mathrm{}}{}}[a_k\mathrm{cos}(k\theta )+b_k\mathrm{sin}(k\theta )].$$
(4)
Using eq.(2), we obtain now for $`G(\theta )`$

$$G(\theta )=\frac{a_0\beta ^2}{2}+\underset{k=1}{\overset{\mathrm{}}{}}[a_k(\beta ^2k^2)\mathrm{cos}(k\theta )+b_k(\beta ^2k^2)\mathrm{sin}(k\theta )].$$
(5)
The coefficients in the Fourier series must obey the condition $`G(\theta )0`$ for all $`\theta `$ to obtain a positive surface mass density. After constructing a fit, this is easy to check a posteriori.

As an example, let us consider an elliptical surface mass density with a flat rotation curve ($`\beta =1`$), which has

$$G(\theta )=A(\mathrm{cos}^2\theta +q^2\mathrm{sin}^2\theta )^{1/2}.$$
(6)
Since $`G(\theta )=G(\theta )`$ is an even function, all $`b_k`$ must vanish ($`b_k=0`$). Since $`G(\theta )=G(\pi \theta )`$, all $`a_k`$ with $`k`$ equal to an odd integer must also vanish. So, the only terms that remain in the Fourier series are the even integer cosine coefficients (cf. also Appendix B). Figure 1 shows how few terms are really required for accurate approximation. For an E3 galaxy, a Fourier series truncated after the first three non-vanishing terms already gives an excellent approximation to the true equidensity contours. Even for a highly flattened E7 galaxy, only the first seven non-vanishing terms are needed. This demonstrates that the Fourier expansions can give good results with few terms.

The advantage of the Fourier expansion is that it makes it very easy to solve the lens equation. For this reason, we list the Fourier coefficients for the popular elliptic density and potential models with flat rotation curves in Appendix B (although we do not need to use them in the main body of this paper).

### 2.3 The Image Positions and Fluxes

The lens equation relates the position of the source ($`\xi ,\eta `$) to the positions of the image ($`x,y`$) through the derivative of the lensing potential (e.g., Schneider, Ehlers & Falco 1992)

$$\xi =x\frac{\varphi }{x},\eta =y\frac{\varphi }{y}.$$
(7)
Using polar coordinates in the image plane, the lens equation can be recast as

$`\xi `$ $`=`$ $`r\mathrm{cos}\theta r^{\beta 1}[\beta \mathrm{cos}\theta F(\theta )\mathrm{sin}\theta F^{}(\theta )],`$ (8)
$`\eta `$ $`=`$ $`r\mathrm{sin}\theta r^{\beta 1}[\beta \mathrm{sin}\theta F(\theta )+\mathrm{cos}\theta F^{}(\theta )],`$ (9)
where we can express $`F(\theta )`$ and $`F^{}(\theta )`$ by Fourier series. Let us assume that we have measured the image positions $`(r_{\mathrm{}},\theta _{\mathrm{}})`$ (in polar coordinates) of a gravitational lens and inserted them into the lens equation. Immediately, we notice that all unknown quantities (i.e. $`\xi `$, $`\eta `$, and all $`a_i`$ and $`b_i`$) enter linearly into the lens equation, except for $`\beta `$.

Let us now add in the constraint that the magnification ratios of the images are equal to the observed values. Suppose that the ratio of the magnification of the $`\mathrm{}`$th image to the $`k`$th is measured to be $`f_k\mathrm{}`$, then

$$f_k\mathrm{}=\frac{detJ(r_{\mathrm{}},\theta _{\mathrm{}})}{detJ(r_k,\theta _k)},$$
(10)
where $`detJ`$ is given in eq. (3). We observe that the coefficients $`a_i`$ and $`b_i`$ enter in simple mixed quadratic form into the equations. However, for the special case of a flat rotation curve ($`\beta =1`$), all unknown quantities enter linearly into the equation for the magnification relation.

In other words, in the astrophysically important flat rotation curve case, the Fourier coefficients describing the lensing galaxy are related to the observables by a simple matrix equation. This makes the problem of fitting to both the image positions and the flux ratios a straightforward matter of matrix inversion.

A number of recent investigations have exploited the well-known linearity of the lens equation to expand the potential in basis functions (e.g., Keeton 2001). Such methods can be a powerful way of exploring a large set of models, although the drawback is that some of the parameter space is often associated with the creation of additional images. The need to check for this can undermine the utility of starting with such an algorithm. It is worthwhile emphasising how our work differs from such expansions. First, the fundamental and new contribution of our paper is that the flux ratio equation is linear for a galaxy with a flat rotation curve. Therefore, in our algorithm both the lens equations and the flux ratio equations are simultaneously solved. Second, the caustic network is straightforwardly available using the work in Section 2.1, and so it is very easy to check for the creation of unwanted additional images.

## 3 Fitting Gravitational Lenses

### 3.1 A Singular Value Decomposition Problem

First, we insert the Fourier series of $`F(\theta )`$ into eqs. (8)–(9) to establish the linear equations for fitting for $`\beta =1`$. Introducing a useful notation, we write for the lens equation

$`\xi (r,\theta )`$ $`=`$ $`r\mathrm{cos}\theta {\displaystyle \frac{a_0}{2}}\alpha _0(\theta ){\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}[a_k\alpha _k(\theta )+b_k\beta _k(\theta )],`$ (11)
$`\eta (r,\theta )`$ $`=`$ $`r\mathrm{sin}\theta {\displaystyle \frac{a_0}{2}}\widehat{\alpha }_0(\theta ){\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}[a_k\widehat{\alpha }_k(\theta )+b_k\widehat{\beta }_k(\theta )],`$ (12)
where $`a_k`$ and $`b_k`$ are the Fourier coefficients of the unknown function $`F(\theta )`$, which describes the angular part of the potential. For the coefficients $`\alpha _k`$ and $`\beta _k`$, we find

$`\alpha _k(\theta )`$ $`=`$ $`\mathrm{cos}\theta \mathrm{cos}(k\theta )+k\mathrm{sin}\theta \mathrm{sin}(k\theta ),`$
$`\beta _k(\theta )`$ $`=`$ $`\mathrm{cos}\theta \mathrm{sin}(k\theta )k\mathrm{sin}\theta \mathrm{cos}(k\theta ),`$ (13)
and for $`\widehat{\alpha }_k`$ and $`\widehat{\beta }_k`$, we find

$`\widehat{\alpha }_k(\theta )`$ $`=`$ $`\mathrm{sin}\theta \mathrm{cos}(k\theta )k\mathrm{cos}\theta \mathrm{sin}(k\theta ),`$
$`\widehat{\beta }_k(\theta )`$ $`=`$ $`\mathrm{sin}\theta \mathrm{sin}(k\theta )+k\mathrm{cos}\theta \mathrm{cos}(k\theta ),`$ (14)
for all $`k0`$. The lens equations (8)–(9) must hold at the positions of the $`n`$ images, given in polar coordinates as ($`r_{\mathrm{}},\theta _{\mathrm{}}`$) with $`\mathrm{}=1,\mathrm{},n`$. This provides $`2n`$ constraints on the unknown source position $`(\xi ,\eta )`$ and the unknown Fourier coefficients.

Now let us consider the magnification of the $`\mathrm{}`$th image. For the flat rotation curve case ($`\beta =1`$), we obtain

$$detJ=\frac{1}{p_{\mathrm{}}\mu _{\mathrm{}}}=1\frac{G(\theta _{\mathrm{}})}{r_{\mathrm{}}},$$
(15)
where $`p_{\mathrm{}}`$ is the parity of the image $`\mathrm{}`$. It is only the flux ratios that can be related to the observational data, not the magnifications themselves. For the flux ratio of image $`\mathrm{}`$ compared to image $`k`$, we obtain

$$f_k\mathrm{}=\frac{p_k\mu _k}{p_{\mathrm{}}\mu _{\mathrm{}}}=\frac{1G(\theta _{\mathrm{}})/r_{\mathrm{}}}{1G(\theta _k)/r_k}.$$
(16)
The flux ratio may be positive or negative depending on the combination of parities of the images. For flux ratios, we always compare to a reference image, which without loss of generality we take to be the first image. As we take ratios, we only need to identify images of the same parity, which in a quadruplet lens are usually on opposite sides of the lensing galaxy. Clearing the fractions, we can write the equation in a similar form to the lens equation

$$(f_1\mathrm{}1)r_1r_{\mathrm{}}=\frac{a_0}{2}\gamma _0(\theta _{\mathrm{}})+\underset{k=1}{\overset{\mathrm{}}{}}[a_k\gamma _k(\theta _{\mathrm{}})+b_k\delta _k(\theta _{\mathrm{}})],$$
(17)
with

$`\gamma _k(\theta _{\mathrm{}})`$ $`=`$ $`(1k^2)[f_1\mathrm{}r_{\mathrm{}}\mathrm{cos}(k\theta _1)r_1\mathrm{cos}(k\theta _{\mathrm{}})]`$
$`\delta _k(\theta _{\mathrm{}})`$ $`=`$ $`(1k^2)[f_1\mathrm{}r_{\mathrm{}}\mathrm{sin}(k\theta _1)r_1\mathrm{sin}(k\theta _{\mathrm{}})].`$ (18)
This provides a further $`n1`$ linear constraints on the unknowns.

Now let us introduce a compact notation to facilitate matters. We set $`\alpha _k\mathrm{}\alpha _k(\theta _{\mathrm{}})`$, where the first index denotes the Fourier component and the second index the image number. Exactly similar definitions can be made for $`\beta _k\mathrm{},\widehat{\alpha }_k\mathrm{},\widehat{\beta }_k\mathrm{},\gamma _k\mathrm{}`$ and $`\delta _k\mathrm{}`$. We are now ready to set up the matrix equation for the fitting to the positions of the images and the flux ratios. Let us write

$$C𝒙=𝒅$$
(19)
where

$$𝒅=(r_{\mathrm{}}\mathrm{cos}\theta _{\mathrm{}},\mathrm{},r_{\mathrm{}}\mathrm{sin}\theta _{\mathrm{}},\mathrm{},(f_1\mathrm{}1)r_1r_{\mathrm{}},\mathrm{})^T$$
(20)
is a vector of $`3n1`$ observables. Here, $`\mathrm{}`$ runs from $`1`$ to $`n`$ for a lens system with $`n`$ images. There is always one less equation for the flux constraints than for the image position constraints. The vector

$$𝒙=(\xi ,\eta ,a_0/2,a_1,b_1,a_2,b_2,\mathrm{}.)^T$$
(21)
contains the unknown quantities, namely the source coordinates ($`\xi ,\eta `$), and the Fourier coefficients. There are in principle infinitely many such coefficients, but in practice the Fourier series is usually terminated so that $`𝒙`$ is a vector of $`3n1`$ components as well. Finally, for the matrix $`C`$, we obtain

$$C=\left(\begin{array}{cccccccc}1& 0& \alpha _0\mathrm{}& \alpha _1\mathrm{}& \beta _1\mathrm{}& \alpha _2\mathrm{}& \beta _2\mathrm{}& \mathrm{}\\ .& .& .& .& .& .& .& \mathrm{}\\ 0& 1& \widehat{\alpha }_0\mathrm{}& \widehat{\alpha }_1\mathrm{}& \widehat{\beta }_1\mathrm{}& \widehat{\alpha }_2\mathrm{}& \widehat{\beta }_2\mathrm{}& \mathrm{}\\ .& .& .& .& .& .& .& \mathrm{}\\ 0& 0& \gamma _0\mathrm{}& \gamma _1\mathrm{}& \delta _1\mathrm{}& \gamma _2\mathrm{}& \delta _2\mathrm{}& \mathrm{}\\ .& .& .& .& .& .& .& \mathrm{}\end{array}\right)$$
(22)
Provided we can measure all $`n`$ image positions and $`n1`$ flux ratios, then $`C`$ has $`3n1`$ rows. The number of columns is in principle arbitrarily large, although in practice it usually makes sense to terminate the Fourier series so that $`C`$ is a $`3n1\times 3n1`$ matrix.

As it stands, the matrix $`C`$ is always singular $`(detC=0)`$. There exists a null space for the matrix $`C`$, i.e., there exists at least one non-vanishing vector $`𝒙_0`$ which satisfies the equation $`C𝒙_0=\mathrm{𝟎}`$. In fact, the matrix $`C`$ has (at least) a two-dimensional null space. Fortunately, the null subspace can be easily constructed. Since $`\alpha _1\mathrm{}=1`$, $`\beta _1\mathrm{}=0`$, $`\widehat{\alpha }_1\mathrm{}=0`$, $`\widehat{\beta }_1\mathrm{}=1`$, $`\gamma _1\mathrm{}=0`$ and $`\delta _1\mathrm{}=0`$ for all $`\mathrm{}`$, it is easy to verify that the null vector must be of the form $`𝒙_0=\lambda _1𝒗_1+\lambda _2𝒗_2`$ with $`𝒗_1=(1,0,0,1,0,\mathrm{})^T`$ and $`𝒗_2=(0,1,0,0,1,0,\mathrm{})^T`$.

From the physical point of view, the Fourier coefficients $`a_1`$ and $`b_1`$ simply produce a shift or transformation of the source position $`(\xi ,\eta )`$. However, both coefficients do not contribute to the surface mass density for the special case of $`\beta =1`$ since $`a_1\mathrm{cos}\theta `$ and $`b_1\mathrm{sin}\theta `$ are solutions of the homogenous differential equation $`F(\theta )+F^{\prime \prime }(\theta )=0`$ (cf. eq.(2)). We can excise this degeneracy by simply setting $`a_1=b_1=0`$, which corresponds to the choice of origin. Now we can just remove the 4th and 5th column of the matrix $`C`$ and add two higher degree Fourier coefficients to the problem to maintain $`C`$ as a $`3n1\times 3n1`$ matrix. The new matrix is usually non-singular and so a Gauss-Jordan elimination can formally solve the problem. However, it is better to use singular value decomposition (SVD) to solve for the unknown vector $`𝒙`$ (cf. Press et al. 1999). In this way, we are always guaranteed to find a numerically stable solution.

The singular value decomposition of the matrix $`C`$ is

$$C=UWV^T,$$
(23)
where the $`3n1\times 3n1`$ matrices $`U`$ and $`V`$ are orthonormal, and the matrix $`W`$ is a diagonal matrix of singular values ($`w_1,w_2,\mathrm{},w_{3n1}`$). As usual in a singular value decomposition (SVD), any very small singular values are removed (see Press et al. 1989) and this must be accompanied by a corresponding reduction in the number of unknown Fourier coefficients. Typically, any singular value whose ratio to the largest singular value is less than $`10^4`$ is removed. After deletion, let us suppose that $`m`$ singular values remain, where $`m3n1`$. If $`m=3n1`$, then the SVD solution is exact; if singular values are removed, then the solution is usually very good.

Let us remark here that one can add further higher degree Fourier coefficients and then the fitting problem becomes underdetermined. We can obtain a whole space of models all of which would equally well reproduce the observed data. It is clear – even with the severe restriction that we have made to a scale-free model with a flat rotation curve – that the variety of models that fit the data is extremely large. This cautions us against making rash statements based on goodness-of-fit to a single model.

### 3.2 Uncertainties

It is of modest interest to provide a model that fits the observational data exactly, as the uncertainties in the flux ratios may be substantial. It is of much greater interest to be able to find a set of models that can reproduce the data within the errors.

The relative image positions are usually known to extremely high accuracy. For lenses observed with the Hubble Space Telescope, the error in the relative astrometry is only $`3`$ mas. From the point of view of modelling, however, the important quantities are the relative positions of the images with respect to the centre of the mass distribution of the lensing galaxy. These coordinates are less accurately determined. The centre of the light distribution of the lensing galaxy is sometimes poorly known (especially if it is at high redshift). In any case, there is no guarantee that the center of the light corresponds to the centre of the mass. In our modelling, we always reproduce the relative image positions exactly, but we allow for an error of no more than 50 mas in the position of the center of the lens.

For the flux ratios, we distinguish between the radio and the optical data. The flux measurements in the radio may be affected by scintillation or scatter-broadening or free-free absorption. An isolated example of radio microlensing is known (Koopmans & de Bruyn 2000), but the effects of microlensing are not generally believed to be substantial in the radio. Whenever we use radio data, we assume that the typical uncertainty in the flux ratio is no more than 5 per cent. However, there are much greater problems in the optical bands. Each image might be affected by differential extinction or by microlensing. Microlensing can cause an image to be below or above its average magnification for decades (see Witt & Mao 1994 for an extensive discussion of the effects of microlensing). The recent years have seen monitoring campaigns of both Q 2237+030 and Q 0957+561, which provide irrefutable evidence for microlensing (Colley et al. 2000; Wozniak et al. 2002; Schmidt et al. 2002). There are at least a further five multiply imaged quasars for which evidence exists implicating microlensing in the optical (HE 1104-1805, PG 1115+080, H 1413+117, B 0218+357, B 1600+434). Therefore, the optical flux ratios are more uncertain and we assume that they have an error of $`<10`$ per cent.

### 3.3 Position Angle of the Lensing Galaxy Axes

The axes of the local right ascension and declination coordinates may be rotated by an arbitrary angle with respect to the major and minor axes of the lensing galaxy. This is implicitly contained in the derived Fourier coefficients.

We can expand the Fourier series as

$`{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}[a_k\mathrm{cos}(k\theta )+b_k\mathrm{sin}(k\theta )]`$ (24)
$`=`$ $`{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}[a_k^R\mathrm{cos}(k\theta +k\theta _R)+b_k^R\mathrm{sin}(k\theta +k\theta _R)],`$
where $`\theta _R`$ is the rotation angle and $`a_k^R`$ and $`b_k^R`$ are the Fourier coefficients in the rotated coordinate system. It is easy to verify that

$`a_k^R`$ $`=`$ $`a_k\mathrm{cos}k\theta _R+b_k\mathrm{sin}k\theta _R`$
$`b_k^R`$ $`=`$ $`a_k\mathrm{sin}k\theta _R+b_k\mathrm{cos}k\theta _R.`$ (25)
If the lensing galaxy has a major axis, the following condition must be satisfied:

$$G^{}(\theta _R)0G^{}(\theta _R+180^{}).$$
(26)
If the lensing galaxy looks regular, it probably has, in addition to a major axis, an approximate reflection symmetry $`G(\theta )G(\theta )`$, or equivalently $`F(\theta )F(\theta )`$ . Therefore, in the rotated system, we expect $`b_k^R0`$ for all $`k`$. In other words, the rotated system is defined by

$$\underset{k=1}{\overset{n}{}}(b_k^R)^2=\mathrm{minimum}.$$
(27)
Solving this equation yields the rotation angle of the major axis of the lensing galaxy

$$2\underset{k=1}{\overset{n}{}}ka_kb_k\mathrm{cos}(2k\theta _R)=\underset{k=1}{\overset{n}{}}k(a_k^2b_k^2)\mathrm{sin}(2k\theta _R).$$
(28)
This equation must in general be solved numerically. For galaxy-like mass distributions, the dominant higher-order terms in the Fourier expansion are those with $`k=2`$ and then the solution is analytic, namely

$$\mathrm{tan}2\theta _R=\frac{b_2}{a_2}.$$
(29)
For applications to individual lenses, we note that the position angle as conventionally defined by astronomers (North through East) is $`90^{}\theta _R`$.

### 3.4 Mass inside the Einstein Ring

Evans & Witt (2001) showed that the mass inside the Einstein ring must be very close to the mass inside the critical curve for the special case of $`\beta =1`$, and can be written as

$$M_EM_{\mathrm{crit}.\mathrm{curve}}=\frac{1}{2\pi }_0^{2\pi }G^2(\theta )𝑑\theta .$$
(30)
We therefore deduce that

$$M_E\frac{a_0^2}{4}+\underset{i=2}{\overset{N}{}}(k^21)^2(a_k^2+b_k^2).$$
(31)
To apply this to a real lens, we have to scale the result by the factor $`\mathrm{\Sigma }_{\mathrm{crit}}D_\mathrm{d}^2`$ containing the normalized surface mass density which is is given by

$$\mathrm{\Sigma }_{\mathrm{crit}}=\frac{c^2D_\mathrm{s}}{4\pi GD_\mathrm{d}D_{\mathrm{ds}}}.$$
(32)
Here, $`D_\mathrm{d},D_\mathrm{s}`$ and $`D_{\mathrm{ds}}`$ are the distance to the deflector, the distance to the source and the distance from the deflector to the source respectively. Given the redshift of the lens and the source, this (32) enables the projected mass within the Einstein ring to be estimated.

## 4 AN APPLICATION: ANOMALOUS FLUX RATIOS

In this section, we consider three of the lens systems for which substructure has been claimed either by Metcalf & Zhao (2002) or by Dalal & Kochanek (2002) or by Chiba (2002).

### 4.1 Q 2237+030 (The Einstein Cross)

The data on the four images of Q 2237+030 (the Einstein Cross) are listed in Table 1. Much of this comes from the CASTLES survey <sup>1</sup><sup>1</sup>1http://cfa-www.harvard.edu/castles/. Q 2237+030 is an unusual lens because so much more information is available for the lensing galaxy than is customary. In particular, the redshift of the lens is so small ($`z=0.039`$) that the galaxy’s light distribution can be measured (e.g., Wyithe et al. 2002). The lens is a face-on barred galaxy (e.g., Schmidt, Webster & Lewis 1998). The optical fluxes are known to be affected by microlensing (Irwin et al. 1989; Wambsganss, Schneider & Paczyński 1990; Wozniak et al. 2000), which causes each image to differ from the average magnification. Typically, we expect that the images are very likely below the average magnification if they are in a quiescent state (for years), while the images are very likely above the average magnification if they are in a more active state (cf. Witt & Mao 1994).

It is often argued that the radio and mid-infrared fluxes are more reliable because the source emitting region is more extended and therefore less affected by microlensing. If this is the case, then gravitational lens models should use the radio and mid-infrared flux ratios in preference to optical. The flux ratios in the radio were measured by Falco et al. (1996) with a signal-to-noise of $`24`$. The mid-infrared flux ratios have recently become available thanks to Agol, Jones & Blaes (2000). The flux ratios in radio and mid-infrared are in good agreement, but differ from the optical flux ratios.

Figure 2 shows the critical curves for four models of Q 2237+030. From eq. (15), the critical curve is also an equidensity contour when the rotation curve is flat. As the models are scale-free, the equidensity contours retain the same shape independent of position, although the size of course varies. In each case, the dashed curve describes a model that exactly reproduces the image positions and flux ratios. It is found by solution of the matrix equation (19) using singular value decomposition (SVD). Figure 3 shows the corresponding caustics of the exact solutions of the critical curves (dashed lines). The more irregular critical curves (I and H band data) correspond to caustics with over-folds like swallowtails and butterflies. However, in each case, the source position still lies inside the four image region. Unwanted extra images are not produced.

We assume that the measured flux ratios in the optical band are in error by $`<10`$ per cent due to microlensing. The relative positions of the image are always maintained within the formal errors (given in Table 1), but the center of the lensing galaxy may be displaced by up to 50 mas. We find the model for which

$$b_2^2+\underset{i3}{}[a_i^2+b_i^2]$$
(33)
is minimised. This suppresses Fourier components higher than the bar-like $`m=2`$ component. In other words, out of all the solutions that reproduce the data, we choose the one that most looks like a galaxy. The models so produced are shown in bold lines in Figure 2; the details of the solutions are listed in Table 2. Finally, for two of the bold line solutions, we show how the image positions and the flux ratios are reproduced in Table 3. All the observable quantities are extremely well reproduced by the model, including the flux ratios.

Two points emerge clearly from this. First, the equidensity contours implied by reproducing the exact radio or mid-infrared fluxes are already in good agreement. The position angle of the lensing galaxy as inferred from eq. (26) using the radio or mid-infrared solutions is $`67.3^{}`$ (measured from North through East). For comparison, Trott & Webster (2002) report that the position angle of the bar in the lensing galaxy is $`39^{}`$, while the major axis of the galaxy in its outer parts is $`77^{}`$. Second, the equidensity contours implied by exactly reproducing the I and the H band optical flux ratios are distorted. However, there are more regular solutions once an allowance is made for the effects of microlensing on the flux ratios.

Distortions of isophotes from elliptical form are generally given in terms of Bender et al.’s (1989) coefficients. These are usually calculated for early-type galaxies and have absolute values less than 0.05 (see e.g., Binney & Merrifield 1998, section 4.3). In fact, barred galaxies, such as the lens in Q 2237+030, may have substantially larger deviations, which moreover depend on position along the major axis. The method of calculating these coefficients for densities given as Fourier expansions is elaborated in Appendix C; we quote the results here. For the radio contours, Bender’s coefficients are $`a_3^B/a_0=0.008`$ and $`a_4^B/a_0=0.005`$; for the mid-infrared, they are $`a_3^B/a_0=0.001`$ and $`a_4^B/a_0=0.011`$. These values are well within the observed range, emphasising that our solutions for the shapes of the equidensity contours are realistic. By comparison, for the distorted I band contours, the coefficients are $`a_3^B/a_0=0.012`$ and $`a_4^B/a_0=0.034`$.

All this provides strong evidence for the point of view advanced in Agol et al. (2000) that the radio and mid-infrared fluxes are trustworthy and that the optical flux ratios are discrepant because of microlensing. All the data are internally consistent and once reasonable error bars are placed on the flux ratios, there is no need to invoke substructure to explain the data (c.f. Metcalf & Zhao 2002).

### 4.2 PG 1115+080

Let us now examine the lens system PG 1115+080, which has been claimed to provide evidence for substructure by Dalal & Kochanek (2002), Metcalf & Zhao (2002) and Chiba (2002). The observational data on PG 1115+080 are listed in Table 4. The I, V and H band fluxes are provided from the CASTLES survey. There is also earlier data on the I band fluxes thanks to Kristian et al. (1993, hereafter K93).

There is already evidence that the flux ratios may come with substantial uncertainties. For example, Foy, Bonneau & Blazit (1985) found component A2 of PG 1115+080 was brighter than A1 in 1984 March-May, whereas most subsequent investigators have found A2 to be fainter than A1. The conclusion advanced by Christian, Crabtree & Waddell (1987) is that the images in PG 1115+080 are being affected either by microlensing or by time-delayed intrinsic variations in the quasar or both. Witt, Mao & Schechter (1995) showed that the likelihood of observing microlensing in images A1 and A2 is much higher than in C and D. Notice too that the CASTLES and the Kristian et al (1993) I band fluxes in Table 4 are not consistent within the stated error bars, the discrepancies almost certainly being due to microlensing. Radio emission would be less affected by microlensing, but PG 1115+080 is radio-quiet and so the optical fluxes have to be used.

Impey et al. (1998) have carried out deep Hubble Space Telescope imaging of PG 1115+080 in the infrared. They found that the primary lens is a nearly round elliptical galaxy with an ellipticity of less than 0.07. The lens is part of a group of galaxies, whose effects must be included in successful models. The group is modelled by a uniform external shear with a magnitude of $`0.103`$ and a direction of $`65.8^{}`$ measured from west to north, following Schechter et al. (1997). The extension of our SVD method to the case of uniform shear is sketched in Appendix D. We implement our fitting procedure, deleting two singular values for numerical stability and finding the solution that minimises (33), but still reproduces the data within our assumed errors. The panels of Figure 4 show the equidensity contous of the resulting models. Again, the relative astrometry of the images is recovered to within the formal error ($`4`$ mas), the lens centre may be misligned by up to 50 mas, while the optical flux ratios are reproduced to within 10 per cent. The comparison between model predictions and observed quantities for the H and I band (K93) solutions are given in Table 6. The contours implied by reproducing the H or I band (K93) fluxes are in good agreement and look plausibly like the projected mass distribution of a very round early-type galaxy, in accordance with Impey et al.’s (1998) imaging. There is some slight distortion in the equidensity contours that reproduce the V or I band (CASTLES) flux ratios. The caustic networks are shown in Figure 5 to demonstrate that unwanted extra images are not formed.

Again, the shapes of the equidensity contours in Figure 4 do exhibit small deviations from pure elliptical form, which can be quantified by Bender et al.’s (1989) coefficients. Using the method of Appendix C, we find $`a_3^B/a_0=0.009`$ and $`a_4^B/a_0=0.004`$ for the H band and $`a_3^B/a_0=0.015`$ and $`a_4^B/a_0=0.004`$ for the I band (K93) contours. There is some data on the shapes of dark haloes of spiral galaxies from numerical simulations (e.g., Heyl, Hernquist & Spergel 1994), which supports the notion that $`|a_3|<0.05`$ and $`|a_4|<0.05`$. So, the amplitudes of the higher order terms in our solutions for the lens galaxy do lie easily within the range expected for dark haloes of spiral galaxies.

It is clear from the smooth models that we have constructed that there is no really compelling evidence for substructure in this lens at the moment, despite abundant claims to the contrary (Dalal & Kochanek 2002; Metcalf & Zhao 2002; Chiba 2002; Kochanek & Dalal 2003). There are smooth and realistic models that do fit the data within the likely uncertainties.

### 4.3 B 1422+231

Patnaik et al. (1992) discovered the gravitational lens B 1422+231 in a radio survey. Table 7 lists some of the available data on the images, largely provided by the CASTLES group. It was realised by both Hogg & Blandford (1994) and Kormann, Schneider & Bartelmann (1994) that it is difficult to reproduce the flux ratios of this lens system with a simple model. Hogg & Blandford used a mixture of isothermal spheres and point masses whereas Kormann et al. used a model with external shear. The latter strategy has been followed by most modern modellers.

Using our SVD method, it was impossible to reproduce the reported errors in the relative astrometry, together with the assumed uncertainties in the flux ratios (5 per cent for the radio, 10 per cent for the optical) at least with physical models. Unphysical models are readily identified as the equidensity contours are self-intersecting. Accordingly, for this lens only, the errors in the relative astrometry were inflated by a factor of 10 to $`30`$ mas. Figure 6 now shows the results of our fitting procedure in four flux bands, but none of the results can be said to be at all satisfactory. The equidensity contours are strongly distorted from those expected in galaxy-like models. Just as alarming as the shape of the contours is the fact that the fitted parameters do not have consistent values from one band to the next.

The conclusion to be drawn by the failure of the fitting procedure is that – even taking into account observational uncertainties – this lens system requires additional ingredients. Possibly, as Mao & Schneider (1998) were the first to suggest, substructure may be causing a dramatic effect on the flux ratios. Possibly, as Kundic et al. (1997) suggest, the effects of external shear caused by nearby galaxies and clusters of galaxies may also need to be included in the analysis.

## 5 DISCUSSION

The main aim of this paper is to present a new, simple approach to fitting a model to the data on a gravitational lens. Given the image positions and the flux ratios (if available), a model that exactly fits that data can be constructed by a simple matrix inversion (preferably by singular value decomposition, available in standard numerical libraries). All the fitting parameters enter linearly into the equation. No complicated non-linear $`\chi ^2`$-fitting needs to be done.

The lensing galaxy is always assumed to be scale-free with a flat rotation curve. There is already a substantial body of evidence from fitting of lenses that early-type galaxies have nearly isothermal profiles (Kochanek 1995; Muñoz et al. 2001; Cohn et al. 2001; Rusin et al. 2002). In particular, the absence of central images strongly constrains the lensing potential of early-type galaxies to be isothermal or steeper and the size of any core region to be small (e.g., Rusin & Ma 2001; Evans & Hunter 2002). There is both stellar dynamical (e.g., Gerhard et al. 2001) and X-ray (Fabbiano 1997) evidence that early-type galaxies have flat rotation curves out to at least 4 effective radii. Hence, our assumption of a flat rotation curve seems exactly what the data require.

Unusually, our algorithm allows full flexibility as regards the angular structure of the lensing potential. Earlier fitting procedures have allowed the lensing potential to have different radial structure (for example, different power-law profiles), but usually only simple forms of angular structure (for example, constant external shear). In our models, the radial structure is always fixed as isothermal, but the shape of the density contours is completely arbitrary. The lens model has a flexible number of free parameters, and of course includes the isothermal elliptic density and potential models which have a distinguished history in gravitational lens modelling. The fit allows direct deduction of the mass contour of the lensing galaxy. Although we have presented examples based on quadruplet lens systems, the algorithm can be adapted for a two image or a higher image model. A higher image model might be needed, for example, in Q 0957+561, where 10 sub-image positions can be detected in the radio band. Another good application might be MG 0414+0534, where some substructure of the images are detected (cf. Trotter, Winn & Hewitt 2000 and their fitting approach). Both lenses seem heavily distorted by external shear. A shear term can be added to our algorithm. If the shear is known, then the lens and flux ratio equations remain linear.

The flux ratios may enter directly into the fit depending on whether they are available or trustworthy. We stress that flux ratios often do not provide good model constraints. The reason for this is easy to understand. The flux ratios depend on the second derivatives of the lensing potential, whereas time delays and image positions are proportional to the potential and its first derivative respectively. Therefore, flux ratios are particularly sensitive to graininess in the gravitational potential. This often manifests itself as microlensing. The flux ratios are also affected by the uncertain differential extinction corrections that must be applied to each image. Radio and mid-infrared fluxes are generally more reliable than optical. However, effects such as scintillation, scatter-broadening and free-free absorption may affect the radio fluxes in some of the lenses. If flux ratios are incorporated into a fit, it is vital that there is a realistic treatment of the errors.

Our fitting algorithm permits the construction of a whole class of degenerate models all of which can fit the image position and flux ratios exactly. In addition, the models are also degenerate in the time delay since the delay depends only on the distance of the image positions to the centre of the lensing galaxy (cf. Witt, Mao & Keeton 2000). Such degeneracy is easy to understand pictorially. Given a Fermat time delay surface, we need only keep the values, the derivatives and the second derivatives of the surface at the image positions fixed. We then have enormous freedom to move the surface (subject only to the constraints that no additional images are introduced and that no negative mass density is produced).

We have used our new fitting algorithm to examine critically some of the claims made recently for anomalous flux ratios. The procedure used by both Dalal & Kochanek (2002) and Metcalf & Zhao (2002) was to show that a family of simple models did not reproduce the flux ratios. However, it does not follow from this that substructure is necessary; it may just be that the simple model did not have enough flexibility to provide a good match. For example, of the five lenses studied by Metcalf & Zhao (2002), they themselves concluded that one (MG 0414+0534) could be satisfactorily explained by a smooth model. In this paper, we have demonstrated that the data on two others (PG 1115+080, Q 2237+030) are consistent with smooth galaxy-like models, especially when a realistic treatment of the errors in the flux ratios is incorporated.

We note that the fraction of mass in substructure expected in galaxy haloes is very small ($`<5`$ per cent) and it occurs overwhelmingly in the outlying portions (e.g., Moore et al. 1999; Ghigna et al. 2000). Substructure evolves as it is subjected to tides, impulsive collisions and dynamical friction. Tidal disruption becomes important as soon as the mean density of the galaxy is equal to the density of the substructure. Similarly, the dynamical friction timescale scales with the square of the galactocentric radius (see e.g., Binney & Tremaine 1987). So, both tides and dynamical friction are efficient at erasing substructure in the inner parts. This effect is referred to as “anti-biasing” by Ghigna et al. (2000). So, on dynamical grounds, little substructure is projected within the Einstein radius. Taking the example of PG 1115+080, the projected Einstein radius is just $`3.6`$ kpc. The fraction of mass in substructure projected within a cylinder of radius $`3.6`$ kpc is clearly much lower than the global fraction of $`5`$ per cent, which pertains to the entire dark halo of total extent $`200`$ kpc. It is crucial that lensing calculations do not assume an everywhere uniform fraction of substructure in the halo, as this does not take into account the “anti-biased” spatial distribution of substructure and therefore necessarily over-emphasises the importance of the effects of substructure on flux ratios. Notice that an interesting consequence of the spatial distribution is that anomalous flux ratios are more likely to occur for lenses with larger angular separation, as these have larger Einstein radii.

The real question at issue is the following. Suppose a simple lens model (such as a simple isothermal ellipsoid plus shear) does not adequately fit the data on positions and flux ratios. What can be legitimately deduced? The difficulty is that there are many modifications of the simple model that remove the discrepancy. One of these is substructure, as pointed out by Dalal & Kochanek and Metcalf & Zhao. As shown in this paper, another is higher order multipoles in the lensing mass, such as diskiness, boxiness, lopsidedness and barredness (see also Möller, Hewett & Blain 2003, who make a similar point). It is crucial to develop techniques to distinguish between higher order quadrupoles on the one hand and substructure on the other. From our modelling, we tend to agree that the substructure candidate B 1422+231 originally pointed out by Mao & Schneider (1998) is strong. Using a novel application of the cusp relation, Keeton, Gaudi & Petters (2002) have argued that B 2045+265 and RX J 0911+0551 may be two further good candidates.

It is an outstanding problem to predict – for different cosmologies – how many quadruplets may have anomalous flux ratios. Accurate calculations will become possible only when distributions of mass and position of substructure become available from high resolution simulations. It will be interesting to see whether the numbers of lenses for which substructure is currently being claimed are compatible with cold dark matter theories.

## Acknowledgments

We acknowledge with gratitude the remarkable service provided to the community by the CASTLES website (C.S. Kochanek, E.E. Falco, C. Impey, J. Lehar, B. McLeod, H.-W. Rix). NWE thanks the Royal Society for financial support and HJW thanks the sub-Department of Theoretical Physics, Oxford for their hospitality during working visits. NWE acknowledges Ian Browne, Leon Koopmans and Jerry Ostriker for a number of interesting discussions. Shude Mao and Olaf Wucknitz are thanked for a helpful comments on the draft manuscript.

## Appendix A Full Solution by the Variation of Parameters

Observationally speaking, it is the angular structure of the density contours $`G(\theta )`$ that is more accessible. It is sometimes useful to be able to generate the corresponding angular structure of the lensing potenial $`F(\theta )`$. This can be done by solving eq (2) using the method of variation of the parameters (e.g., Bronshtein & Semendyayev 1998, section 3.3.1.3.4; Evans & Witt 2001) to give

$`F(\theta )`$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}(\beta \theta )}{\beta }}\left[C_1+{\displaystyle _0^\theta }G(\vartheta )\mathrm{cos}(\beta \vartheta )d\vartheta \right]`$ (34)
$``$ $`{\displaystyle \frac{\mathrm{cos}(\beta \theta )}{\beta }}\left[C_2+{\displaystyle _0^\theta }G(\vartheta )\mathrm{sin}(\beta \vartheta )d\vartheta \right].`$
This establishes the relation between the potential and the surface mass density in terms of two constants $`C_1`$ and $`C_2`$, whose values will be given shortly.

Therefore, the deflection angle has components

$`\varphi _x`$ $`=`$ $`r^{\beta 1}\mathrm{sin}(\beta 1)\theta \left[C_1+{\displaystyle _0^\theta }G(\vartheta )\mathrm{cos}(\beta \vartheta )d\vartheta \right]`$
$`r^{\beta 1}\mathrm{cos}(\beta 1)\theta \left[C_2+{\displaystyle _0^\theta }G(\vartheta )\mathrm{sin}(\beta \vartheta )d\vartheta \right],`$
$`\varphi _y`$ $`=`$ $`r^{\beta 1}\mathrm{cos}(\beta 1)\theta \left[C_1+{\displaystyle _0^\theta }G(\vartheta )\mathrm{cos}(\beta \vartheta )d\vartheta \right]`$ (35)
$`+r^{\beta 1}\mathrm{sin}(\beta 1)\theta \left[C_2+{\displaystyle _0^\theta }G(\vartheta )\mathrm{sin}(\beta \vartheta )d\vartheta \right].`$
Another expression can be given for the deflection angle, namely

$`\varphi _x`$ $`=`$ $`{\displaystyle \frac{\beta x}{r^{2\beta }}}F(\theta ){\displaystyle \frac{y}{r^{2\beta }}}F^{}(\theta ),`$
$`\varphi _y`$ $`=`$ $`{\displaystyle \frac{\beta y}{r^{2\beta }}}F(\theta )+{\displaystyle \frac{x}{r^{2\beta }}}F^{}(\theta ).`$ (36)
By substituting the Fourier series into eqs. (35)–(36), we obtain expressions for the constants

$`C_1`$ $`=`$ $`{\displaystyle \underset{k}{}}kb_k,`$
$`C_2`$ $`=`$ $`\beta {\displaystyle \frac{a_0}{2}}{\displaystyle \underset{k}{}}\beta a_k.`$ (37)
This completes the solution written down in eq. (34). Given the angular structure of the density $`G(\theta )`$, we can find that of the potential $`F(\theta )`$, and vice versa.

## Appendix B Fourier Coefficients of Elliptical Distributions

Let us consider the Fourier expansion of the elliptical distribution in the flat rotation curve case ($`\beta =1`$)

$$F(\theta )=A(\mathrm{cos}^2\theta +q^2\mathrm{sin}^2\theta )^{1/2}.$$
(38)
The Fourier coefficients are

$`a_k`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _0^{2\pi }}F(\theta )\mathrm{cos}(k\theta )𝑑\theta `$ (39)
$`b_k`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _0^{2\pi }}F(\theta )\mathrm{sin}(k\theta )𝑑\theta .`$ (40)
However, the only non-vanishing Fourier coefficients are $`a_{2n}`$ where $`n`$ is an integer. The lowest-order terms are

$`a_0`$ $`=`$ $`{\displaystyle \frac{4A}{\pi }}K`$
$`a_2`$ $`=`$ $`{\displaystyle \frac{4A}{\pi (1q^2)}}[(1+q^2)K2q^2E]`$
$`a_4`$ $`=`$ $`{\displaystyle \frac{4A}{3\pi (1q^2)^2}}[(3q^2+1)(q^2+3)K8q^2(1+q^2)E]`$
$`a_6`$ $`=`$ $`{\displaystyle \frac{4A}{15\pi (1q^2)^3}}[(1+q^2)(15+98q^2+15q^4)K`$
$`2q^2(23+82q^2+23q^4)E]`$
$`a_8`$ $`=`$ $`{\displaystyle \frac{4A}{105\pi (1q^2)^4}}[32q^2(1+q^2)(11+74q^2+11q^4)E`$
$`+(105+1436q^2+3062q^4+1436q^6+105q^8)K],`$
where $`K`$ and $`E`$ denote the following integrals

$`K`$ $`=`$ $`{\displaystyle _0^{\pi /2}}{\displaystyle \frac{d\vartheta }{(1(1q^2)\mathrm{sin}^2\vartheta )^{1/2}}}`$
$`E`$ $`=`$ $`{\displaystyle _0^{\pi /2}}𝑑\vartheta (1(1q^2)\mathrm{sin}^2\vartheta )^{1/2}.`$ (42)
For elliptical potential models, this gives the Fourier series for the angular function $`F(\theta )`$. Note that the expressions become more and more cumbersome with increasing order. However, the Fourier series converges (numerically) extremely rapidly.

For elliptical density models, exactly the same coefficients hold good but for the Fourier expansion for $`G(\theta )`$, the shape function in the density. Of course, the Fourier coefficients of $`G(\theta )`$ are directly related to the Fourier coefficients of the potential contour $`F(\theta )`$ via eqs. (4) and (5).

## Appendix C Calculation of the Isophotal Distortions

Most giant elliptical galaxies have isophotes that are very well approximated by ellipses. The same is of course not true for disk or barred or interacting galaxies. Data from numerical simulations of halo formation are sparse, but suggest that dark haloes may also be well-approximated by ellipses (e.g., Heyl et al. 1994). So, it is useful to be able to estimate the size of the distortions of our fits from pure ellipses. This is done by computing Bender et al.’s (1989) coefficients, defined as

$`a_k^B`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _0^{2\pi }}\left(r_\mathrm{i}(\theta )r_\mathrm{e}(\theta )\right)\mathrm{cos}(k\theta ),`$ (43)
$`b_k^B`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _0^{2\pi }}\left(r_\mathrm{i}(\theta )r_\mathrm{e}(\theta )\right)\mathrm{sin}(k\theta ),`$ (44)
where $`r_\mathrm{i}`$ is the polar equation of the isophote and $`r_\mathrm{e}`$ is the best fitting ellipse (see also Binney & Merrifield 1998). In particular, a positive (or negative) $`a_4^B`$ means that the isophotes are disk-like (or boxy).

From eqn (5), we deduce the equation of the isophote as

$$r_\mathrm{i}=\frac{a_0}{2}+\underset{k=2}{\overset{\mathrm{}}{}}[a_k(1k^2)\mathrm{cos}(k\theta )+b_k(1k^2)\mathrm{sin}(k\theta )].$$
(45)
Here, we have specialised to the case of a flat rotation curve by setting $`\beta =1`$. As shown in Section 3.3, we can transform to the coordinate system aligned with the major and minor axis and so we can set $`b_2=0`$.

Now we must determine the best fitting ellipse, which is the one which reproduces the same Fourier coefficients $`a_0`$ and $`a_2`$ as our isophote. In other words, we solve for the amplitude $`A`$ and axis ratio $`q`$ in eqn (38),

$`a_0`$ $`=`$ $`{\displaystyle \frac{4A}{\pi }}K,`$
$`3a_2`$ $`=`$ $`{\displaystyle \frac{4A}{\pi (1q^2)}}[(1+q^2)K2q^2E].`$ (46)
These equations do not have analytic solutions, but the roots can be easily obtained numerically. Having solved for $`A`$ and $`q`$, our isophote $`r_\mathrm{i}`$ and the best fitting ellipse $`r_\mathrm{e}`$ have Fourier series which differ only for coefficients $`a_k`$ and $`b_k`$ with $`k3`$.

It is now straightforward to derive formulae for Bender et al.’s coefficients. Tabulated values are available in the literature only for the $`a_3^B`$ and $`a_4^B`$, so we confine explicit formulae to these cases:

$`a_3^B`$ $`=`$ $`8a_3,`$ (47)
$`a_4^B`$ $`=`$ $`15a_4{\displaystyle \frac{4A}{3\pi (1q^2)^2}}[(3q^2+1)(q^2+3)K8q^2(1+q^2)E]`$
Numerical values for the lens galaxies in Q 2237+030 PG 1115+080 are given in the main text. In both cases, the amplitudes of $`a_3^B`$ and $`a_4^B`$ are within the expected ranges, confirming that our solutions are realistic.

## Appendix D Inclusion of External Shear

In this appendix, we give the equations incorporating the effects of external shear. The lens equation becomes

$$\xi =x+\gamma _1x+\gamma _2y\frac{\varphi }{x},\eta =y+\gamma _2x\gamma _1y\frac{\varphi }{y}.$$
(48)
The vector of observables becomes

$`𝒅`$ $`=`$ $`((1+\gamma _1)r_{\mathrm{}}\mathrm{cos}\theta _{\mathrm{}}+\gamma _2r_{\mathrm{}}\mathrm{sin}\theta _{\mathrm{}},\mathrm{},`$ (49)
$`(1\gamma _1)r_{\mathrm{}}\mathrm{sin}\theta _{\mathrm{}}+\gamma _2r_{\mathrm{}}\mathrm{cos}\theta _{\mathrm{}},\mathrm{},`$
$`(f_1\mathrm{}1)r_1r_{\mathrm{}}(1\gamma _1^2\gamma _2^2),\mathrm{})^T,`$
where $`\mathrm{}`$ runs from $`1`$ to $`n`$ for a lens system with $`n`$ images. The matrix $`C`$ retains the same form as given in eq (22), but the coefficients $`\gamma _k\mathrm{}`$ and $`\delta _k\mathrm{}`$ become

$`\gamma _k(\theta _{\mathrm{}})`$ $`=`$ $`(1k^2)[f_1\mathrm{}r_{\mathrm{}}\mathrm{cos}(k\theta _1)W(\theta _1)r_1\mathrm{cos}(k\theta _{\mathrm{}})W(\theta _{\mathrm{}})]`$
$`\delta _k(\theta _{\mathrm{}})`$ $`=`$ $`(1k^2)[f_1\mathrm{}r_{\mathrm{}}\mathrm{sin}(k\theta _1)W(\theta _1)r_1\mathrm{sin}(k\theta _{\mathrm{}})W(\theta _{\mathrm{}})],`$
with

$$W(\theta )=1+\gamma _1\mathrm{cos}2\theta _{\mathrm{}}+\gamma _2\mathrm{sin}2\theta _{\mathrm{}}.$$
(50)
The vector of unknowns is unchanged from eq (21). For a given external shear, the problem remains linear and can be solved by SVD.


format: text

text:
# An asymptotic expansion for the discrete harmonic potential

## 1. Introduction

The discrete harmonic potential $`a`$ of a random walk $`R`$, starting from $`0`$, on a lattice $`Z`$ can be most easily defined using

$$a(z):=\underset{n=0}{\overset{\mathrm{}}{}}P_n(0)P_n(z)$$
(1)
where $`P_n(z)`$ is the probability of $`R`$ to be at $`z`$ at the $`n`$th step. It is easy to conclude from (1) that $`\mathrm{\Delta }a=\delta _{\{0\}}`$ where $`\mathrm{\Delta }`$ is the discrete Laplacian for the walk $`R`$. Less obvious are the facts that $`a`$ is positive, and that these two properties determine $`a`$ up to the addition of a positive constant. “Positive” may be replaced by “with sub-linear growth” in dimension greater than $`1`$.

The discrete harmonic potential is an interesting quantity, strongly related to the discrete Green function $`g(x,y)`$: on the one hand, $`a(x)=g(x,x)`$ where $`g`$ is the Green function corresponding to the set $`\{0\}`$ \[S76, P11.6, page 118\]. On the other hand, $`g`$ for arbitrary sets can be calculated if $`a`$ is known — when the set is finite there is an explicit formula, \[S76, T14.2, page 143\].

Thus it becomes interesting to calculate or estimate $`a`$. Except in dimension $`1`$ where $`a(x)=|x|`$, the problem is far from trivial. Take as an example the harmonic potential of the regular random walk on $`^2`$ which has a logarithmic nature. The estimate of $`a`$ most common in the literature is

$$a(z)=\frac{2}{\pi }\mathrm{log}|z|+\lambda +O(|z|^2)$$
(2)
where $`\lambda `$ is some number, which can be expressed using Euler’s $`\gamma `$ constant, $`\lambda =\frac{2}{\pi }\gamma +\frac{1}{\pi }\mathrm{log}8`$. The earliest proof of (2) we are aware of is in Stöhr \[S49\]. A more accessible proof of a weaker result, giving an error estimate of $`o(1)`$, can be found in Spitzer \[S76, P12.3, page 124\]. Spitzer’s proof (which is not specific to the lattice $`^2`$), like Stöhr’s, relies on the fact that the Fourier transform $`\widehat{P_n}`$ has an explicit formula, which can be summed to give a “pseudo-Fourier” representation of $`a`$ as the explicit integral

$$a(z)=\frac{1e^{i\xi ,z}}{1\widehat{P_1}(\xi )}𝑑\xi ;$$
(3)
indeed, this is how he proves that the sum in (1) converges in the first place. Fukai and Uchiyama used this technique to find the next order approximation — in the case of the simple random walk on $`^2`$ it is $`\frac{\mathrm{Re}z^4}{6\pi |z|^6}`$ — and to show that an asymptotic polynomial expansion exists. See \[FU96\] for the two dimensional case and \[U98\] for the case $`d3`$.

The purpose of this note is to give two alternative approaches to the computation of $`a`$, which allow to get asymptotics of arbitrary precision. For example, here are the first few terms for the regular random walk on $`^2`$:

$`a(z)`$ $`=`$ $`{\displaystyle \frac{2}{\pi }}\mathrm{log}|z|+\lambda \mathrm{Re}({\displaystyle \frac{z^4}{6\pi |z|^6}}+`$ (4)
$`+{\displaystyle \frac{3z^4}{20\pi |z|^8}}+{\displaystyle \frac{5z^8}{24\pi |z|^{12}}}+`$
$`+{\displaystyle \frac{51z^8}{56\pi |z|^{14}}}+{\displaystyle \frac{35z^{12}}{36\pi |z|^{18}}}+`$
$`+{\displaystyle \frac{217z^8}{160\pi |z|^{16}}}+{\displaystyle \frac{45z^{12}}{4\pi |z|^{20}}}+{\displaystyle \frac{1925z^{16}}{192\pi |z|^{24}}}+`$
$`+O(|z|^{10})).`$
The first approach is a brute-force one which starts from (1): $`P_n`$ may be represented as a sum of multinom coefficients, and plugging in Stirling’s expansion, and some elementary algebra gives the representation (4). Unfortunately, the algebra involved is too long. In effect, for all but the first coefficient (the $`\frac{\mathrm{Re}z^4}{6\pi |z|^6}`$), it is simply not practical to do the calculations by hand. We demonstrate this approach in two ways: we use it to prove that a polynomial approximation exists, reproducing the results of \[FU96, U98\] in the case the random walk is bounded (Fukai and Uchiyama show these results under the weaker assumption that each step of the random walk has $`K+d+ϵ`$ moments, where $`K`$ is the required precision in the expansion, $`d`$ is the dimension and $`ϵ>0`$). See theorem 1 on page 1; and we use it to get some explicit constants and high order expansions for other walks. See section 4.

The second approach is to approximate the discrete Laplacian $`\mathrm{\Delta }`$ using an appropriate differential operator. For example, for the regular random walk on $`^2`$ we have, from the Taylor expansion,

$`(\mathrm{\Delta }f)(z)`$ $`=`$ $`\frac{1}{4}\left(f(z+1)+f(z1)+f(z+i)+f(zi)\right)f(z)`$ (5)
$`=`$ $`{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{^2f}{x^2}}+{\displaystyle \frac{^2f}{y^2}}\right)+{\displaystyle \frac{1}{48}}\left({\displaystyle \frac{^4f}{x^4}}+{\displaystyle \frac{^4f}{y^4}}\right)+\mathrm{}`$
In fact, $`\mathrm{\Delta }=\frac{1}{2}(\mathrm{cosh}\frac{}{x}+\mathrm{cosh}\frac{}{y})𝐈`$, but we will not use this representation. Given that the expansion of $`a`$ has an appropriate polynomial form, the expansion (5) allows to get many equations on the coefficients in (4). These equations are of course insufficient — the equation $`\mathrm{\Delta }a=\delta _{\{0\}}`$ does not determine $`a`$ uniquely, and the condition that $`a`$ is positive or slowly increasing is hard to encode into the differential equations. However these differential equations allow to prove many of the apparent properties of (4), most notably that the coefficient of $`z^l|z|^k`$ for $`l<\frac{1}{2}k`$ is always zero. This is theorem 2. Unfortunately, it is still not clear why all the coefficients are negative.

Returning to the case of the regular random walk on $`^2`$, it can be seen that the equation $`\mathrm{\Delta }a=\delta _{\{0\}}`$ combined with the $`\frac{\pi }{2}`$-rotational symmetry gives that the values of $`a`$ are uniquely defined by their values on one diagonal, say $`\{m+im\}_{m=0}^{\mathrm{}}`$, and vice versa, any arbitrary sequence on the diagonal can be extended to a symmetric function $`f`$ with $`\mathrm{\Delta }f=\delta _{\{0\}}`$ using a very simple recursion. On the other hand, it turns out that the values of $`a`$ on the diagonal can be calculated explicitly using the integral (3) to get

$$a(m+im)=\frac{4}{\pi }\left(1+\frac{1}{3}+\frac{1}{5}+\mathrm{}+\frac{1}{2m1}\right)$$
(6)
\[S76, chapter 15\]. Indeed we have just sketched an algorithm for calculating any value of $`a(z)`$: this is the McCrea-Whipple algorithm \[MW40\]111Can be found in Spitzer as well: \[S76, chapter 15, page 148\]. As is now evident, all the values are of the form $`n+\frac{1}{\pi }q`$, $`n`$ and $`q`$. Interestingly, when not on the diagonal, $`n`$ and $`q`$ increase exponentially, even though the result is only logarithmic in size222No, this is not a very efficient way to calculate $`\pi `$, it requires at least $`k^3`$ operations to get a precision of $`k`$ digits…. As for the coefficients of (4), it turns out — perhaps unsurprisingly — that (6) can be used to complete the few coefficients that cannot be deduced from (5). This is theorem 3.

There is a third approach to the problem which we shall not discuss at length as it seems inferior to both former ones. Roughly it goes as follows (for the case of the regular random walk on $`^2`$). “Guess” that the solution is approximately $`\frac{2}{\pi }\mathrm{log}|z|`$. We discretize to $`^2`$ in the natural way: define

$$f(z):=\{\begin{array}{cc}\frac{2}{\pi }\mathrm{log}|z|\hfill & z0\hfill \\ 1\hfill & z=0\hfill \end{array}.$$

Actually, it is not necessary to guess the exact value $`\frac{2}{\pi }`$, other constants not too different also work. A calculation can show that

$$\underset{z^2}{}|(\mathrm{\Delta }f\delta _{\{0\}})(z)|<1.$$
(7)
This allows us to write a series of corrections of $`f`$ as follows: $`f_1=f`$ and then

$$f_n=f_{n1}(\mathrm{\Delta }f_{n1}\delta _{\{0\}})f$$

where $``$ is the convolution on $`^2`$. Since $`\mathrm{\Delta }(gh)=g\mathrm{\Delta }h`$ we get

$$\mathrm{\Delta }f_n\delta _{\{0\}}=(\mathrm{\Delta }f_{n1}\delta _{\{0\}})(\mathrm{\Delta }f\delta _{\{0\}})$$

and hence (7) gives that the $`f_n`$’s converge333The convergence rate is obviously exponential, but one might also define $`f_n=f_{n1}(\mathrm{\Delta }f_{n1}\delta _{\{0\}})f_{n1}`$ and get double exponential convergence. in the $`L^1`$ norm to an $`a`$ which satisfies $`\mathrm{\Delta }a=\delta _{\{0\}}`$ and $`|a|C\mathrm{log}|z|+C`$ so it is the harmonic potential up to an additive constant. Furthermore, the fact that $`\mathrm{\Delta }f\delta _{\{0\}}=O(|z|^4)`$ can be used to derive similar estimates for $`a`$. However, the best approximation we were able to get from that method is

$$a(z)=\frac{2}{\pi }\mathrm{log}|z|+\lambda \mathrm{Re}\frac{z^4}{6\pi |z|^6}+O(|z|^4\mathrm{log}|z|)$$
(8)
and getting this approximation was not significantly easier than the brute force approach. Moreover, this approach did not allow to get actual values for any of the constants — only that some constants exist, and their values were derived from (6).

We have described the contents of this paper except for the last section. That section contains the results of some computer-aided simulations and a calculation of an *explicit* constant in the $`O()`$ in (2).

We wish to thank Greg Lawler for some useful remarks.

### 1.1. Standard definitions

A lattice is a discrete additive subgroup of $`^n`$. The dimension of a lattice $`Z`$ (denoted by $`dimZ`$) is the dimension of the linear span of $`Z`$ as a linear subspace of $`^n`$. A basis for a $`d`$-dimensional lattice is a set $`\{e_1,\mathrm{},e_d\}`$ such that $`Z=e_1+\mathrm{}+e_d`$. $`Z^{}`$ is called a sublattice of $`Z`$ if it is a subgroup of $`Z`$, and its index is the size of $`Z/Z^{}`$. The volume of (the cell of) a lattice $`Z`$, denoted by $`\mathrm{vol}Z`$, is the volume of $`^d/Z`$ (the discreteness of $`Z`$ allows us to select a measurable set of representatives). Alternatively it can be defined as

$$\underset{r\mathrm{}}{lim}\frac{|B(r,0)|}{\mathrm{\#}(B(r,0)Z)}$$

where $`B(r,0)`$ is a ball of radius $`r`$ around $`0`$ and $`|B(r,0)|`$ is its volume.

A random walk $`R`$ on a lattice $`Z`$ is a probabilistic process $`\{R_0,R_1,\mathrm{}\}`$, $`R_iZ`$, such that $`R_0=0`$ and $`R_iR_{i1}`$ has a distribution independent of $`i`$, and such that for every $`zZ`$ there exists some $`n`$ such that $`(z=R_n)>0`$.444This definition is a bit restrictive. For example, a $`1`$-dimensional random walk which always goes to the right is not a “random walk on a lattice”. It is bounded if $`R_iR_{i1}`$ is bounded. We will typically confuse the process $`R`$ with the distribution of any step (e.g. with the distribution of $`R_1`$) so that we can use notations such as $`𝔼R`$ comfortably. The dimension of $`R`$ is the dimension of the lattice $`Z`$, and is denoted by $`dimR`$.

The drift of a random walk $`R`$ is $`𝔼R`$. The random walk is balanced if its drift is zero.

The discrete Laplacian $`\mathrm{\Delta }`$ of a random walk $`R`$ is an operator on functions on the lattice $`Z`$ defined by

$$(\mathrm{\Delta }f)(z)=𝔼f(z+R_1)f(z).$$

Functions with $`\mathrm{\Delta }f0`$ are called (discretely) harmonic.

The continuous Laplacian on $`^d`$, $`\frac{^2}{z_1^2}+\mathrm{}+\frac{^2}{z_d^2}`$ will be denoted by $`\mathrm{\Delta }_C`$. Functions with $`\mathrm{\Delta }_Cf0`$ will be called continuously harmonic.

$`v`$ will always refer to the $`L^2`$ norm. Similarly, for a matrix $`A`$, $`A`$ refers to its norm as an operator from $`L^2`$ to $`L^2`$.

## 2. The direct approach

The aim of this section is to prove the following theorem.

###### Theorem 1.

Let $`a`$ be the harmonic potential of a $`d`$-dimensional balanced bounded random walk on a lattice $`Z`$. Let $`K>0`$ be some integer. Then there exists a constant $`\tau `$, a linear map $`A`$, a polynomial $`Q=Q_K`$ in $`d`$ variables and an integer $`L=L_K`$ such that

$$a(z)=\tau \mu _d(Az)+\frac{Q(z)}{Az^L}+O(z^K)$$
(9)
where

$$\mu _d(z):=\{\begin{array}{cc}z^{2d}\hfill & d2\hfill \\ \mathrm{log}z\hfill & d=2\hfill \end{array}.$$

$`\tau `$ and $`A`$ have explicit formulas — see the comments after the end of the proof (page 2 below). The expression $`Q(z)/Az^L`$ is $`\lambda +O(Az^{1d})`$ (this is classical, of course, but will also follow from the proof below). The number $`\lambda `$ is of course also walk-specific.

In the case $`d=1`$ the equation $`\mathrm{\Delta }a=\delta _{\{0\}}`$ deteriorates into a simple recursion formula and the theorem is nothing but an exercise (and $`Q\lambda `$). Thus we concentrate on the case $`d>1`$.

###### Lemma 1.

Let $`p_1,\mathrm{},p_k]0,1[`$ such that $`p_i=1`$. Let $`K`$. Then there exists a $`Q=Q_{p,K}`$ such that

$`B_p(n,p_1n+w_1,\mathrm{},p_kn+w_k)=`$ (10)
$`{\displaystyle \frac{1}{(2\pi n)^{(k1)/2}\sqrt{p_1\mathrm{}p_k}}}e^{{\scriptscriptstyle w_i^2/2p_in}}Q(n,w_1,\mathrm{},w_k)+O(n^K)`$
for every $`w_1,\mathrm{},w_k`$ such that $`w_i=0`$. $`B`$ here is the standard multinom

$$B_p(n,m_1,\mathrm{},m_k)=\frac{n!}{m_1!\mathrm{}m_k!}p_1^{m_1}\mathrm{}p_k^{m_k},m_1+\mathrm{}+m_k=n$$

and $`Q`$ satisfies that $`n^LQ`$ is a polynomial for some $`L=L_K`$ and $`Q(n,w)=1+o(1)`$ uniformly in $`wC\sqrt{n\mathrm{log}n}`$.

The proof of theorem 1 is replete with these rational functions $`Q`$ obeying the condition above. Notice that this simply means that

$$Q(n,w_1,\mathrm{},w_k)=\underset{i,\stackrel{}{j}}{}c_{i\stackrel{}{j}}n^iw_1^{j_1}\mathrm{}w_k^{j_k}$$

such that $`c_{0\stackrel{}{0}}=1`$ and otherwise $`c_{i\stackrel{}{j}}0`$ only if $`Li<0`$ and $`0j_1+\mathrm{}+j_k<2i`$. It turns out (and this can be deduced from the proofs below with some care) that $`L_K=4K`$ is always sufficient. However we will have no use for this fact.

###### Proof of lemma 1.

Consider Stirling’s series,

$$n!=\frac{n^n}{e^n}\sqrt{2\pi n}\mathrm{exp}\left(Q_1(n)+O(n^K)\right)$$

where $`n^{K1}Q_1`$ is some polynomial with $`Q_1(n)=o(1)`$ as $`n\mathrm{}`$. This gives

$`B_p(n,p_1n_1+w_1,\mathrm{},p_kn_k+w_k)={\displaystyle \frac{1}{(2\pi n)^{(k1)/2}\sqrt{p_1\mathrm{}p_k}}}{\displaystyle \underset{i=1}{\overset{k}{}}}{\displaystyle \frac{1}{\sqrt{1+\frac{w_i}{p_in}}}}`$
$`\mathrm{exp}\left({\displaystyle \underset{i=1}{\overset{k}{}}}(p_in+w_i)\mathrm{log}\left(1+{\displaystyle \frac{w_i}{p_in}}\right)+Q_1(n){\displaystyle \underset{i=1}{\overset{k}{}}}Q_1(p_in+w_i)+O(n^K)\right)`$
Assume for the moment that $`w<C\sqrt{n\mathrm{log}n}`$ where $`C`$ will be fixed later. In this case we can insert monomials such as $`\frac{w_i^{2K+1}}{n^{2K+1}}`$ into $`O(n^K)`$. Thus, each of the factors $`\frac{1}{\sqrt{1+\frac{w_i}{p_in}}}`$ equals $`Q_2(n,w_i)+O(n^K)`$ with $`n^{2K}Q_2`$ some polynomial with $`Q_2=1+o(1)`$ ($`Q_2`$ depends on $`p_i`$, of course). The term inside the exponent evaluates to $`\frac{w_i^2}{2p_in}+Q_3(n,\stackrel{}{w})+O(n^K)`$ with $`n^{2K}Q_3`$ a polynomial satisfying $`Q_3=o(1)`$ uniformly in $`w<C\sqrt{n\mathrm{log}n}`$ as $`n\mathrm{}`$. This proves the case $`w<C\sqrt{n\mathrm{log}n}`$. In the case $`w>C\sqrt{n\mathrm{log}n}`$ both sides of (10) are $`O(n^K)`$ which follows easily from the fact that in this case, for $`C`$ sufficiently large

$$\mathrm{exp}\left(\frac{w_i^2}{2p_in}\right)\mathrm{exp}\left(C\mathrm{log}n\right)Cn^K.\mathit{}$$

###### Lemma 2.

Let $`v^d`$ and $`A\mathrm{GL}(^d)`$. Then

$$\underset{z^d}{}\mathrm{exp}\left(\frac{1}{n}A(z+v)^2\right)=\frac{(\pi n)^{d/2}}{|A|}+O(e^{cn})$$

where the $`O()`$ is uniform in $`v`$. ($`|A|`$ here is the determinant of $`A`$).

Furthermore, if $`P`$ is any polynomial and $`y^d`$ then

$$\underset{z^d+y}{}\mathrm{exp}\left(\frac{1}{n}A(z+v)^2\right)P(z)=n^{d/2}Q(n)+O(e^{cn}v^{\mathrm{deg}P})$$
(11)
where $`Q`$ is a polynomial. The coefficients of $`Q`$ are fixed polynomial functions (depending on $`A`$) of the coefficients of $`P`$ and of $`v`$. $`Q`$ is independent of $`y`$ and the $`O()`$ is uniform in $`v`$ and $`y`$.

The use of the parameter $`n`$ in the formulation of the lemma is a bit artificial. Later, when we shall use the lemma, $`n`$ will be an integer number (the number of steps), but this fact is not needed for the formulation or proof of the lemma.

###### Proof.

We use the $`d`$-dimensional Poisson summation formula, $`f(\stackrel{}{n})=\widehat{f}(2\pi \stackrel{}{n})`$ (the Fourier transform being defined by $`\widehat{f}(\xi )=f(x)e^{ix,\xi }𝑑x`$) for $`f(z)=\mathrm{exp}\left(\frac{1}{n}A(z+v)^2\right)(z+v)^\stackrel{}{r}`$ where $`\stackrel{}{r}=(r_1,\mathrm{},r_d)\{0,1,\mathrm{}\}^d`$ and $`(z+v)^\stackrel{}{r}`$ is short for $`(z_1+v_1)^{r_1}\mathrm{}(z_d+v_d)^{r_d}`$. A simple calculation shows that

$$\widehat{f}(\xi )=e^{iv,\xi }i^r\frac{d^r}{d\xi ^\stackrel{}{r}}\frac{(\pi n)^{d/2}}{|A|}\mathrm{exp}\left(\frac{1}{4}n\left(A^1\right)^{}\xi ^2\right)$$
(12)
where $`r=r_1+\mathrm{}+r_d`$, and as usual $`\frac{d^r}{d\xi ^\stackrel{}{r}}=\frac{d^{r_1}}{d\xi _1^{r_1}}\mathrm{}\frac{d^{r_d}}{d\xi _d^{r_d}}`$. This means that as $`n\mathrm{}`$ we have

$$\underset{\xi 2\pi ^d}{}\widehat{f}(\xi )=\widehat{f}(0)+O\left(n^{r+d/2}A^1^{d+2r}e^{n/4A}\right)$$

(note that $`|A|^1A^1^d`$). The exponent inside the $`O()`$ kills of course all the other factors. The case $`\stackrel{}{r}=\stackrel{}{0}`$ immediately gives the first part of the lemma. For the second part, put $`\xi =0`$ in (12) and get that $`\widehat{f}(0)0`$ only when $`r`$ is even, and that $`\widehat{f}(0)=n^{d/2}Q_\stackrel{}{r}(n)`$ where $`Q_\stackrel{}{r}`$ depends only on $`A`$. Writing $`P=_\stackrel{}{r}c_\stackrel{}{r}(z+v)^\stackrel{}{r}`$ we get $`Q=_\stackrel{}{r}c_\stackrel{}{r}Q_\stackrel{}{r}`$ which concludes the lemma: the fact that $`c_\stackrel{}{r}`$ are polynomial in the coefficients of $`P`$ and in $`v`$ gives the same for the coefficients of $`Q`$; the fact that $`c_\stackrel{}{r}=O(v^{\mathrm{deg}P})`$ gives the same in the $`O()`$ in (11). The result does not depend on $`y`$ (except in the $`O()`$ error) since $`f(n)`$ does not depend on $`v`$ — in effect the dependence on $`v`$ appears only via the dependence of $`c_\stackrel{}{r}`$ on $`v`$. ∎

###### Lemma 3.

Let $`R`$ be a $`d`$-dimensional balanced bounded random walk on the lattice $`Z`$. Let $`K`$. Then there exists a sublattice $`Y\times Z`$, a constant $`\tau ^{}`$, a quadratic function $`B`$ and a $`Q`$ such that

$$P_n(z)=\{\begin{array}{cc}\tau ^{}n^{d/2}e^{B(z)/n}Q(n,z)+O(n^K)\hfill & (n,z)Y\hfill \\ 0\hfill & (n,z)Y\hfill \end{array}.$$
(13)
with $`n^LQ`$ a polynomial in $`n`$ and in the coordinates of $`z`$, and $`Q(n,z)=1+o(1)`$ uniformly in $`zC\sqrt{n\mathrm{log}n}`$.

Assume for simplicity that $`Z^d`$. A simple comparison of the expectations of the left and right side of (13) shows that $`B(z)=\frac{1}{2}M^1z,z`$, where $`M`$ is the correlation matrix of $`R`$, $`M_{ij}=𝔼R,e_iR,e_j`$, where $`e_i`$ are unit vectors. Notice that the expression $`M^1z,z`$ is independent of the coordinate system chosen.

The lattice $`Y`$ is clearly $`d+1`$ dimensional since otherwise we would have $`dimR<d`$. Thus $`YZ`$ is $`d`$-dimensional and $`\mathrm{vol}YZ`$ is finite. This, with the requirement $`P_n=1`$ allows to calculate $`\tau ^{}`$ and get

$$\tau ^{}=\frac{\mathrm{vol}(YZ)}{(2\pi )^{d/2}|detM|^{1/2}}.$$
(14)
We will not be using these equalities in the proof, though.

###### Proof.

By applying a linear transformation we may assume $`R`$ is a random walk on $`^d`$. Assume $`R`$ moves from $`z`$ to $`z+v_i`$ with probability $`p_i`$ for $`i=1,\mathrm{},N+1`$. Let

$$W=\{(n,w_1,\mathrm{},w_{N+1}):p_in+w_ii\text{ and }w_i=0\}$$

which is an $`(N+1)`$-dimensional lattice. For an $`n`$ let $`k_i`$ be the number of times (between $`0`$ and $`n1`$) that $`R`$ moved from $`z`$ to $`z+v_i`$. We use lemma 1 with some $`K_2`$ which will be fixed later. This gives us

$$(k_i=np_i+w_i)=Cn^{N/2}\mathrm{exp}\left(\frac{1}{n}B_2(w)\right)Q_2(n,w_1,\mathrm{},w_N)+O(n^{K_2})$$
(15)
where $`B_2`$ is a strictly positive quadratic form, $`C`$ is some constant and $`n^LQ_2`$ is a polynomial with $`Q_2=1+o(1)`$ uniformly for $`wC\sqrt{n\mathrm{log}n}`$. Now, for a point $`z^d`$ we can write

$`(R_n=z)`$ $`={\displaystyle \underset{wW(n,z)}{}}(k_i=np_i+w_i)`$ (16)
$`W(n,z)`$ $`:=\{w:(n,w)W{\displaystyle w_iv_i}=z\}.`$
Since $`W(n,z)`$ are intersections of a lattice with a translated subspace, we get that $`W(n,z)=\mathrm{}`$ for $`(n,z)`$ outside some lattice $`Y`$, and for $`(n,z)Y`$ we get that $`W(n,z)`$ are translations of one fixed lattice $`X`$. Since we know $`Y`$ is $`d+1`$ dimensional (see the comment just after the statement of this lemma) we get that $`X`$ is $`(Nd)`$-dimensional. The calculation of the sum in (16) will be done of course with lemma 2 but we need some preparations first. By applying a linear transformation to the $`N`$-dimensional space of $`w`$’s we may assume that

$$W(n,z)=z_1e_1+\mathrm{}+z_de_d+e_{d+1}+\mathrm{}+e_N+y_n$$

where $`y_n\mathrm{span}\{e_{d+1},\mathrm{},e_N\}`$. Combining (15), (16) and the linear transformation we get

$$(R_n=z)=\underset{wW(n,z)}{}Cn^{N/2}\mathrm{exp}\left(\frac{1}{n}B_3(w)\right)Q_3(n,w)+O(n^{K_3})$$
(17)
where $`K_3=K_2N+d`$. This last error estimate follows easily from the following obvious estimates:

$$\underset{\begin{array}{c}wW(n,z)\\ p_in+w_i[0,n]\end{array}}{}O(n^{K_2})=O(n^{K_3}),\underset{\begin{array}{c}wW(n,z)\\ p_in+w_i[0,n]\end{array}}{}e^{B_3(w)/n}Q_3(n,w)=O(e^{cn})$$

From now on it would be easier to replace $`B_3(w)`$ with $`Aw^2`$ for some $`A\mathrm{GL}(^N)`$ which is possible since $`B_3`$ is also strictly positive. We denote $`x=(0,\mathrm{},0,w_{d+1},\mathrm{},w_N)`$ so that $`w=x+z`$ and we can write

$$Aw^2=Ax^2+Az^2+2x,A^{}Az.$$

Since $`A^{}A`$ is strictly positive we get that its lower right $`Nd`$ minor is also strictly positive and thus invertible. Therefore there exists some $`v=(0,\mathrm{},0,v_{d+1},\mathrm{},v_N)`$ such that

$$(A^{}Av)_i=(A^{}Az)_i,i=d+1,\mathrm{},N.$$

This of course implies $`x,A^{}Av=x,A^{}Az`$ and then

$$Aw^2=Az^2Av^2+A(x+v)^2.$$
(18)
The importance of (18) is that its left part, $`Az^2Av^2`$, is constant for $`wW(n,z)`$, while $`\mathrm{exp}\left(\frac{1}{n}A(x+v)^2\right)`$ can be summed by using lemma 2 on the last $`Nd`$ coordinates. Notice that it depends on the $`z`$’s only via $`v`$, which is linear in $`z`$. In other words, we can now define our quadratic form $`B`$:

$$B(z):=Az^2Av(z)^2$$

and apply lemma 2. We get

$`(\text{17})`$ $`=Cn^{N/2}e^{B(z)/n}{\displaystyle \underset{x^{Nd}+y_n}{}}\mathrm{exp}\left({\displaystyle \frac{1}{n}}A(x+v)^2\right)Q_3(n,z,x)+O(n^{K_3})`$
$`=Cn^{d/2}e^{B(z)/n}Q(n,z)+O(n^{K_3})+O(e^{cn}v^{\mathrm{deg}Q_3})`$
We may already remark that picking $`K_2=K+Nd`$ will give us a precision of $`O(n^K)`$ as required. Obviously $`O(e^{cn}v^{\mathrm{deg}Q_3})=O(n^K)`$ for any $`K`$ as $`v`$ is linear in $`z`$ and $`zC\sqrt{n\mathrm{log}n}`$. To say something about $`Q`$, we need to consider $`Q_3`$ as a polynomial in $`x`$ with coefficients which are polynomials in $`z`$ and negative powers in $`n`$, i.e.

$$Q_3=\underset{\stackrel{}{r}}{}c_\stackrel{}{r}(n,z)x^\stackrel{}{r}$$

where $`n^Lc_\stackrel{}{r}`$ are polynomials. Lemma 2 tells us that $`n^LQ`$ is a polynomial in $`n`$ and that its coefficients are polynomial in the coefficients of $`Q_3`$ and in $`v`$ (which is linear in $`z`$). Thus we get that $`n^LQ`$ is polynomial in both $`n`$ and $`z`$. This concludes the lemma: the only claim we haven’t shown is that $`Q=1+o(1)`$ which can be deduced, for example, from the central limit theorem555Actually, it is also possible to deduce it from the arguments above with a little care as to the exact dependence between $`P`$, $`Q`$ and $`v`$ in lemma 2.. ∎

###### Remark.

For specific walks it is possible to prove lemma 3 using much simpler techniques. For example, for the simple random walk on $`^2`$, we have

$$P_n(x+iy)=4^n\underset{k=0}{\overset{n}{}}\left(\genfrac{}{}{0pt}{}{n}{k}\right)\left(\genfrac{}{}{0pt}{}{k}{(kx)/2}\right)\left(\genfrac{}{}{0pt}{}{nk}{(nky)/2}\right)$$

where we make the notational convention that $`\left(\genfrac{}{}{0pt}{}{n}{\alpha }\right)=0`$ if $`\alpha `$ is not an integer. This formula allows to get lemma 3 from lemma 1 more-or-less directly, with no need to go through lemma 2. This works for other walks as well. However it does not seem to work for the very interesting case of the random $`6`$-walk on the triangular lattice, a walk which can be considered no less natural than the random walk on $`^2`$. We mean here a two dimensional random walk that goes at each step with probability $`\frac{1}{6}`$ from $`z`$ to one of the points $`z+\omega ^i`$, $`i=0,\mathrm{},5`$ where $`\omega =\frac{1}{2}+i\frac{\sqrt{3}}{2}=\sqrt[6]{1}`$.

###### Lemma 4.

Let $`q`$ and $`1jq`$. For every $`s>1`$ and for $`r\mathrm{}`$ we have the superpolynomial estimate

$$F_{s,j,q}(r):=\underset{\begin{array}{c}nj(q)\\ nj\end{array}}{}n^se^{r/n}=\frac{\mathrm{\Gamma }(s1)}{qr^{s1}}+O(r^K)K;$$
(19)
while for $`s=1`$ we have

$$F_{1,j,q}(r):=\underset{\begin{array}{c}nj(q)\\ nj\end{array}}{}\frac{1}{n}\left(e^{r/n}1\right)=q^1\mathrm{log}r+C+O(r^K)K.$$

###### Proof.

Again we use the (one dimensional) Poisson summation formula, this time for the function

$$f(x)=\{\begin{array}{cc}x^se^{1/x}\hfill & x>0\hfill \\ 0\hfill & x0\hfill \end{array}.$$

Simple Fourier operations show that

$$F_{s,j,q}=\frac{1}{qr^{s1}}\underset{\xi }{}e^{2\pi ij\xi /q}\widehat{f}\left(2\pi \frac{r\xi }{q}\right).$$

Since $`fC^{\mathrm{}}`$ we get that $`\widehat{f}(\xi )=O(\xi ^K)`$ for all $`K`$ and hence $`F=\frac{r^{1s}}{q}\widehat{f}(0)+O(r^K)`$. This concludes the case $`s>1`$. For the case $`s=1`$ we notice that the $`O()`$ in (19) is uniform in $`s(1,2]`$ since it depends only on the $`L^1`$ norm of the $`K`$’th derivative of $`f`$. Thus we need only to calculate

$$\underset{s1}{lim}\frac{\mathrm{\Gamma }(s1)}{r^{s1}}\zeta (s)$$

where $`\zeta `$ is Riemann’s function. Writing $`\mathrm{\Gamma }=\frac{1}{s}+C+o(1)`$ and $`\zeta (s)=\frac{1}{s1}+C+o(1)`$ and a little use of l’Hôpital rule will give $`\mathrm{log}r+C`$ and conclude the lemma. ∎

###### Remark.

The exact asymptotics of the error are harder to figure out. This is of course a question about the asymptotics of $`\widehat{f}`$. A much better estimate than $`\widehat{f}(\xi )\xi ^K`$ can be had by changing the path of integration of $`x^se^{1/x+ix\xi }`$ from $`[0,\mathrm{})`$ to $`[0,1+i][1+i,\mathrm{}+i)`$ or to its conjugate (depending on the sign of $`\xi `$). On the segment $`[0,\xi ^{1/2}(1+i)]`$ we have $`e^{1/x}=O(e^{c\sqrt{\xi }})`$. On the reminder we have $`e^{ix\xi }=O(e^{c\sqrt{\xi }})`$. So we get $`\widehat{f}(\xi )=O(e^{c\sqrt{\xi }})`$.

In the other direction, there are known precise results connecting the “depth” of the zero of $`f`$ at $`0`$, i.e. the rate of convergence of $`f`$ to zero, with the property that $`\widehat{f}`$ is in some weighted $`H^p`$ space. See theorem 1.1 in \[MS02\]. We omit the details. This shows, for example, that $`\widehat{f}(\xi )O(e^{\xi ^{1/2+ϵ}})`$ for any $`ϵ>0`$.

###### Proof of theorem 1.

Let $`B`$, $`\tau ^{}`$ and $`Y`$ be the quadratic form, constant and lattice given for $`R`$ and $`K+1`$ by lemma 3. Denoting

$$\stackrel{~}{a}(z):=\underset{n:(n,z)Y}{}\tau ^{}n^{d/2}P_n(z)$$

we see that $`\stackrel{~}{a}(z)a(z)`$ is a constant (indeed, the precaution of adding $`\tau ^{}n^{d/2}`$ is only necessary in the case $`d=2`$ in order to make the sum converge). The fact that $`R`$ is bounded shows that $`(R_n=z)`$ is zero for $`n<cz`$ and of course $`_{n<cz}e^{B(z)/n}=O(e^{cz})`$. Therefore we can use (13) to get

$$\stackrel{~}{a}(z)=\underset{n:(n,z)Y}{}\tau ^{}n^{d/2}\left(e^{B(z)/n}Q_2(n,z)1\right)+\underset{n>cz}{}O(n^{K1})+O(e^{cz})$$

Lemma 4 allows to calculate the left sum, while the right one is obviously $`O(z^K)`$. This concludes the proof of theorem 1. ∎

We see now that the quadratic form of theorem 1 is the same as that of lemma 3, i.e. $`B(z)=\frac{1}{2}M^1z,z`$ where $`M`$ is the correlation matrix of $`R`$. Thus $`A=M^{1/2}`$ — of course, this is not uniquely defined, but the quadratic form $`Az^2`$ is. As for $`\tau `$, the factor $`\mathrm{vol}YZ`$ of lemma 3 with the term $`\frac{1}{q}`$ of lemma 4 gives $`\mathrm{vol}Z`$ (this is easy to see) and we get

$$\tau =\frac{\mathrm{vol}Z}{(2\pi )^{d/2}|detM|^{1/2}}\{\begin{array}{cc}2\hfill & d=2\hfill \\ \mathrm{\Gamma }(\frac{1}{2}d1)\hfill & d2\hfill \end{array}$$

(the factor $`2`$ in dimension $`2`$ comes from the fact that we formulated $`\mu _2=\mathrm{log}z`$ rather than $`\mathrm{log}z^2`$).

## 3. The differential approach

We shall demonstrate the differential approach by proving a few results that would seem quite difficult using the direct approach only.

For simplicity of notation we shall always assume that the correlation matrix of $`R`$ is constant. Thus we call a random walk *spherical* if it is bounded, balanced, and $`𝔼R,e_iR,e_j=C\delta (i,j)`$ (again, $`e_i`$ are unit vectors).

###### Theorem 2.

Let $`R`$ be a $`2`$-dimensional spherical random walk on $``$. Then

$$a(z)=\alpha \mathrm{log}|z|+\lambda +\mathrm{Re}\left(\alpha _{lk}\frac{z^l}{|z|^k}\right)+O(|z|^K)$$
(20)
and $`\alpha _{lk}=0`$ whenever $`l<\frac{1}{2}k`$.

The sum here is over a finite collection of $`l0`$ and $`k0`$ with the condition $`1lkK`$.

###### Proof.

The complex representation (20) is nothing but a restatement of theorem 1: writing $`x=\frac{1}{2}(z+\overline{z})`$ and $`y=\frac{1}{2i}(z\overline{z})`$ we get a representation in the form

$$\alpha _{ijk}\frac{z^i\overline{z}^j}{|z|^k}.$$

but of course if both $`i0`$ and $`j0`$ then we can cancel one of them with the $`|z|^k`$ factor. The fact that the end result is real means that $`\alpha _{0jk}=\overline{\alpha _{j0k}}`$ and we get (20). Thus what we need to prove is the fact that $`\alpha _{lk}=0`$ whenever $`l<\frac{k}{2}`$. We change the representation again, writing

$$a(z)=\beta \mathrm{log}(z\overline{z})+\underset{k+l>K}{}\beta _{kl}z^k\overline{z}^l+O(|z|^K)$$
(21)
and now we have to show that $`\beta _{kl}=0`$ whenever both $`k`$ and $`l`$ are negative (of course, $`\beta _{kl}=0`$ when $`k+l>0`$). Now, the first thing to notice is that we may differentiate (21) formally. If $`R`$ moves to $`v_i`$ with probability $`p_i`$ then

$`(\mathrm{\Delta }f)(z)`$ $`={\displaystyle p_if(z+v_i)}f(z)=`$
$`={\displaystyle \underset{n=1}{\overset{K1}{}}}{\displaystyle \frac{1}{n!}}{\displaystyle \underset{i}{}}p_i{\displaystyle \frac{^n}{^nv_i}}f(z)+O\left({\displaystyle \frac{1}{K!}}{\displaystyle \underset{i}{}}p_i\underset{\eta [0,1]}{\mathrm{max}}\left|{\displaystyle \frac{^K}{^Kv_i}}f(z+\eta v_i)\right|\right)`$
where we used the one-dimensional Taylor expansion for the functions $`f(z+tv_i)`$ for each $`i`$ (the constant in the $`O()`$ above is of course simply $`1`$). Denote the $`n`$th differential operator in the sum by $`D_n`$ and denote the maximum inside the $`O()`$ by $`M_K`$. Clearly,

$$D_nz^k\overline{z}^l=O(|z|^{k+ln})M_nz^k\overline{z}^l=O(|z|^{k+ln}).$$
(22)
Now, $`\mathrm{\Delta }a=0`$ with (21) and (22) show that

$$\left(\underset{n=1}{\overset{K1}{}}D_k\right)\left(\beta \mathrm{log}z\overline{z}+\underset{k+l>K}{}\beta _{kl}z^k\overline{z}^l\right)=O(|z|^K).$$
(23)
This is exactly what we mean when we say that (21) can be differentiated formally.

Now to the actual calculation. Notice that $`D_1=0`$ (because $`R`$ is balanced) and that $`D_2`$ is, up to multiplication with a constant, the continuous Laplacian, because $`R`$ is spherical. Therefore

$$D_2(z^k\overline{z}^l)=2\gamma klz^{k1}\overline{z}^{l1}$$
(24)
for some constant $`\gamma 0`$. Assume that we know $`\beta _{kl}`$ for all $`k+l>m`$, and that $`\beta _{kl}=0`$ if both are negative. We write

$$\left(\underset{n=1}{\overset{m}{}}D_n\right)\left(\beta \mathrm{log}z\overline{z}+\underset{k+l>m}{}\beta _{kl}z^k\overline{z}^l\right)=\delta _{kl}z^k\overline{z}^l$$

we get that $`\delta _{kl}=0`$ if both $`k`$ and $`l`$ are negative. Moreover, (23) for $`K>m+2`$ gives equations for $`\beta _{kl}`$ with $`k+l=m`$, namely $`\beta _{kl}=\frac{1}{2\gamma kl}\delta _{k1,l1}`$. This shows inductively that $`\beta _{kl}=0`$ when they are both negative, and the theorem is proven. ∎

The case $`l=0`$ in theorem 2 has some potential uses for precise estimates of hitting probabilities, so it seems worth of special mention:

###### Corollary.

For every spherical random walk on $`^2`$, the high-order expansion of the harmonic potential $`a`$, considered as a function of a continuous variable, satisfies the superpolynomial estimate

$$_0^{2\pi }a(Re^{it})𝑑t=A\mathrm{log}R+B+O(R^K)K.$$

###### Remark.

All the above discussion has analogues in dimension $`>2`$, using the notions of spherical harmonics. Spherical harmonics are functions on $`S^{d1}`$ which are eigenvalues of the angular part of the continuous Laplacian (which is the angular momentum squared operator). It turns out that each is a restriction to $`S^{d1}`$ of a polynomial on $`^d`$ and we shall use this representation. For a general introduction to the topic, see \[SW71, sect. 4.2\].

* In two dimensions, a factor $`P/|z|^\alpha `$ (where $`P`$ is a homogeneous polynomial) may be decomposed into a sum of the terms $`z^k\overline{z}^l`$, $`k+l=\mathrm{deg}P\alpha `$. The high dimension analogue is a unique decomposition into a sum of terms $`Q/z^\beta `$ where $`Q`$ is a homogeneous (continuously) harmonic polynomial, $`\mathrm{deg}Q\beta =\mathrm{deg}P\alpha `$.
* The analogue of the fact that only $`z^k`$ and $`\overline{z}^l`$ are harmonic is the fact that a term $`Q/z^\beta `$ is harmonic if and only if $`\beta =2\mathrm{deg}Q+d2`$. This follows from the following general formula:

$$\mathrm{\Delta }_C\left(\frac{Q}{z^\beta }\right)=\frac{z^2\mathrm{\Delta }_CQ+(\beta ^2\beta (2\mathrm{deg}Q+d2))Q}{z^{\beta +2}}$$
(25)
and from the harmonicity of $`Q`$ (remember that $`\mathrm{\Delta }_C`$ is the continuous Laplacian). (25) also gives the analogue of (24).
* The analogue of “$`k`$ and $`l`$ are both negative” is “$`\beta <2\mathrm{deg}Q+d2`$”, or, equivalently, “$`\mathrm{\Delta }_C^n(Q/z^\beta )`$ is never zero, for any $`n`$”.
* The analogue of theorem 2 is the following: if $`R`$ is a spherical random walk and

$$a(z)=\frac{\tau }{z^{d2}}+\lambda +\frac{Q(z)}{z^\beta }+O(z^K)$$

then

$$\mathrm{\Delta }_C^K\left(\frac{Q(z)}{z^\beta }\right)=0$$

We omit the proof, which is a similar induction.
* The corollary to theorem 2 is also true in higher dimensions. This follows from the fact that the terms $`Q/z^\beta `$ for different $`\beta `$’s are orthogonal in the sphere measure:

$$_{S^{d1}}Q_1Q_2=0Q_1,Q_2\text{ harmonic, }\mathrm{deg}Q_1\mathrm{deg}Q_2$$

and in particular $`Q`$ is orthogonal to $`1`$ if $`\mathrm{deg}Q>0`$.

###### Theorem 3.

The coefficients of the high-order expansion of the harmonic potential $`a`$ of the regular random walk on $`^2`$ can be calculated from the differential equations and (6).

###### Proof.

There is almost nothing to prove here. We saw during the proof of the previous theorem that (23) allows to calculate all coefficients except those where $`k=0`$ or $`l=0`$ (i.e. the coefficients of the harmonic summands). There is only one of these (and its conjugate) in every level $`k+l=n`$. Therefore the coefficient of $`z^n`$ can be determined from the coefficients of $`z^{n1}\overline{z},z^{n2}\overline{z}^2,\mathrm{}`$ and from (6). ∎

The proof of theorem 3 is also the second algorithm for getting (4).

###### Theorem 4.

1. The term of order $`k`$ in any walk is $`P/z^\alpha `$ where $`\alpha 4k+2d`$ and $`P`$ is a polynomial of degree $`\alpha k3k+2d`$.
2. If $`D_30`$, in particular if $`R`$ is reversible666A reversible (or symmetric) walk is a random walk satisfying $`R=R`$. In graph language, a reversible random walk is a walk on a simple graph, while a non-reversible walk is a walk on a directed graph. Many of the claims of this paper need slight alterations for a non-reversible walk. For example, we have $`\stackrel{~}{\mathrm{\Delta }}a=\delta _{\{0\}}`$ where $`\stackrel{~}{\mathrm{\Delta }}`$ is the Laplacian of the time-reversed walk., then the term of order $`k`$ is $`P/z^\alpha `$ where $`\alpha 3k+2d`$ and $`P`$ is a polynomial of degree $`\alpha k2k+2d`$.

###### Proof.

Take as an example clause 2. We may assume that the walk is spherical. In two dimensions, the theorem is equivalent to showing that only terms of the form $`z^k\overline{z}^l`$, $`0k\frac{1}{3}|l|`$ (or $`0l\frac{1}{3}|k|`$) appear. On the one hand, $`k`$ can only increase (by $`1`$) in the step of taking $`D_2^1`$ (see the proof of theorem 2), and this operation also increases $`l`$ by $`1`$. On the other hand, the step of taking $`D_4+D_5+\mathrm{}`$ decreases $`l`$ by at least $`4`$, and the theorem follows. In higher dimensions the same argument works using the decomposition into spherical harmonics discussed before theorem 3. ∎

We note that clause 1 of theorem 4 is already mentioned in \[FU96, U98\].

Theorem 4 can be generalized to cases where higher symmetries exist. For example, for the $`6`$-walk on the triangular lattice, $`D_3D_50`$ and $`D_4=\frac{9}{4}\mathrm{\Delta }_C^2`$ which allow to get that the term of order $`|z|^k`$ is $`P/|z|^\nu `$ where $`\nu \frac{5}{2}k`$ and $`\mathrm{deg}P\frac{3}{2}k`$.

###### Theorem 5.

The harmonic potential of a reversible random walk contains only terms of even order in even dimension and only terms of odd order in odd dimension.

###### Proof.

Since $`D_{2k+1}0`$ for all $`k`$, we only need to handle the terms we cannot calculate, namely the (continuously) harmonic terms. However, a harmonic term of these orders must be $`P/z^\alpha `$ with $`\mathrm{deg}P`$ odd, and since $`a(z)=a(z)`$ such terms cannot appear in the expansion of $`a`$. ∎

## 4. Explicit constants

This note is a spin-off from a different work \[K\] where it was necessary to have an *explicit* value for the constants in the $`O()`$ in (2). How should one go about to calculate something like that? Generally, one would try to get some explicit constant in the $`O`$ of an approximation of one order higher, and then use the McCrea-Whipple algorithm to check the first few values.

Which approach would be best to get an explicit constant in the $`O()`$? The differential approach doesn’t seem to give any explicit constants whatsoever. Therefore it was necessary to use the direct approach. A program was written to handle all the algebra and the calculations of the explicit constants in the $`O()`$.777The program is several thousands of lines long. It is available from the first author for examination, further development etc. After careful selection of parameters, the program got (4) with an error term which is smaller than $`10^{25}|z|^{12}`$ for any $`|z|>1200`$. The program doing the algebra ran for half an hour on a PC. This may or may not be enough — it depends on the difference between the actual maximum and the asymptotic value. However, this difference is $`>0.01`$ (see below) and $`10^{25}|z|^{12}<10^5|z|^2`$ for $`|z|>1200`$ so this quantity is indeed negligible.

Interestingly, here the third approach (the one discussed briefly on page 1) shows a slight advantage over the other two. The third approach allows to get an explicit constant in the $`O()`$ in (8) with much less effort — a few trivial programs to find optimal parameters were all that was necessary. Contrariwise to the simplicity of the programs, the estimates are not as good and it was necessary to run the McCrea-Whipple algorithm for $`|z|<6000`$, which takes a few days.

The maximum of $`|z|^2\left(a(z)\frac{2}{\pi }\mathrm{log}|z|\lambda \right)`$ was found at $`z=3`$ and so the constant implicit in the $`O()`$ in (2) is

$$\frac{432+9\mathrm{log}72+18\gamma }{\pi }1530.06882$$

Notice that this is only slightly larger than the asymptotics $`\frac{1}{6\pi }0.05305`$. Naturally, explicit constants can be calculated for higher order estimates as well.

### 4.1. Other simulations

We add below two additional $`a`$’s. The first is for a “drunk king’s walk” — a random walk on $`^2`$ going at each step to each of the 8 neighboring lattice points with equal probability. Here are the first few coefficients:

$`a(z)`$ $`={\displaystyle \frac{4}{3\pi }}\mathrm{log}|z|+\lambda +\mathrm{Re}({\displaystyle \frac{z^4}{9\pi |z|^6}}+{\displaystyle \frac{11z^4}{90\pi |z|^8}}{\displaystyle \frac{5z^8}{36\pi |z|^{12}}}{\displaystyle \frac{167z^8}{252\pi |z|^{14}}}+`$
$`+{\displaystyle \frac{35z^{12}}{54\pi |z|^{18}}}{\displaystyle \frac{1673z^8}{2160\pi |z|^{16}}}+{\displaystyle \frac{15z^{12}}{2\pi |z|^{20}}}{\displaystyle \frac{1925z^{16}}{288\pi |z|^{24}}}+\mathrm{})`$
so in this case it is no longer true that all the coefficients have a constant sign. Another interesting case is the non-reversible walk on the directed triangular lattice going at each step with probability $`\frac{1}{3}`$ to $`z+1`$, $`z+\omega `$ or $`z+\omega ^2`$ where $`\omega =\frac{1}{2}+i\frac{\sqrt{3}}{2}=\sqrt[3]{1}`$. Here are the first few coefficients:

$`a(z)`$ $`={\displaystyle \frac{\sqrt{3}}{\pi }}\mathrm{log}|z|+\lambda +{\displaystyle \frac{\sqrt{3}}{\pi }}\mathrm{Re}({\displaystyle \frac{z^3}{6|z|^4}}+{\displaystyle \frac{z^6}{12|z|^8}}+{\displaystyle \frac{z^3}{18|z|^6}}+{\displaystyle \frac{5z^9}{54|z|^{12}}}+{\displaystyle \frac{17z^6}{135|z|^{10}}}+`$
$`{\displaystyle \frac{35z^{12}}{216|z|^{16}}}+{\displaystyle \frac{19z^9}{54|z|^{14}}}+{\displaystyle \frac{7z^{15}}{18|z|^{20}}}+{\displaystyle \frac{85z^6}{1134|z|^{12}}}+{\displaystyle \frac{98z^{12}}{81|z|^{18}}}+{\displaystyle \frac{385z^{18}}{324|z|^{24}}}+\mathrm{})`$
We see that the phenomenon of constant sign is not specific to the regular random walk (we checked this up to coefficients of order $`O(|z|^{28})`$). Do note that in this case all coefficients are positive while in the $`^2`$ case all coefficients are negative.