text:
# Lensing Properties of Scale-Free Galaxies

## 1 Introduction

There are now some sixty or so multiply imaged quasars known. Most of these exhibit double or quadruple images, although some triples and rings are also known (see Warren et al. 1999, Ratnatunga, Griffiths, Ostrander 1999, Wisotzki et al. 1999 and Tonry & Kochanek 1999 for some recent examples, as well as Pospieszalka et al. 1999 for details of the gravitational lensing database which maintains a list of candidates). This paper is mainly concerned with the sixteen or so quadruple systems. In propitious circumstances, the observables of such systems are the location of the lens, the positions of the four images, their parities and their fluxes. Lensing invariants are constructed directly from the lensing observables and remain invariant against changes in the source or lens configuration (provided a caustic is not crossed). A classical problem in the theory of gravitational lensing is the construction of a lens model that reproduces the observed properties of the images. One of the advantages of lensing invariants is that they provide simple tests by which we can determine quickly whether a given set of images can be produced by a particular lens.

Witt & Mao’s (1995) paper provides the first example of a lensing invariant. They considered a binary lens of two point masses, which has either three or five images according to whether the source lies outside or inside the caustic (e.g., Schneider, Ehlers & Falco 1992, section 8.3). Their main result is that the sum of the signed magnifications of the five images is unity when the source lies within the caustic, namely

$$\underset{i=1}{\overset{5}{}}\mu _ip_i=1.$$
(1)
Here, $`\mu _i`$ is the absolute magnification of the image, while $`p_i`$ is the parity. This result holds good irrespective of the position of the source, provided it remains within the central caustic. Subsequently, Rhie (1997) showed that the sum of the signed magnifications of the images of the $`N`$ point mass lens is also equal to unity, provided the source location yields the maximum number of $`N^2+1`$ images. These are remarkably simple results, bearing in mind that the lens equation that provides the locations of the images must be solved numerically. Further examples of lensing invariants were found by Dalal (1998) and Dalal & Rabin (2000).

Witt & Mao’s (2000) paper provides a study of lensing potentials stratified on similar concentric ellipses and falling off like a power-law with index $`\beta `$ suitable for modelling quadruple lenses. They showed that for the point mass ($`\beta =0`$) and isothermal ($`\beta =1`$) cases, there is an invariant

$$\underset{i=1}{\overset{4}{}}\mu _ip_i=\frac{2}{2\beta }=B.$$
(2)
They introduced higher-order moments and showed that they too had simple results in the point mass and isothermal cases

$$\underset{i=1}{\overset{4}{}}\mu _ip_ix_i=B\xi ,\underset{i=1}{\overset{4}{}}\mu _ip_iy_i=B\eta .$$
(3)
Here, $`(x_i,y_i`$) are the coordinates of the $`i`$th image, whereas $`(\xi ,\eta `$) is the position of the source. Witt & Mao (2000) argued on the basis of numerical experimentation that (2) and (3) remain good approximations for all other values of $`\beta `$. They suggested the use of interpolation to extend the result to all the power-law potentials. Witt & Mao’s (2000) paper is especially important to us here, as our results extend the work that they began.

In §2, we exploit contour integral methods to provide a new way of finding the invariants of a gravitational lens system. The idea is very simple. The images are the roots of the lens equation. So, by a judicious choice, the observables at each image can be made to equal the residue of a complex integrand. The lensing invariants, which are sums over the images, then become equal to a contour integral by Cauchy’s theorem. Although our method is more general, this paper concentrates on applications to the elliptic power-law potentials. In §3 we evaluate invariants for the case in which $`B=2/(2\beta )`$ is a positive integer. Our complex integrand is then single-valued. The simplest cases are the point mass ($`\beta =0`$, $`B=1`$) and the isothermal ($`\beta =1`$, $`B=2`$) lenses noted by Witt & Mao (2000). The invariants for $`B=3`$ include the contributions from a weak fifth image, while those for integer $`B>3`$ include small contributions from spurious images, which are non-physical solutions of the lens equation.

In §4, we progress to the case of non-integer $`B`$ when the contributions from a branch cut must be added to the earlier results. Remarkably this gives us exact expressions for the lensing invariants for the four bright images for all $`B`$, because the added contributions cancel out the effects of the fifth and the spurious images. Thus, our contour integral method provides not only a simpler but also a more powerful way to determine the lensing invariants than the earlier algebraic methods. Astrophysical applications of the invariants are presented in §5. They are useful because they can short-cut the modelling process, telling us quickly whether an assumed lensing potential (in this case the elliptic power-law potentials) can reproduce the data. We pay particular attention to the Einstein Cross gravitational lens system (G2237+0305) as our example. Finally, §6 discusses our conclusions and presents plans for future work.

## 2 The Contour Integral Method

### 2.1 The Lens Equation

At outset, we shall assume that the lensing potential is stratified on similar concentric ellipses

$`\psi =\{\begin{array}{cc}\frac{A}{\beta }(x^2+y^2q^2)^{\beta /2}\hfill & \text{if }0<\beta <2,\hfill \\ & \\ \frac{A}{2}\mathrm{log}(x^2+y^2q^2)\hfill & \text{ if }\beta =0.\hfill \end{array}`$ (7)
This family of models was introduced into gravitational lensing by Blandford & Kochanek (1987) and subsequently studied by others (e.g., Kassiola & Kovner 1993; Witt 1996; Witt & Mao 1997, 2000; Evans & Wilkinson 1998). Here, $`q`$ is the axis ratio of the projected equipotentials and is chosen to satisfy $`0<q1`$ without loss of generality. The parameter $`\beta `$ controls the radial profile of the potential while $`A`$ measures the depth of the well. The point mass (Schwarzschild lens) has $`\beta =0`$ and $`q=1`$, while the isothermal sphere has $`\beta =1`$ and $`q=1`$. The gravitational convergence (or projected density) is

$$\kappa =\frac{A}{2q^2}\frac{[1+q^2(\beta 1)]x^2+[1+q^2(\beta 1)]y^2}{(x^2+y^2q^2)^{2\beta /2}}.$$
(8)
This is positive definite provided $`\beta >1q^2`$. Three-dimensional analogues of these models are already familiar in galactic dynamics as the power-law models (Evans 1993, 1994). They are the only flattened and reasonably realistic galaxy models known, for which the self-gravitation equations can be solved to find simple two-integral distribution functions. For example, they have been used in modelling the nearby elliptical galaxy M32 (van der Marel et al. 1994), the inner parts of the Galactic bulge (Evans & de Zeeuw 1994), as well as the dark halo of our own Galaxy in the interpretation of the microlensing results (Alcock et al. 1997).

For the geometrically thin lens with projected potential (7), the paths of photons are given by (e.g., Schneider, Ehlers & Falco 1992, section 5.1)

$`\xi =x+\gamma _1x+\gamma _2y{\displaystyle \frac{Ax}{(x^2+y^2q^2)^{1\beta /2}}},\eta =y+\gamma _2x\gamma _1y{\displaystyle \frac{Ayq^2}{(x^2+y^2q^2)^{1\beta /2}}},`$ (9)
where ($`\xi ,\eta `$) are Cartesian coordinates of the source. Here, $`\gamma _1`$ and $`\gamma _2`$ allow for a constant external shear in an arbitrary direction. The external shear may be produced by the gravity field of a cluster or by the effects of nearby galaxies and usually lies between $`0<|\gamma |<0.3`$ (Keeton, Kochanek & Seljak 1997). We shall always assume that the lens is not circularly symmetric (that is, either $`q1`$ or $`\gamma _10`$ or $`\gamma _20`$). The lensing properties of circularly symmetric lenses are qualitatively different from those of more realistic flattened lenses. For example, the tangential caustic degenerates to a point in the case of circular symmetry (e.g., Schneider et al. 1992, section 8.1). So, circularly symmetric lenses can give misleading results and they are not widely used in modelling.

A formulation in terms of complex numbers has been used before to ease calculations in lensing theory (e.g., Bourassa, Kantowski & Norton 1973, Bourassa & Kantowski 1975, Witt 1990). However, we shall find it helpful to use somewhat different definitions, namely

$$\zeta =\xi +\mathrm{i}q\eta ,z=x+\mathrm{i}y/q.$$
(10)
The lens equation (9) becomes

$$\zeta =\frac{1}{2}[1+q^2+\gamma _1(1q^2)]z+\frac{1}{2}[1q^2+\gamma _1(1+q^2)+2\mathrm{i}q\gamma _2]\overline{z}\frac{Az}{(z\overline{z})^{1\beta /2}}.$$
(11)
It is convenient to introduce a variable $`t`$ such that:

$$t=\frac{A}{(z\overline{z})^{1/B}},B=\frac{2}{2\beta },$$
(12)
where the bar represents complex conjugation, and henceforth to set the constant $`A`$ equal to unity. The dependence of our results on $`A`$ can be recovered by dividing all powers of $`z`$ and $`\zeta `$ and their conjugates by the same powers of $`A^{B/2}`$. Then, the lens equation can be written in matrix form as

$$\left(\begin{array}{c}\zeta \hfill \\ \overline{\zeta }\hfill \end{array}\right)=\left(\begin{array}{cc}P\hfill & \hfill Q\\ \overline{Q}\hfill & \hfill P\end{array}\right)\left(\begin{array}{c}z\hfill \\ \overline{z}\hfill \end{array}\right),$$

with

$$P=\frac{1}{2}[1+q^2+\gamma _1(1q^2)]t=P_0t,Q=\frac{1}{2}[1q^2+\gamma _1(1+q^2)+2\mathrm{i}q\gamma _2],$$
(13)
from which it follows, by the standard algorithms of matrix inversion, that

$$\left(\begin{array}{c}z\hfill \\ \overline{z}\hfill \end{array}\right)=\frac{1}{P^2|Q|^2}\left(\begin{array}{cc}P\hfill & \hfill Q\\ \overline{Q}\hfill & \hfill P\end{array}\right)\left(\begin{array}{c}\zeta \hfill \\ \overline{\zeta }\hfill \end{array}\right).$$
(14)
By forming $`t^B=1/(z\overline{z})`$, we see that the solutions of the lens equation also satisfy the equation

$$K(t;\zeta ,\overline{\zeta })t^B(P\zeta Q\overline{\zeta })(P\overline{\zeta }\overline{Q}\zeta )\left[P^2|Q|^2\right]^2=0.$$
(15)
Let us designate $`K(t,\zeta ,\overline{\zeta })=0`$ as the imaging equation. Its real and positive solutions provide the possible image positions $`t`$ (or equivalently $`z`$), given the source location $`\zeta `$ and its conjugate $`\overline{\zeta }`$. The imaging equation must have at least as many solutions as the original lens equation (9). It may have more; that is, there may be solutions of (15) that do not correspond to true images and which we will call spurious roots.

For example, when $`B`$ is an integer $`N2`$, then the imaging equation is a quartic with four roots each of which corresponds to one of the images of a quadruple lens. When $`B=N>2`$, then the imaging equation is a polynomial of degree $`N+2`$. In the case $`B=3`$, the roots all correspond to the five images of a quintuple lens. However, when $`B>3`$, there are always roots of the imaging equation that are spurious, as there are never more than five images. This is shown in Appendix A.

### 2.2 The Magnification Invariant

Mathematically speaking, the lens equation defines a map from the lens plane to the source plane (see Schneider et al. 1992, chapter 5). The Jacobian matrix of this mapping is denoted by $`J`$, where

$$J=\left(\begin{array}{cc}\frac{\zeta }{z}|_{\overline{z}}\hfill & \hfill \frac{\zeta }{\overline{z}}|_z\\ \frac{\overline{\zeta }}{z}|_{\overline{z}}\hfill & \hfill \frac{\overline{\zeta }}{\overline{z}}|_z\end{array}\right)J^1=\left(\begin{array}{cc}\frac{z}{\zeta }|_{\overline{\zeta }}\hfill & \hfill \frac{z}{\overline{\zeta }}|_\zeta \\ \frac{\overline{z}}{\zeta }|_{\overline{\zeta }}\hfill & \hfill \frac{\overline{z}}{\overline{\zeta }}|_\zeta \end{array}\right).$$

With our choice of source and image variables, the reciprocal of the determinant of $`J`$ corresponds physically to the signed magnification of an image via

$$\frac{q^2}{\mathrm{det}J}|_{x_i,y_i}=\mu _ip_i,$$
(16)
where $`\mu _i`$ is the absolute value of the magnification and $`p_i`$ is the parity of the image located at $`(x_i,y_i)`$. For certain positions, $`\mathrm{det}J`$ may vanish and the magnification is infinite. These are the critical points and lines. The caustics are the images of the critical points and curves under the lens mapping (9).

By differentiating the expression for $`t`$, we obtain:

$`{\displaystyle \frac{B}{t^{1+B}}}{\displaystyle \frac{t}{\zeta }}|_{\overline{\zeta }}`$ $`=`$ $`{\displaystyle \frac{z}{\zeta }}|_{\overline{\zeta }}\overline{z}+{\displaystyle \frac{\overline{z}}{\zeta }}|_{\overline{\zeta }}z,`$ (17)
$`{\displaystyle \frac{B}{t^{1+B}}}{\displaystyle \frac{t}{\overline{\zeta }}}|_\zeta `$ $`=`$ $`{\displaystyle \frac{z}{\overline{\zeta }}}|_\zeta \overline{z}+{\displaystyle \frac{\overline{z}}{\overline{\zeta }}}|_\zeta z.`$ (18)
Recalling the properties of the Jacobian matrix $`J`$ and its inverse $`J^1`$, it follows that

$`{\displaystyle \frac{\zeta }{z}}|_{\overline{z}}{\displaystyle \frac{t}{\zeta }}|_{\overline{\zeta }}+{\displaystyle \frac{\overline{\zeta }}{z}}|_{\overline{z}}{\displaystyle \frac{t}{\overline{\zeta }}}|_\zeta `$ $`=`$ $`{\displaystyle \frac{t^{1+B}}{B}}\overline{z},`$ (19)
$`{\displaystyle \frac{\zeta }{\overline{z}}}|_z{\displaystyle \frac{t}{\zeta }}|_{\overline{\zeta }}+{\displaystyle \frac{\overline{\zeta }}{\overline{z}}}|_z{\displaystyle \frac{t}{\overline{\zeta }}}|_\zeta `$ $`=`$ $`{\displaystyle \frac{t^{1+B}}{B}}z.`$ (20)
By multiplying (19) by $`\overline{\zeta }_{\overline{z}}`$ and (20) by $`\overline{\zeta }_z`$ and subtracting, we construct the determinant of the Jacobian matrix and so establish

$$\mathrm{det}J\frac{t}{\zeta }|_{\overline{\zeta }}=\frac{t^{1+B}}{B}(\frac{\overline{\zeta }}{\overline{z}}|_z\overline{z}\frac{\overline{\zeta }}{z}|_{\overline{z}}z).$$
(21)
However, equation (15) provides us with the requirement that $`K(t;\zeta ,\overline{\zeta })=0`$ and so, using a standard formula in the theory of partial differentiation, we have

$$\frac{1}{\mathrm{det}J}=\frac{B}{t^{1+B}(\frac{\overline{\zeta }}{\overline{z}}|_z\overline{z}\frac{\overline{\zeta }}{z}|_{\overline{z}}z)}\frac{\frac{K}{\zeta }|_{\overline{\zeta },t}}{\frac{K}{t}|_{\zeta ,\overline{\zeta }}}.$$
(22)
Now comes a remarkable simplification. For the ellipsoidally stratified potentials, we find by direct evaluation that

$$\frac{\overline{\zeta }}{\overline{z}}|_z\overline{z}\frac{\overline{\zeta }}{z}|_{\overline{z}}z=\frac{(P^2+|Q|^2)\overline{\zeta }2P\overline{Q}\zeta }{P^2|Q|^2}.$$
(23)
A similar factor occurs on partially differentiating $`K(t;\zeta ,\overline{\zeta })`$; in other words

$$\frac{K}{\zeta }|_{\overline{\zeta },t}=t^B[(P^2+|Q|^2)\overline{\zeta }2P\overline{Q}\zeta ].$$
(24)
So, equation (22) simplifies to

$$\frac{1}{\mathrm{det}J}=\frac{B}{t}\frac{P^2|Q|^2}{\frac{K}{t}|_{\zeta ,\overline{\zeta }}}.$$
(25)
Let us recall that the images correspond to simple zeros of $`K(t;\zeta ,\overline{\zeta })`$ and the residue of $`1/K`$ at a simple zero is just $`1/K_t`$. This gives us a contour integral representation of the magnification invariant, which is defined as the sum of the signed magnifications of the images

$$\underset{\mathrm{images}}{}\mu _ip_i=\frac{q^2B}{2\pi \mathrm{i}}_𝒞\frac{dt}{t}\frac{P^2|Q|^2}{K(t;\zeta ,\overline{\zeta })},$$
(26)
where $`𝒞`$ is a contour in the complex lens plane enclosing only the simple poles corresponding to visible images. It is an invariant, as it is a pure number which depends only on the structure of the lensing potential or, equivalently, the mass distribution.

### 2.3 The Configuration Invariants

The higher order invariants are sums of powers of the coordinates of the images weighted with their signed magnifications. We shall call them the configuration invariants. From the lens equation (9), we deduce that

$`x(t;\zeta ,\overline{\zeta })`$ $`=`$ $`{\displaystyle \frac{q^2[(1\gamma _1tq^2)\xi \gamma _2\eta ]}{P^2|Q|^2}},`$ (27)
$`y(t;\zeta ,\overline{\zeta })`$ $`=`$ $`{\displaystyle \frac{q^2[(1+\gamma _1t)\eta \gamma _2\xi ]}{P^2|Q|^2}}.`$ (28)
The higher order configuration invariants are the sums over the images of products of the signed magnifications with the position coordinates. The $`m`$th moment with respect to $`x`$ and $`n`$th moment with respect to $`y`$ is given by the contour integral

$$\underset{\mathrm{images}}{}\mu _ip_ix_i^my_i^n=\frac{q^2B}{2\pi \mathrm{i}}_𝒞\frac{dt}{t}\frac{P^2|Q|^2}{K(t;\zeta ,\overline{\zeta })}[x(t;\zeta ,\overline{\zeta })]^m[y(t;\zeta ,\overline{\zeta })]^n.$$
(29)
Witt & Mao (2000) introduced some higher order moments of $`x`$ and $`y`$ for the cases when the integers $`m`$ and $`n`$ are positive and they showed that the results for the lensing potentials (7) were very simple for the point mass ($`\beta =0`$) and isothermal ($`\beta =1`$) cases.

It is also possible to derive moments with respect to the complex image coordinates. The $`m`$th moment with respect to $`z`$ and the $`n`$th moment with respect to $`\overline{z}`$ is given by the contour integral

$$\underset{\mathrm{images}}{}\mu _ip_iz_i^m\overline{z}_i^n=\frac{q^2B}{2\pi \mathrm{i}}_𝒞\frac{dt}{t}\frac{[P\zeta Q\overline{\zeta }]^m[P\overline{\zeta }\overline{Q}\zeta ]^n}{K(t;\zeta ,\overline{\zeta })[P^2|Q|^2]^{m+n1}}$$
(30)
For positive moments, the results are of course just a recasting of (29), but for reciprocal moments, the results are different. In both cases (29) and (30), the contour $`𝒞`$ is chosen to enclose the poles corresponding to the visible images, but to exclude other poles in the complex $`t`$-plane such as those at the zeros of $`P^2|Q|^2`$ when $`m+n>1`$.

## 3 Invariants for integer $`B`$

To begin with, let us study elliptic power-law potentials for which the exponent $`B=2/(2\beta )`$ is an integer. This may strike the reader as narrow in its scope, but it is not really so. As we will see in §4, results which are exactly true for these cases are also good approximations for the four bright images of all the power-law galaxies if the source is sufficiently close to the galactic center. The simplest cases are $`B=1`$, or $`\beta =0`$, when the potential approximates that of a point mass plus shear, and $`B=2`$, or $`\beta =1`$, when the potential is isothermal-like and represents a galaxy or cluster with a flat rotation curve. In both these cases, the imaging equation (15) is a quartic and has four roots. When $`B=3`$ or $`\beta =4/3`$, so that the convergence or projected mass density falls off like $`r^{2/3}`$, the imaging equation (15) is a quintic and has five roots and there is a weak fifth image. The cases $`B=1,2`$ and $`3`$ are special because there are then no spurious roots of the imaging equation when the source lies within the central caustic. There are spurious roots of the imaging equation for $`B>3`$, which contribute to the invariants which we calculate. However their contributions are small and many of them cancel, with the consequence that the results remain useful even when there are spurious roots. Cases of integer $`B`$ avoid the need for a branch cut to define the $`t^B`$ term in the imaging equation (15) as a single-valued function in the complex plane.

### 3.1 The Magnification Invariant

As the images are given by the zeros of the imaging equation, we have

$$\underset{\mathrm{images}}{}\mu _ip_i=\frac{q^2B}{2\pi \mathrm{i}}_𝒞\frac{dt}{t}\frac{P^2|Q|^2}{K(t;\zeta ,\overline{\zeta })}.$$
(31)
The integrand has simple poles at the zeros of $`K(t;\zeta ,\overline{\zeta })`$, as well as at the origin. The contour $`𝒞`$ is traversed in an anti-clockwise direction and encloses the poles corresponding to images.

Suppose the integrand is taken over a large circle at infinity $`𝒞_{\mathrm{}}`$. As $`t\mathrm{}`$, the integrand falls off at least as fast as $`O(t^3)`$, so we have;

$$_𝒞_{\mathrm{}}\frac{dt}{t}\frac{P^2|Q|^2}{K(t;\zeta ,\overline{\zeta })}=0.$$
(32)
Distorting the contour so it encloses the zeros of $`K(t;\zeta ,\overline{\zeta })`$ and the pole at the origin, we have

$$\underset{\mathrm{images}}{}\mu _ip_i+\frac{q^2B}{2\pi \mathrm{i}}_{𝒞_ϵ}\frac{dt}{t}\frac{P^2|Q|^2}{K(t;\zeta ,\overline{\zeta })}=0,$$
(33)
where $`𝒞_ϵ`$ is a small circle of radius $`ϵ`$ enclosing the origin in the anticlockwise direction. As $`t0`$, then $`K(t;\zeta ,\overline{\zeta })(P_0^2|Q|^2)^2,`$ where $`P_0=P(t=0)=\frac{1}{2}[1+q^2+\gamma _1(1q^2)]`$ so that

$$\underset{\mathrm{images}}{}\mu _ip_i=\frac{Bq^2}{P_0^2|Q|^2}=\frac{B}{1\gamma _1^2\gamma _2^2}.$$
(34)
This remarkably simple result was first obtained by Witt & Mao (2000), who deduced it for the $`B=1`$ and $`2`$ cases and used interpolation to extend it as an approximation to all $`B`$. In fact, the result (34) is exact for $`B=3`$ as well, provided the contribution from the fifth de-magnified image is included.

### 3.2 The Configuration Invariants: First and Second Moments

The first moments can be obtained in a similar way. When $`m=1,n=0`$ the integrand (30) still has poles only at the zeros of the imaging equation and the origin. Evaluating the residue at the origin, we find

$$\underset{\mathrm{images}}{}\mu _ip_iz_i=\frac{B[P_0\zeta Q\overline{\zeta }]}{q^2(1\gamma _1^2\gamma _2^2)^2}$$
(35)
and hence

$$\underset{\mathrm{images}}{}\mu _ip_ix_i=\frac{B((1\gamma _1)\xi \gamma _2\eta )}{(1\gamma _1^2\gamma _2^2)^2},\underset{\mathrm{images}}{}\mu _ip_iy_i=\frac{B((1+\gamma _1)\eta \gamma _2\xi )}{(1\gamma _1^2\gamma _2^2)^2}.$$
(36)
Witt & Mao (2000) give these results in the case of on-axis shear only ($`\gamma _2=0`$), while Dalal & Rabin (2000) give them for the case $`B=2`$ only. Algorithms using the real analysis tend to become cumbersome with off-axis shear. It is part of the power of the contour integral method that results can be obtained easily for general shear. When there is no off-axis shear, the $`x`$ (or $`y`$) moments of the signed magnifications of the images are related only to the $`\xi `$ (or $`\eta `$) offset of the lens.

The second moments are given by $`m+n=2`$. The integrand (30) now has additional simple poles at the two real points $`t_1`$ and $`t_2`$, where

$$t_1=P_0+|Q|,t_2=P_0|Q|,$$
(37)
where $`P^2|Q|^2=0`$. Provided that $`\gamma _1^2+\gamma _2^2<1`$, then $`P_0>0`$ and $`P_0^2|Q|^2=q^2(1\gamma _1^2\gamma _2^2)>0`$, and both $`t_1`$ and $`t_2`$ are positive with $`t_1>t_2`$ and $`t_2<1`$. The residues at these locations give rise to constant terms independent of the source position in the configuration invariants. To show this, let us consider the evaluation of one of the second moments ($`m=1,n=1`$) in a little detail. At large radii, the integrand falls off at least as fast as $`O(t^5)`$. If evaluated about a circle at infinity $`𝒞_{\mathrm{}}`$, the contour integral (30) vanishes. Distorting the contour so it encloses the zeros of $`K(t;\zeta ,\overline{\zeta })`$ and the poles at $`t=0,t_1`$ and $`t_2`$, we have

$$\underset{\mathrm{images}}{}\mu _ip_iz_i\overline{z}_i+\frac{q^2B}{2\pi \mathrm{i}}_{𝒞_ϵ+𝒞_{t_1}+𝒞_{t_2}}\frac{dt}{t}\frac{[P\zeta Q\overline{\zeta }][P\overline{\zeta }\overline{Q}\zeta ]}{[P^2|Q|^2]K(t;\zeta ,\overline{\zeta })}=0.$$
(38)
Here, $`𝒞_ϵ`$, $`𝒞_{t_1}`$ and $`𝒞_{t_2}`$ are small circles of radius $`ϵ`$ enclosing the simple poles at $`t=0,t_1`$ and $`t_2`$ respectively. By evaluating the residues, we find;

$$\underset{\mathrm{images}}{}\mu _ip_iz_i\overline{z}_i=\frac{q^2B}{2|Q|}\left\{\frac{1}{t_2^{B+1}}\frac{1}{t_1^{B+1}}\right\}+\frac{B(P_0\zeta Q\overline{\zeta })(P_0\overline{\zeta }\overline{Q}\zeta )}{q^4(1\gamma _1^2\gamma _2^2)^3}$$
(39)
After evaluating the moment of $`z_i^2`$ in a similar manner, we obtain the real second moments:

$`{\displaystyle \underset{\mathrm{images}}{}}\mu _ip_ix_i^2`$ $`=`$ $`{\displaystyle \frac{q^2B}{4|Q|^2}}\left\{{\displaystyle \frac{|Q|\mathrm{Re}(Q)}{t_2^{B+1}}}{\displaystyle \frac{|Q|+\mathrm{Re}(Q)}{t_1^{B+1}}}\right\}+{\displaystyle \frac{B((1\gamma _1)\xi \gamma _2\eta )^2}{(1\gamma _1^2\gamma _2^2)^3}},`$ (40)
$`{\displaystyle \underset{\mathrm{images}}{}}\mu _ip_iy_i^2`$ $`=`$ $`{\displaystyle \frac{q^4B}{4|Q|^2}}\left\{{\displaystyle \frac{|Q|+\mathrm{Re}(Q)}{t_2^{B+1}}}{\displaystyle \frac{|Q|\mathrm{Re}(Q)}{t_1^{B+1}}}\right\}+{\displaystyle \frac{B((1+\gamma _1)\eta \gamma _2\xi )^2}{(1\gamma _1^2\gamma _2^2)^3}},`$ (41)
$`{\displaystyle \underset{\mathrm{images}}{}}\mu _ip_ix_iy_i`$ $`=`$ $`{\displaystyle \frac{\gamma _2q^4B}{4|Q|}}\left\{{\displaystyle \frac{1}{t_2^{B+1}}}+{\displaystyle \frac{1}{t_1^{B+1}}}\right\}+{\displaystyle \frac{B((1+\gamma _1)\eta \gamma _2\xi )((1\gamma _1)\xi \gamma _2\eta )}{(1\gamma _1^2\gamma _2^2)^3}}.`$ (42)
These results simplify considerably when the shear acts only on axis $`(\gamma _2=0)`$, in which case $`|Q|=\mathrm{Re}(Q)`$, $`t_1=1+\gamma _1`$, and $`t_2=q^2(1\gamma _1)`$. The $`x^2`$ (or $`xy`$ or $`y^2`$) moments of the signed magnifications of the images are related only to the $`\xi ^2`$ (or $`\xi \eta `$ or $`\eta ^2`$) offset of the lens. The results (40) - (42) are given by Witt & Mao (2000) for the case of vanishing off-axis shear. For completeness, the third moments are derived in Appendix B.

### 3.3 The Configuration Invariants: Reciprocal Moments

The reciprocal moments can also be found. The first order reciprocal moments are given by $`m=1,n=0`$ or $`m=0,n=1`$. In addition to the poles at the images and the origin, the integrand also has simple poles at either $`t_3`$ or $`t_4`$ given by

$$t_3=q^2\left(1\gamma _1\frac{\gamma _2\eta }{\xi }\right),t_4=1+\gamma _1\frac{\gamma _2\xi }{\eta }.$$
(43)
The real first order reciprocal moments are

$$\underset{\mathrm{images}}{}\frac{\mu _ip_i}{x_i}=C_1+\frac{B}{(1\gamma _1)\xi \gamma _2\eta },\underset{\mathrm{images}}{}\frac{\mu _ip_i}{y_i}=C_2+\frac{B}{(1+\gamma _1)\eta \gamma _2\xi }.$$
(44)
The constants $`C_1`$ and $`C_2`$ are

$`C_1`$ $`=`$ $`{\displaystyle \frac{q^2B\gamma _2^2}{\xi t_4[\xi ^2t_4^B\gamma _2^2]}}={\displaystyle \frac{B\gamma _2^2}{[(1\gamma _1)\xi \gamma _2\eta ][\gamma _2^2+\xi ^{2B}q^{2B}((1\gamma _1)\xi \gamma _2\eta )^B]}},`$ (45)
$`C_2`$ $`=`$ $`{\displaystyle \frac{B\gamma _2^2}{\eta t_5[\eta ^2t_5^B\gamma _2^2]}}={\displaystyle \frac{B\gamma _2^2}{[(1+\gamma _1)\eta \gamma _2\xi ][\gamma _2^2+\eta ^{2B}((1+\gamma _1)\eta \gamma _2\xi )^B]}}.`$ (46)
Sending $`\xi \eta `$, $`\eta \xi `$, $`\gamma _1\gamma _1`$, $`q1`$, then we see that $`C_1C_2`$. This follows as $`C_1`$ is the residue of the integrand at $`t_3`$, and $`C_2`$ at $`t_4`$. Note that $`C_1`$ and $`C_2`$ both vanish when the shear is on-axis only ($`\gamma _2=0`$). In this instance, it is again true that the reciprocal $`x`$ (or $`y`$) moments of the signed magnifications depend only on the reciprocal $`\xi `$ (or $`\eta `$) offset of the source.

The complex reciprocal moments have poles of order $`m`$ or $`n`$ at $`t=t_5=P_0Q\overline{\zeta }/\zeta `$ and/or its conjugate $`t=\overline{t}_5=P_0\overline{Q}\zeta /\overline{\zeta }`$. For $`m=1`$ and $`n=0`$, the residues at $`t=0`$ and $`t=t_5`$ exactly cancel so that

$$\underset{\mathrm{images}}{}\frac{\mu _ip_i}{z_i}=0.$$
(47)
This implies the vanishing of the following two real moments

$$\underset{\mathrm{images}}{}\frac{\mu _ip_ix_i}{x_i^2+y_i^2q^2}=\underset{\mathrm{images}}{}\frac{\mu _ip_iy_i}{x_i^2+y_i^2q^2}=0.$$
(48)
A complication that arises with sufficiently negative $`m+n`$ is that there is an additional contribution from the integration over the large circle at infinity $`C_{\mathrm{}}`$ to be included. This happens when $`(m+n)`$ equals or exceeds the larger of $`B`$ and $`2`$, and hence with second order and higher reciprocal moments. The second order reciprocal moments in complex form have the compact expressions

$$\underset{\mathrm{images}}{}\frac{\mu _ip_i}{z_i\overline{z}_i}=\underset{\mathrm{images}}{}\frac{\mu _ip_i}{(x_i^2+y_i^2q^2)^2}=\{\begin{array}{cc}0,\hfill & B<2,\hfill \\ \multicolumn{2}{c}{}\\ \frac{2q^2}{|\zeta |^21},\hfill & B=2,\hfill \\ \multicolumn{2}{c}{}\\ \frac{Bq^2}{|\zeta |^2},\hfill & B>2,\hfill \end{array}$$
(49)
and

$$\underset{\mathrm{images}}{}\frac{\mu _ip_i}{z_i^2}=\underset{\mathrm{images}}{}\frac{\mu _ip_i}{(x_i^2+y_i^2q^2)^2}\left[x_i^2\frac{y_i^2}{q^2}\frac{2\mathrm{i}x_iy_i}{q}\right]=\{\begin{array}{cc}\frac{Bq^2\overline{Q}t_5^{B1}}{|Q|^2},\hfill & B<2,\hfill \\ \multicolumn{2}{c}{}\\ 2q^2\left[\frac{\overline{\zeta }^2}{|\zeta |^21}\frac{\overline{Q}t_5}{|Q|^2}\right],\hfill & B=2,\hfill \\ \multicolumn{2}{c}{}\\ Bq^2\left[\frac{1}{\zeta ^2}\frac{\overline{Q}t_5^{B1}}{|Q|^2}\right],\hfill & B>2.\hfill \end{array}$$
(50)
There are also moments for which one of $`m`$ and $`n`$ is negative, and the other positive. The simplest is that for $`m=1`$ and $`n=1`$, which gives

$$\underset{\mathrm{images}}{}\frac{\mu _ip_iz_i}{\overline{z}_i}=\frac{BP_0}{\overline{Q}(1\gamma _1^2\gamma _2^2)},$$
(51)
or, in real form,

$`{\displaystyle \underset{\mathrm{images}}{}}{\displaystyle \frac{\mu _ip_i(x_i^2y_i^2q^2)}{x_i^2+y_i^2q^2}}`$ $`=`$ $`{\displaystyle \frac{BP_0\mathrm{Re}(Q)}{(1\gamma _1^2\gamma _2^2)|Q|^2}},`$
$`{\displaystyle \underset{\mathrm{images}}{}}{\displaystyle \frac{\mu _ip_ix_iy_i}{x_i^2+y_i^2q^2}}`$ $`=`$ $`{\displaystyle \frac{\gamma _2q^2BP_0}{2(1\gamma _1^2\gamma _2^2)|Q|^2}}.`$ (52)
### 3.4 Spurious Roots

Now let us confront the problem of the spurious roots of the imaging equation (15) that do not correspond to physical images. When $`B2`$, there are four roots on the real positive $`t`$-axis corresponding to the four images of a quadruple lens. When $`B>2`$, there are $`B+2`$ roots in total. With $`|\zeta |`$ small, the $`B2`$ additional roots appear at large $`t`$ (or at small $`|z|`$). By balancing terms in the imaging equation, we see that these roots $`t_\mathrm{s}`$ satisfy to leading order

$$t_\mathrm{s}=\frac{1}{|\zeta |^{2/(B2)}}\mathrm{exp}\left(\frac{2n\pi \mathrm{i}}{B2}\right)+O(1),$$
(53)
where $`n`$ is an integer so that there are $`B2`$ distinct roots of this form. Only one of them ($`n=0`$) lies on the real positive $`t`$-axis and corresponds to a physical image. In fact, it is the strongly de-magnified fifth image that occurs so close to the centre of the lens that it is often hard to detect. The other $`B3`$ roots ($`0<nB3`$) are spurious and do not provide physical images. Figure 1 illustrates the behaviour of the roots in the complex $`t`$-plane as $`B`$ increases. The first panel shows the case $`B=3`$, with the five roots corresponding to the five images along the real positive $`t`$-axis. The second panel shows the case $`B=4`$, when a spurious root appears at large negative real $`t`$. When $`4<B<6`$, this negative real root splits into a pair with arguments $`\mathrm{exp}\pm 2\pi \mathrm{i}/(B2)`$. In the third panel, at $`B=5`$, we have three roots equally spaced around a circle of large radius in the complex $`t`$-plane. In the fourth panel, the cycle of events begins again with a new spurious root appearing on the negative real axis when $`B=6`$.

Some small part of the invariants that we calculated earlier in this section are due to spurious roots. Their contributions can now be calculated and subtracted. We shall do this for the sum of the signal magnitudes. In fact, we shall subtract also the contribution from the weak fifth image to obtain an approximation for the sum of the signed magnitude of the four bright images. The contribution to the integral (31) from the residue at each root (53) is

$$\mathrm{Residue}(t=t_\mathrm{s})=\frac{Bq^2}{(B2)}|\zeta |^{4/(B2)}\mathrm{exp}\left(\frac{4n\pi \mathrm{i}}{B2}\right)+O(|\zeta |^{6/(B2)}).$$
(54)
For odd values of $`B`$, it is convenient to set $`B=3+2N`$. The spurious roots lie at angles $`\mathrm{exp}(\pm 2\pi \mathrm{i}n/(2N+1))`$ and so

$`{\displaystyle \underset{4\mathrm{images}}{}}\mu _ip_i`$ $`=`$ $`{\displaystyle \frac{B}{(1\gamma _1^2\gamma _2^2)}}{\displaystyle \frac{Bq^2}{(B2)}}|\zeta |^{4/(B2)}`$ (55)
$`\times \left[1+{\displaystyle \underset{n=1}{\overset{N}{}}}\mathrm{exp}\left({\displaystyle \frac{4n\pi \mathrm{i}}{2N+1}}\right)+\mathrm{exp}\left({\displaystyle \frac{4n\pi \mathrm{i}}{2N+1}}\right)\right]+O(|\zeta |^{6/(B2)}).`$
For even integers $`B=2+2N`$, the spurious roots lie at angles $`\mathrm{exp}(\pm n\pi i/N)`$ and so

$`{\displaystyle \underset{4\mathrm{images}}{}}\mu _ip_i`$ $`=`$ $`{\displaystyle \frac{B}{(1\gamma _1^2\gamma _2^2)}}{\displaystyle \frac{Bq^2}{(B2)}}|\zeta |^{4/(B2)}`$ (56)
$`\times \left[2+{\displaystyle \underset{n=1}{\overset{N1}{}}}\mathrm{exp}\left({\displaystyle \frac{4n\pi \mathrm{i}}{N}}\right)+\mathrm{exp}\left({\displaystyle \frac{4n\pi \mathrm{i}}{N}}\right)\right]+O(|\zeta |^{6/(B2)}).`$
Except for the cases $`B=3`$ and $`B=4`$, the terms in square brackets on the right-hand side of (55) and (56) sum to zero, and so there is no contribution to the magnification invariant at this order. This is not the case for $`B=3`$, when the sum in (55) is empty, and we obtain:

$$\underset{4\mathrm{images}}{}\mu _ip_i=\frac{3}{(1\gamma _1^2\gamma _2^2)}3q^2|\zeta |^4+O(|\zeta |^6),$$
(57)
or for $`B=4`$ when the sum in (56) is empty, and we obtain

$$\underset{4\mathrm{images}}{}\mu _ip_i=\frac{4}{(1\gamma _1^2\gamma _2^2)}4q^2|\zeta |^2+O(|\zeta |^4).$$
(58)
In both cases, the corrections match those given by the more general theory which we develop in the next section.

## 4 Invariants for all $`B`$

The general case in which $`B`$ is no longer an integer requires that a branch cut be introduced to make single-valued the $`t^B`$ term in the imaging equation (11). This apparent complication is a blessing in disguise. It allows us to generate invariants which are valid for a continuous range of values of $`B`$. They are infinite series in the source coordinates, whose leading terms are just those that we derived in §3 — namely, eqns (34) - (36), (39) - (42), and (44). Importantly, they are also analytic functions of $`B`$, and hence are useful for all values of $`B`$ in the physically relevant range $`1B<\mathrm{}`$ ($`0\beta <2`$). This section explains why this is so. We first derive the magnification invariant for the case of general $`B<2`$ (§4.1) and then the configuration invariants (§4.2). In §4.3 we show how these invariants could also be derived more laboriously by perturbation theory. Readers happy to take the results on trust can turn directly to the applications (§5).

### 4.1 A Recalculated Magnification Invariant

Suppose that $`B`$ is not an integer and that $`1B2`$, so that there are no spurious roots and only the four bright images. We make $`t^B`$ single-valued and real on the positive real $`t`$-axis with a cut along the negative real axis. We again use the contour integral (31) for the magnification invariant. We enlarge the contour $`C`$ as much as possible — the integrand has no singularities other than those on the cut — until it consists of the infinite loop $`𝒞_{\mathrm{cut}}`$ shown in Figure 2 together with a large circle at infinity. The large circle again contributes nothing, because of the decay of the integrand for large $`t`$. Hence, the contour integral becomes

$$\underset{\mathrm{images}}{}\mu _ip_i=\frac{q^2B}{2\pi \mathrm{i}}_{𝒞_{\mathrm{cut}}}\frac{dt}{t}\frac{P^2|Q|^2}{K(t;\zeta ,\overline{\zeta })}.$$
(59)
Let us write $`K=K_1+K_2`$, where

$$K_1=(P^2|Q|^2)^2=(tt_1)^2(tt_2)^2,K_2=t^B(P\zeta Q\overline{\zeta })(P\overline{\zeta }\overline{Q}\zeta ).$$
(60)
We introduce the following two linear combinations of the complex source coordinates:

$$\lambda =\frac{1}{2}\left(\zeta +\frac{Q\overline{\zeta }}{|Q|}\right),\nu =\frac{1}{2}\left(\zeta \frac{Q\overline{\zeta }}{|Q|}\right).$$
(61)
They have the important properties that

$$\overline{\lambda }\nu +\lambda \overline{\nu }=0,|\lambda |^2+|\nu |^2=|\zeta |^2,$$
(62)
and allow us to write

$$P\zeta Q\overline{\zeta }=\lambda (t_2t)+\nu (t_1t),(P\zeta Q\overline{\zeta })(P\overline{\zeta }\overline{Q}\zeta )=|\lambda |^2(t_2t)^2+|\nu |^2(t_1t)^2.$$
(63)
We expand the denominator of our integrand in an infinite series to get

$$_{𝒞_{\mathrm{cut}}}\frac{dt}{t}\frac{P^2|Q|^2}{K(t;\zeta ,\overline{\zeta })}=\underset{j=0}{\overset{\mathrm{}}{}}_{𝒞_{\mathrm{cut}}}𝑑t\frac{P^2|Q|^2}{tK_1}\left(\frac{K_2}{K_1}\right)^j.$$
(64)
This is valid for sufficiently small $`|\zeta |`$ because both $`|K_1|`$ and $`|K_2|`$ are bounded below on the cut for finite $`t`$ and $`|K_2/K_1|0`$ as $`|t|\mathrm{}`$. The cut is unnecessary for the first $`j=0`$ term of the series (64). It has just a simple pole at $`t=0`$ and so can easily be evaluated. The result is exactly the same as that of §3.1 as the residue calculation is unchanged. In the other terms of the integrand, we set $`t=w\mathrm{exp}(\pm \pi \mathrm{i})`$ on the upper and lower sides of the cut to get

$$\underset{4\mathrm{images}}{}\mu _ip_i=\frac{B}{1\gamma _1^2\gamma _2^2}+Bq^2\underset{j=1}{\overset{\mathrm{}}{}}\frac{\mathrm{sin}(Bj\pi )}{\pi }_0^{\mathrm{}}𝑑ww^{jB1}\frac{[|\lambda |^2(w+t_2)^2+|\nu |^2(w+t_1)^2]^j}{(w+t_1)^{2j+1}(w+t_2)^{2j+1}}.$$
(65)
Binomial expansion of the integrand gives us integrals which are all of the form

$$I_{j,\mathrm{}}(t_1,t_2)=\frac{\mathrm{sin}Bj\pi }{\pi }_0^{\mathrm{}}\frac{w^{jB1}dw}{(w+t_1)^{\mathrm{}+1}(w+t_2)^{2j\mathrm{}+1}},1B2,$$
(66)
for integers $`\mathrm{}`$ in $`[0,2j]`$. The integrals can all be expressed in terms of Gaussian hypergeometric functions $`{}_{2}{}^{}F_{1}^{}`$, using formula \[3.197.1\] of Gradshteyn & Ryzhik (1965) as

$$I_{j,\mathrm{}}=\frac{1}{(2j+1)!}\underset{s=1}{\overset{2j+1}{}}(jBs)\frac{t_2^{jB+\mathrm{}2j1}}{t_1^{\mathrm{}+1}}{}_{2}{}^{}F_{1}^{}(\mathrm{}+1,jB;2j+2;1\frac{t_2}{t_1}).$$
(67)
Appendix C gives closed formulae for some of the low order integrals $`I_{j,\mathrm{}}`$. Binomial expansion of the numerator of the integrand gives the following infinite series for the magnification invariant of the four bright images. We henceforth label this as $`(B)`$ to emphasise that it is an analytic function of $`B`$;

$$(B)=\underset{4\mathrm{images}}{}\mu _ip_i=\frac{B}{1\gamma _1^2\gamma _2^2}+Bq^2\underset{j=1}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{j}{}}\left(\genfrac{}{}{0pt}{}{j}{m}\right)I_{j,2m}|\lambda |^{2m}|\nu |^{2j2m}$$
(68)
The $`j`$th term is a homogeneous polynomial of degree $`2j`$ in the source coordinates. Writing out the first correction term explicitly, we have

$$(B)=\frac{B}{1\gamma _1^2\gamma _2^2}+Bq^2\left[K_{\mathrm{2\hspace{0.17em}0}}^{\mathrm{0\hspace{0.17em}0}}\xi ^2+K_{\mathrm{1\hspace{0.17em}1}}^{\mathrm{0\hspace{0.17em}0}}\xi \eta +K_{\mathrm{0\hspace{0.17em}2}}^{\mathrm{0\hspace{0.17em}0}}\eta ^2\right]+O(|\zeta |^4),$$
(69)
where

$`K_{\mathrm{2\hspace{0.17em}0}}^{\mathrm{0\hspace{0.17em}0}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}[I_{1,2}+I_{1,0}]+{\displaystyle \frac{\mathrm{Re}Q}{2|Q|}}[I_{1,2}I_{1,0}],K_{\mathrm{1\hspace{0.17em}1}}^{\mathrm{0\hspace{0.17em}0}}={\displaystyle \frac{\gamma _2q^2}{|Q|}}[I_{1,2}I_{1,0}],`$
$`K_{\mathrm{0\hspace{0.17em}2}}^{\mathrm{0\hspace{0.17em}0}}`$ $`=`$ $`{\displaystyle \frac{q^2}{2}}[I_{1,2}+I_{1,0}]{\displaystyle \frac{\mathrm{Re}Q}{2|Q|}}[I_{1,2}I_{1,0}],`$ (70)
and, using eqn (C), we find,

$`I_{1,2}I_{1,0}`$ $`=`$ $`{\displaystyle \frac{(B1)}{4|Q|^2}}\left[t_1^{B2}t_2^{B2}(B2)|Q|(t_1^{B3}+t_2^{B3})\right],`$
$`I_{1,2}+I_{1,0}`$ $`=`$ $`{\displaystyle \frac{(B1)}{4|Q|^2}}\left[t_1^{B2}+t_2^{B2}(B2)|Q|(t_1^{B3}t_2^{B3})\right]+{\displaystyle \frac{1}{4|Q|^3}}\left[t_2^{B1}t_1^{B1}\right].`$ (71)
This gives explicit formulae for the first order corrections to the magnification invariant. Once again, the formulae simplify considerably if there is no off-axis shear and $`Q`$ is real.

Although the preceding derivation of the infinite series (68) for the magnification invariant is valid only for $`1B2`$, the series itself is not so limited. As we show in Appendix D, it converges for arbitrarily large values of $`B`$ for sufficiently small $`|\zeta |`$. In fact, the series is the unique analytic continuation of the magnification invariant $`(B)`$, which is an analytic function of $`B`$, and hence can be used for $`B>2`$. Here are two simple tests that give credence to this remarkable result, the justification of which is given in Appendix D. First, for $`B=3`$, $`I_{1,2}=I_{1,0}=0`$ and there is no $`O(|\zeta |^2)`$ term in the series (68). At the next order, we find from equation (67) that $`I_{2,4}=I_{2,2}=I_{2,0}=1`$ and hence

$$(3)=\frac{3}{1\gamma _1^2\gamma _2^2}3q^2|\zeta |^4+O(|\zeta |^6),$$
(72)
in agreement with the result (57) of §3.4. Second, let us set $`B=4`$. Then $`I_{1,2}=I_{1,0}=1`$, and evaluating the $`j=1`$ term of (68) and using (62) gives

$$(4)=\frac{4}{1\gamma _1^2\gamma _2^2}4q^2|\zeta |^2+O(|\zeta |^4).$$
(73)
This matches the approximate result (58) we obtained in §3.4 after subtracting the weak fifth and one spurious image.

### 4.2 Recalculated Configuration Invariants

We now obtain infinite series expansions for the configuration invariants (30), also summed over the four bright images, by manipulating that contour integral in a manner similar to that of the previous section. The integrand now generally has additional finite singularities at one or more of $`t_1`$, $`t_2`$, $`t_5`$, and $`\overline{t}_5`$, and so their residues must be accounted for when they are crossed in the process of enlarging the contour and wrapping it around the cut. The result is:

$`{\displaystyle \underset{4\mathrm{images}}{}}\mu _ip_iz_i^m\overline{z}_i^n`$ $`=`$ $`{\displaystyle \frac{B[P_0\zeta Q\overline{\zeta }]^m[P_0\overline{\zeta }\overline{Q}\zeta ]^n}{q^{2m+2n}[1\gamma _1^2\gamma _2^2]^{m+n+1}}}`$ (74)
$``$ $`q^2B{\displaystyle \underset{\mathrm{poles}}{}}\mathrm{residues}\mathrm{of}\left\{{\displaystyle \frac{[P\zeta Q\overline{\zeta }]^m[P\overline{\zeta }\overline{Q}\zeta ]^n}{tK(t;\zeta ,\overline{\zeta })[P^2|Q|^2]^{m+n1}}}\right\}+S_{m,n}(B),`$
where

$$S_{m,n}(B)=Bq^2\underset{j=1}{\overset{\mathrm{}}{}}\frac{\mathrm{sin}(Bj\pi )}{\pi }_0^{\mathrm{}}\frac{dww^{jB1}[\lambda (w+t_2)+\nu (w+t_1)]^{m+j}[\overline{\lambda }(w+t_2)+\overline{\nu }(w+t_1)]^{n+j}}{[(w+t_1)(w+t_2)]^{2j+m+n+1}}.$$
(75)
The binomial expansion of the integrand is less simple than in §4.1, as $`mn`$ and a double summation is needed. The result is the following homogeneous polynomial of $`\lambda `$, $`\overline{\lambda }`$, $`\nu `$, and $`\overline{\nu }`$:

$$S_{m,n}(B)=Bq^2\underset{j=1}{\overset{\mathrm{}}{}}\underset{\mathrm{}=0}{\overset{2j+m+n}{}}I_{j,\mathrm{},m+n}\underset{k=\mathrm{max}(0,\mathrm{}nj)}{\overset{k=\mathrm{min}(\mathrm{},m+j)}{}}\left(\genfrac{}{}{0pt}{}{m+j}{k}\right)\left(\genfrac{}{}{0pt}{}{n+j}{\mathrm{}k}\right)\lambda ^k\overline{\lambda }^\mathrm{}k\nu ^{m+jk}\overline{\nu }^{n+j+k\mathrm{}}.$$
(76)
The $`j`$th term is now a homogeneous polynomial of degree $`2j+m+n`$ in the source coordinates. A more general integral with an extra subscript (which was zero previously) now appears as a coefficient. It is

$`I_{j,\mathrm{},N}(t_1,t_2)`$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}Bj\pi }{\pi }}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{w^{jB1}dw}{(w+t_1)^{\mathrm{}+1}(w+t_2)^{2j+N\mathrm{}+1}}},0B2`$
$`=`$ $`{\displaystyle \frac{(1)^{N+1}}{(2j+N+1)!}}{\displaystyle \underset{s=1}{\overset{2j+N+1}{}}}(jBs){\displaystyle \frac{t_2^{jB+\mathrm{}2jN1}}{t_1^{\mathrm{}+1}}}{}_{2}{}^{}F_{1}^{}(\mathrm{}+1,jB;2j+2+N;1{\displaystyle \frac{t_2}{t_1}}).`$
Although the integral definition of $`I_{j,\mathrm{},N}`$ is valid only if $`B2`$, the hypergeometric formula is valid for all $`B>1`$. Explicit closed form expressions for all the hypergeometric functions are given in Appendix C. Using these results allows us to extend the formulae given in §3.2 and §3.3 to configuration invariants. The larger $`m+n`$ is, the longer they are, so we shall here quote only the first reciprocal moments, which are

$`{\displaystyle \underset{\mathrm{images}}{}}{\displaystyle \frac{\mu _ip_ix_i}{(x_i^2+y_i^2q^2)}}`$ $`=`$ $`{\displaystyle \frac{Bq^2}{2}}\left[\left[\left(1{\displaystyle \frac{\mathrm{Re}Q}{2|Q|}}\right)I_{1,0,1}+\left(1+{\displaystyle \frac{\mathrm{Re}Q}{2|Q|}}\right)I_{1,1,1}\right]\xi +{\displaystyle \frac{\gamma _2q^3\eta }{|Q|}}[I_{1,1,1}I_{1,0,1}]\right],`$
$`{\displaystyle \underset{\mathrm{images}}{}}{\displaystyle \frac{\mu _ip_iy_i}{(x_i^2+y_i^2q^2)}}`$ $`=`$ $`{\displaystyle \frac{Bq^4}{2}}\left[\left[\left(1+{\displaystyle \frac{\mathrm{Re}Q}{2|Q|}}\right)I_{1,0,1}+\left(1{\displaystyle \frac{\mathrm{Re}Q}{2|Q|}}\right)I_{1,1,1}\right]\eta +{\displaystyle \frac{\gamma _2\xi }{|Q|}}[I_{1,1,1}I_{1,0,1}]\right],`$
where

$`I_{1,0,1}`$ $`=`$ $`{\displaystyle \frac{t_1t_2^{B2}}{4|Q|^2}}\left[(B2)\left({\displaystyle \frac{t_2}{t_1}}\right)+\left({\displaystyle \frac{t_2}{t_1}}\right)^{2B}B+1\right],`$
$`I_{1,1,1}`$ $`=`$ $`{\displaystyle \frac{t_2^{B1}}{4|Q|^2}}\left[1+\left({\displaystyle \frac{t_2}{t_1}}\right)^{1B}\left[B2+(B1)\left({\displaystyle \frac{t_2}{t_1}}\right)\right]\right].`$ (78)
The first order corrections for the first moments are given in Appendix C.

### 4.3 Perturbation Analysis

Invariants, and corrections to them, may also be derived by using perturbation expansions to solve the imaging equation. Its four main roots tend in pairs to $`t_1`$ and $`t_2`$ as the source tends to the center of the lensing galaxy, and the expansions are in integer powers of $`\zeta `$ and $`\overline{\zeta }`$. One pair of images lies at

$$t=t_1\pm t_1^{B/2}|\lambda |+\frac{1}{2}Bt_1^{B1}|\lambda |^2\pm \frac{1}{2}t_1^{3B/2}|\lambda |\left[\frac{B(3B2)|\lambda |^2}{t_1^2}+\frac{|\nu |^2}{4|Q|^2}\right]+O(|\zeta |^4),$$
(79)
and the other at

$$t=t_2\pm t_2^{B/2}|\nu |+\frac{1}{2}Bt_2^{B1}|\nu |^2\pm \frac{1}{2}t_2^{3B/2}|\nu |\left[\frac{B(3B2)|\nu |^2}{t_2^2}+\frac{|\lambda |^2}{4|Q|^2}\right]+O(|\zeta |^4),$$
(80)
where $`\lambda `$ and $`\nu `$ are as defined in equation (61). The signed magnifications for the first pair of roots are, from equation (25),

$`\mu p`$ $`=`$ $`{\displaystyle \frac{q^2B}{4|Q|t_1}}\pm {\displaystyle \frac{Bq^2|\lambda |t_1^{(B/2)1}}{8|Q|}}\left[{\displaystyle \frac{1}{|Q|}}{\displaystyle \frac{B2}{t_1}}\right]`$ (81)
$`+{\displaystyle \frac{Bq^2t_1^{B2}}{8|Q|^2}}\left[(B1)|\lambda |^2\left[1{\displaystyle \frac{(B2)|Q|}{t_1}}\right]{\displaystyle \frac{t_1|\zeta |^2}{2|Q|}}\right],`$
while those for the second pair can be found from them by applying the transformations $`t_1t_2`$, $`\lambda \nu `$, and $`|Q||Q|`$. Summation then gives equation (69) again, as it should.

Finding configuration invariants by perturbation expansion requires the additional effort of expanding for the image positions $`z`$ and $`\overline{z}`$. Though the individual terms are of some intrinsic interest, perturbation analysis is an unnecessarily laborious way of computing invariants now that we have the contour integrals which compute them to all orders.

## 5 Astrophysical Applications

There is inevitable uncertainty in the distribution of mass in lenses, and all methods (whether parametrised fitting or non-parametric modelling) must make strong assumptions as to the mass distribution before they can make progress. The major application of lensing invariants is to short-cut the modelling process. In §5.1, we provide simple tests to determine whether a given lens system is well modelled by an elliptic power-law potential with shear. Next, we use the lensing invariants to develop an algorithm for estimation of the model parameters and apply it to fake data in §5.2 and real data for the Einstein Cross (G2237+0305) in §5.3. The work in §4 has demonstrated that the simple equations for the moments we derived in §3 – namely, eqns (34) - (36), (39) - (42), and (44) – are excellent approximations for all $`B`$. Here, we shall assume that these results are exact.

### 5.1 Tests for Elliptic Power-Law Lenses

From an observational perspective, only the flux ratios of the images and not their absolute magnifications are measured. The flux ratios are defined as $`f_i=\mu _i/\mu _1`$. So, they are normalised to the flux of the reference image $`\mu _1`$, which most observers generally take as the brightest image. For the quadruple lenses, there are three flux ratios $`f_2,f_3`$ and $`f_4`$, with $`f_1=1`$ by definition. So, the moments of the flux ratios are:

$$F_{k,x}=\underset{i=1}{\overset{4}{}}f_ip_ix_i^k,F_{k,y}=\underset{i=1}{\overset{4}{}}f_ip_iy_i^k.$$
(82)
In other words, $`F_0=F_{0,x}=F_{0,y}`$ is just the magnification invariant normalised by the flux of the reference image $`\mu _1`$. Similarly, the higher moments of the flux ratios are given by dividing the configuration invariants by $`\mu _1`$.

The source position is just

$$\xi =(1+\gamma _1)\frac{F_{1,x}}{F_0}+\gamma _2\frac{F_{1,y}}{F_0},\eta =\gamma _2\frac{F_{1,x}}{F_0}+(1\gamma _1)\frac{F_{1,y}}{F_0}.$$
(83)
Substituting these into the first moments (36) and the reciprocal first moments (44), we obtain

$$\frac{F_{1,x}F_{1,x}}{F_0^2}=1,\frac{F_{1,y}F_{1,y}}{F_0^2}=1.$$
(84)
Provided the lensing galaxy can be identified and the position angle of its major axis can be measured, then all quantities in eqn (84) are directly available from the data, So, these equations provide two simple, independent checks that are easy to apply. If they are satisifed, then the lens may be well-represented by an elliptic power-law potential with shear. If the major axis of the lensing galaxy is not known, then one of equations can be used to estimate the position angle and the second used to check that the elliptic power-law potential is a reasonable model.

The main difficulty in application of the lensing invariants is the uncertainty in the observable quantities. The relative positions of the lensing galaxy and the images may be accurate to better than $`1\%`$. In propitious circumstances, the relative fluxes may be good to perhaps $`5\%`$. Unhappily, microlensing may render these fluxes much more uncertain, with the error perhaps even as high as $`50\%`$ in the optical (Schechter 2000). Microlensing is of course valuable, as it sets constraints on the mass fraction in compact objects in the haloes of the lensing galaxies (e.g., Witt, Mao & Schechter 1995). However, for all researchers modelling the smooth, large-scale distribution of mass in the lensing galaxy, microlensing is a serious encumbrance. It is sometimes claimed that microlensing is primarily a difficulty in the optical wavebands and may be circumvented using radio fluxes. However, there have been very recent identifications of microlensing events in the radio (Koopmans & de Bruyn 2000), and so even the radio fluxes may be more uncertain than is widely believed. As benchmark values for our calculations, we take $`5\%`$ and $`15\%`$ errors. In some cases, the radio fluxes may already be accurate to $`5\%`$; in other cases, this accuracy will be achievable in the near future, once the effects of microlensing are better understood and quantified. So, $`15\%`$ is a conservative value for the errors in present-day radio flux measurements.

Figure 3 shows histograms of the distribution of the ratios (84) in 400 cases, allowing for $`5\%`$ and $`15\%`$ errors in the flux measurements. In both cases, the distributions are peaked around unity, but there is substantial scatter. Even for datasets with flux errors of $`5\%`$, the test should be used in the form

$$\frac{F_{1,x}F_{1,x}}{F_0^2}1\pm 2,\frac{F_{1,y}F_{1,y}}{F_0^2}1\pm 2.$$
(85)
for practical application. If the flux errors are as high as $`15\%`$, the tests become much less reliable. If a lens system passes this test (as for example the Einstein Cross does), then it is worthwhile to proceed to estimation of the lens parameters.

### 5.2 Monte Carlo Simulations

There are three galaxy model parameters, namely $`A,q`$ and $`\beta `$. In the general case, there are two unknown components of shear $`\gamma _1`$ and $`\gamma _2`$. Finally, there are three more unknowns, namely the source offsets $`\xi `$ and $`\eta `$ and the position angle $`\mathrm{\Theta }`$ of the lensing galaxy. It is possible to solve formally for the model parameters, given the magnification invariant and the first, second and third moments. However, this is not a useful way too proceed as small errors in the observables amplify into larger errors in the first, second and third moments. This is exacerbated by the plus and minus pattern imposed by the parities. The fitting of lens models to data normally proceeds by the minimisation of the differences between the observed and true positions (and flux ratios, if desired) using a downhill simplex routine to search parameter space (e.g., Kochanek 1991). The lensing invariants can be fruitfully exploited to reduce the dimensionality of the parameter space in which the search proceeds.

As an example, let us consider a simple case. Suppose that we assume that the lensing galaxy has negligible external shear. First, it is straightfoward to find the source offsets

$$\xi =\frac{F_0}{F_{1,x}},\eta =\frac{F_0}{F_{1,y}}.$$
(86)
From the reciprocal second moments, it is then easy to see that the flattening of the potential is

$$q^2=\frac{F_{1,y}}{F_{1,x}}\sqrt{\frac{F_{2,x}F_0F_{1,x}^2}{F_{2,y}F_0F_{1,y}^2}}.$$
(87)
From the second moments, we obtain

$$q^{2B}=\frac{F_{2,x}F_0F_{1,x}^2}{F_{2,y}F_0F_{1,y}^2}.$$
(88)
from which $`B`$ can be estimated. There remain an unknown position angle $`\mathrm{\Theta }`$ of the lensing galaxy and an unknown potential normalisation, which can be obtained by minimizing

$$\chi ^2(A,\mathrm{\Theta })=\underset{\mathrm{images}}{}\left[\xi x_i(1\frac{A}{(x_i^2+y_i^2q^2)^{1\beta /2}})\right]^2+\left[\eta y_i(1\frac{Aq^2}{(x_i^2+y_i^2q^2)^{1\beta /2}})\right]^2$$
(89)
Figure 3 shows the results of Monte Carlo simulations using fake data. Taking $`A=1`$, $`B=2`$ and $`q=0.75`$, the lens equation is first solved to find the true image positions and fluxes for random source offsets. A random rotation is performed to take the image positions from the intinsic coordinate system of the lensing galaxy to the observer’s coordinate system on the plane of the sky. The fluxes are adusted by random errors of $`5\%`$ and $`15\%`$ to give the simulated data. This is fed into the algorithm and the random position angle $`\mathrm{\Theta }`$ and the parameters of the lensing galaxy $`A,B`$ and $`q`$ are recovered. Figure 4 shows histograms of 200 such trials. We see that if the relative fluxes are accurate to $`5\%`$ or even $`15\%`$, then there is good recovery of the model parameters $`q`$ and $`A`$. The parameter controlling the fall-off the potential $`B`$ seems more vulnerable to errors, as there are tails to the distribution for the case of $`15\%`$ errors. Tests also show that essentially no useful information can be extracted if the fluxes are only accurate to $`50\%`$.

### 5.3 The Einstein Cross

Table 1 presents the available astrometric data on the Einstein Cross quadruple system from Crane et al. (1991). The astrometry is provided in terms of the relative right ascension and declination (in arcseconds) with respect to the center of the lensing galaxy. The radio data on the fluxes of the four images A,B,C and D is given, rather than the optical data which is known to be strongly effected by microlensing (e.g., Ostensen et al. 1994). The bulge of the lensing galaxy has a measured position angle of $`68^{}`$ on the sky (Rix et al. 1992). If the misalignment of the source is not too great, it is natural to expect the images of negative parity are located roughly parallel to the major axis, those of positive parity are roughly perpendicular to the major axis. Accordingly, we reckon that the negative parity images are C and D, whilst the positive parity ones are A and B.

To fit an elliptical power-law model in the presence of external on-axis shear is straightforward. We use the relations

$$\xi =\frac{F_0(1+\gamma _1)}{F_{1,x}},\eta =\frac{F_0(1\gamma _1)}{F_{1,y}},q^2\left(\frac{1\gamma _1}{1+\gamma _1}\right)=\frac{F_{1,y}}{F_{1,x}}\sqrt{\frac{F_{2,x}F_0F_{1,x}^2}{F_{2,y}F_0F_{1,y}^2}},$$
(90)
to minimise the $`\chi ^2`$

$$\chi ^2(\gamma _1,A,\beta ,\mathrm{\Theta })=\underset{\mathrm{images}}{}\left[\xi x_i(1+\gamma _1\frac{A}{(x_i^2+y_i^2q^2)^{1\beta /2}})\right]^2+\left[\eta y_i(1\gamma _1\frac{Aq^2}{(x_i^2+y_i^2q^2)^{1\beta /2}})\right]^2$$
(91)
over $`\gamma _1`$, $`A`$, $`\mathrm{\Theta }`$. The results for the model parameters are given in Table 2. First, let us note that the recovered position angle agrees well with the observed position angle of the major axis of the lensing galaxy. Second, only mild on-axis shear is required. The flattening of the lensing galaxy is significant, but not unreasonable. Third, the lensing galaxy is not close to isothermal, and there are no good fits with $`\beta =1`$. As all the images are at almost the same projected distance, the density fall-off, and hence the parameter $`\beta `$, is not well-constrained. The total mass enclosed within the images is more robust and comparable to the values found by previous investigators (e.g., Rix et al. 1992). Table 3 shows how well the best fit reproduces the data. The positions and the relative brightnesses of the four images are reasonably well reproduced. There are mild discrepancies in the magnifications of images B and D, which are both somewhat fainter in our model. The computational advantage of the lensing invariants (90) is that the $`\chi ^2`$ minimisation can proceed in fewer dimensions, with obvious gains in accuracy and speed.

Note that Witt & Mao (2000) give a cautionary example of a model with a flat rotation curve and with a projected density distribution stratified on ellipses. Although superficially similar to the elliptic power-law models, the lensing invariants are quite different. For example the sum of the signed magnifications is $`2.8`$ rather than 2. The magnifications depend on the second derivatives of the lensing potential and so are more susceptible to the details of the modelling than the image positions. Real lenses certainly do exhibit diverse behaviour. Some – for example, B1422+231 where the lens is both a galaxy and a cluster or B1608+656 where the lens is a pair of interacting galaxies – are not likely to be well described by any simple models. Nonetheless, much of the optical depth is in isolated giant elliptical galaxies or in spirals dominated by dark haloes. It is reasonable to model such lenses with simple potentials and the algorithms we have developed have a real rôle to play.

## 6 Conclusions

This paper has presented a powerful new method for obtaining lensing invariants. Such invariants are calculable from the observables and they remain unchanged as the source of radiation is moved within the central caustic. The magnification invariant is the sum of the signed magnifications of the images, whereas the configuration invariants are sums over powers (either positive or negative) of the image positions multiplied by the signed magnifications. The idea behind the method is to use Cauchy’s theorem and the residue calculus to recast sums over the images as contour integrals.

We have illustrated the method by application to a simple family of galaxy models, the power-law models in which the gravitational potential is scale-free and stratified on similar concentric ellipsoids. They may be embedded in an external shear field of arbitrary orientation. The convergence $`\kappa `$ of the models falls off with distance like a power-law with index $`2/B`$. The lensing invariants are exact for the models with $`B=1`$ (the point mass case), $`B=2`$ (an isothermal galaxy with projected potential $`\psi r`$) and $`B=3`$ (a galaxy with $`\psi r^{4/3}`$). Some of the lensing invariants for the power-law models were previously calculated by Witt & Mao (2000) for the case of on-axis shear. Our results complement theirs by extending the invariants into the régime of arbitrary shear, by giving general and more extensive formulae, as well as by providing some entirely new invariants (the reciprocal moments).

In their discovery papers, Witt & Mao (1995, 2000) used direct elimination methods to find the first examples of lensing invariants. We believe that our contour integral method offers considerable advantages over this. All that is needed is the identification of the poles and the branch points of the lensing equation and the calculation of the residues using Cauchy’s theorem. Our method is simpler than the multi-dimensional residue calculus proposed by Dalal & Rabin (2000). Nevertheless, the idea motivating Dalal & Rabin’s work – namely, that sums over images may be replaced by contour integrals – is the same as that behind our paper. Only the familiar single-variable calculus is required in our formulation of the problem, but it has allowed us to go beyond the polynomial equations and meromorphic functions studied by Dalal & Rabin, and to handle branch points. In turn, this led us to our discovery of full series expansions for the lensing invariants of the four bright sources, unsullied by spurious roots of the imaging equation.

We have shown how the lensing invariants may be exploited in modelling by using Monte Carlo simulations to demonstrate the feasibility of the recovery of parameters, such as the flattening and the normalisation of the potential well. Most lens fitting is conventionally done by $`\chi ^2`$ minimisation. The lensing invariants play an important rôle by reducing the dimensionality of the space in which the global minimum is sought, with consequent gains in accuracy and speed. The power-law models with shear are sensible ones to use when the lens is a reasonably isolated spiral or elliptical galaxy. To illustrate this, we have used the lensing invariants to provide a good fit to the astrometry and photometry of the Einstein Cross (G2237+0305).

Of course, caution is needed in application of the lensing invariants identified in this paper, as they only apply to a particular family of galaxy models. The pressing problem is to calculate the lensing invariants for more general potentials, such as galaxy models with cores or density distributions stratified on ellipses (e.g. Kassiola & Kovner 1993, Barkana 1998). This will lead to a greater understanding of how the invariants can constrain the mass distribution for real multiply-lensed systems. We believe that our contour integral method can be adapted to these cases and are presently pursuing further investigations.

We are indebted to Hans Witt for useful conversations and much encouragement. NWE thanks the Royal Society for financial support. CH thanks the sub-department of Theoretical Physics, University of Oxford for hospitality during a stimulating six-month sabbatical visit, and the Florida State University for awarding that sabbatical. His work is supported in part by NSF through grant DMS-9704615.

## Appendix A Images and Caustics

When the source is sufficiently close to the center of the lens and $`|\zeta |`$ small, then there are four images if $`B2`$ and five images if $`B>2`$. To see this, we split $`K`$ into two components $`K=K_1+K_2`$, as defined in eqn (60). We recall that, provided that $`\gamma _1^2+\gamma _2^2<1`$, both $`t_1`$ and $`t_2`$ (defined in eqn (37)) are positive with $`t_1>t_2`$ and $`t_2<1`$. As Figure 5 shows, the graph of $`K_1`$, with its two zero minima at $`t=t_1`$ and $`t=t_2`$, has two nearby pairs of intersections with the positive $`K_2`$, which is small when $`|\zeta |`$ is small. There is a fifth image for large positive $`t`$ when $`B>2`$ because, however small $`|\zeta |`$, the asymptotic $`t^{B+2}|\zeta |^2`$ growth of $`K_2`$ eventually overtakes the $`t^4`$ growth of $`K_1`$. This fifth root is given by the $`n=0`$ case of equation (53).

The number of images changes when pairs of roots of the imaging equation merge and become double. Both $`K`$ and $`K_t`$ vanish at double roots, and then

$$\frac{1}{K_2}\frac{K_2}{t}+\frac{1}{K_1}\frac{K_1}{t}=\frac{B}{t}+\frac{2}{P|Q|}+\frac{2}{P+|Q|}\frac{\zeta }{P\zeta Q\overline{\zeta }}\frac{\overline{\zeta }}{P\overline{\zeta }\overline{Q}\zeta }=0.$$
(A1)
Because each denominator is linear in $`t`$, equation (A1) gives a quartic in $`t`$, degenerating to a cubic in the special case of $`B=2`$. Moreover, the equation is independent of $`|\zeta |`$, and depends only on the angular argument $`\varphi `$ of the complex source position $`\zeta =|\zeta |e^{\mathrm{i}\varphi }`$. Hence the source positions which give double images are found by looking for positive roots of equation (A1) along each specific angular direction. We demonstrate below that, with our restrictions on the lens parameters, eqn (A1) has one real root in $`[t_2,t_1]`$ for all $`B`$, one real root in $`[t_1,\mathrm{})`$ for $`B>2`$, and no others. Once $`t`$ is found the imaging equation gives $`|\zeta |^2=[P^2|Q|^2]^2/[t^B(PQe^{2\mathrm{i}\varphi })(P\overline{Q}e^{2\mathrm{i}\varphi })]`$. The root in $`[t_2,t_1]`$ corresponds to the one point on the tangential caustic in the source direction $`\varphi `$, and the root in $`[t_1,\mathrm{})`$ for $`B>2`$ corresponds to the one point on the radial caustic in that direction. These roots must generally be found numerically, but there are four special directions for which there are analytical formulas. These are the directions for which $`Qe^{2\mathrm{i}\varphi }`$ is real. When $`e^{2\mathrm{i}\varphi }=Q/|Q|`$, then $`t=t_2`$ where $`P=|Q|`$ is a triple root of the quartic, and gives a cusp at $`|\zeta |=2|Q|/t_2^{B/2}`$ on the tangential caustic. The remaining root is $`t=Bt_1/(B2)`$, which for $`B>2`$ gives a point $`|\zeta |=(2/B)[(B2)/Bt_1]^{(B2)/2}`$ on the radial caustic. When $`e^{2\mathrm{i}\varphi }=Q/|Q|`$, then $`t=t_1`$ where $`P=|Q|`$ is a triple root of the quartic, and $`t=Bt_2/(B2)`$ is the fourth root. The former gives a cusp on the tangential caustic at which $`|\zeta |=2|Q|/t_1^{B/2}`$ for $`B2`$, and also for $`B>2`$ if $`t_1`$ is the smaller of the roots. In the latter case, the fourth root gives a point $`|\zeta |=(2/B)[(B2)/Bt_2]^{(B2)/2}`$ on the radial caustic. However, if $`Bt_2/(B2)<t_1`$, then the latter point is on the tangential caustic, while $`|\zeta |=2|Q|/t_1^{B/2}`$ on the radial one. Depending on the parameters, either caustic may lie wholly inside the other, or they may intersect as in the case shown in Figure 6.

In the general case, $`Q\overline{\zeta }/\zeta =|Q|e^{\mathrm{i}\theta }`$ for some angle $`\theta `$ which is not an integer multiple of $`\pi `$. To find the roots of equation (A1) for $`1<\mathrm{cos}\theta <1`$, we rewrite that equation as

$$\frac{B}{t}=\frac{2(P|Q|\mathrm{cos}\theta )}{P^2+|Q|^22P|Q|\mathrm{cos}\theta }+\frac{2}{tt_1}+\frac{2}{tt_2}=\frac{2P^3+6P^2|Q|\mathrm{cos}\theta 6P|Q|^2+2|Q|^3\mathrm{cos}\theta }{(P^2|Q|^2)(P^2+|Q|^22P|Q|\mathrm{cos}\theta )}.$$
(A2)
Graphs of its two sides as functions of $`t`$ are shown in Figure 7. The right hand side has vertical asymptotes at $`t=t_2`$ and $`t=t_1`$, and a zero in between. This is the only zero of the numerator because its derivative with respect to $`t`$, which is the negative of that with respect to $`P`$, is $`6[P|Q|\mathrm{cos}\theta ]^2+6|Q|^2\mathrm{sin}^2\theta `$ and hence always positive. The two sides of equation (A2) have different signs and hence no roots in $`(0,t_2))`$. They have at least one root in $`(t_2,t_1)`$, and another in $`(t_1,\mathrm{})`$ for $`B>2`$ because the right hand side, which decays as $`2/t`$ is ultimately the smaller for large $`t`$. To prove that these are the only possibilities, we rewrite equation (A2) as a quartic in two different ways, and apply Descartes’ rule of signs (Henrici 1974). We also reduce the number of parameters to three by working with the variable $`\sigma =(tt_1)/|Q|=1P/|Q|`$. The first form of the quartic is

$`(B2)\sigma ^4`$ $`+`$ $`2[(B2)(2+\mathrm{cos}\theta )(p_0+\mathrm{cos}\theta )]\sigma ^3`$ (A3)
$`+`$ $`2(1+\mathrm{cos}\theta )[3(Bp_03)\sigma ^2+2(B3p_05)\sigma 4(1+p_0)]=0,`$
where $`p_0=P_0/|Q|>1`$. Let $`c_n`$ denotes the coefficient of $`\sigma ^n`$. Each coefficient is manifestly either zero or negative for $`B2`$, and hence there can then be no roots with $`\sigma >0`$ i.e., $`t>t_1`$. For $`B>2`$, $`c_4>0`$, $`c_0<0`$, and $`c_2c_1=2(1+\mathrm{cos}\theta )(B+1+3p_0)>0`$ and so $`c_2>c_1`$. We now show that, regardless of the sign of $`c_2`$, there can be only one sign change from positive to negative in the sequence of coefficients $`(c_4,c_3,c_2,c_1,c_0)`$. If $`c_20`$, then $`Bp_0+3`$ and $`c_34+2p_0(1+\mathrm{cos}\theta )>0`$, and the single sign change occurs somewhere between $`c_2`$ and $`c_0`$. If $`c_2<0`$, there is a single sign change somewhere between $`c_4`$ and $`c_2`$, depending on the sign of $`c_3`$. By Descartes’ rule of signs, the single sign change in the sequence of coefficients means that the quartic has exactly one root in $`t_1>1`$.

To show that there is one and only one root of the quartic in $`(t_1,t_2)`$, we use a bilinear transformation to the variable $`\tau =(tt_2)/(t_1t)=(|Q|P)/(|Q|+P)`$ which maps the interval $`t_2<t<t_1`$ to the whole of $`\tau >0`$. Equation (A2) remains a quartic, and becomes

$`(1+\mathrm{cos}\theta )[(1+p_0)\tau ^4+(B+p_01)\tau ^3](1\mathrm{cos}\theta )[(1+p_0B)\tau +(p_01)]=0.(\mathrm{A4})`$
Now $`c_4>0`$, $`c_3>0`$, $`c_2=0`$, and $`c_0<0`$, when $`c_n`$ now denotes the coefficient of $`\tau ^n`$. There is a single sign change from positive to negative in the coefficient sequence regardless of the sign of $`c_1`$, and hence, by Descartes’ rule of signs, always one root for $`t`$ in $`(t_2,t_1)`$. This root gives a point on the tangential caustic.

## Appendix B Third Order Moments

Witt & Mao (2000) give the third order moments in the case of on-axis shear only. The third order moments, derived from the residues of simple poles at $`t=0`$, and double poles at $`t=t_1`$ and $`t=t_2`$ are

$`{\displaystyle \mu _ip_ix_i^3}`$ $`=`$ $`{\displaystyle \frac{1}{4}}(R_1+3R_2)+{\displaystyle \frac{B[(1\gamma _1)\xi \gamma _2\eta ]^3}{(1\gamma _1^2\gamma _2^2)^4}},`$ (B1)
$`{\displaystyle \mu _ip_ix_i^2y}`$ $`=`$ $`{\displaystyle \frac{q}{4}}(I_1+I_2)+{\displaystyle \frac{B[(1\gamma _1)\xi \gamma _2\eta ]^2[(1+\gamma _1)\eta \gamma _2\xi ]}{(1\gamma _1^2\gamma _2^2)^4}},`$ (B2)
$`{\displaystyle \mu _ip_ix_iy_i^2}`$ $`=`$ $`{\displaystyle \frac{q^2}{4}}(R_2R_1)+{\displaystyle \frac{B[(1\gamma _1)\xi \gamma _2\eta ][(1+\gamma _1)\eta \gamma _2\xi ]^2}{(1\gamma _1^2\gamma _2^2)^4}},`$ (B3)
$`{\displaystyle \mu _ip_iy_i^3}`$ $`=`$ $`{\displaystyle \frac{q^3}{4}}(3I_2I_1)+{\displaystyle \frac{B[(1+\gamma _1)\eta \gamma _2\xi ]^3}{(1\gamma _1^2\gamma _2^2)^4}}.`$ (B4)
The $`R_1`$, $`R_2`$, $`I_1`$, and $`I_2`$ terms are all linear in the source coordinates. They are the real and imaginary terms of the $`m=3`$, $`n=0`$ and $`m=2`$, $`n=1`$ configuration moments respectively, and are

$`R_1`$ $`=`$ $`{\displaystyle \frac{q^2B}{4|Q|^3}}\{\{{\displaystyle \frac{1}{t_1^{B+1}}}[1{\displaystyle \frac{(B+1)|Q|}{t_1}}]{\displaystyle \frac{1}{t_2^{B+1}}}[1+{\displaystyle \frac{(B+1)|Q|}{t_2}}]\}[\mathrm{Re}(Q)\xi \gamma _2q^2\eta ]`$
$`\left\{{\displaystyle \frac{1}{t_1^{B+1}}}\left[2+{\displaystyle \frac{(B+1)|Q|}{t_1}}\right]+{\displaystyle \frac{1}{t_2^{B+1}}}\left[2{\displaystyle \frac{(B+1)|Q|}{t_2}}\right]\right\}`$
$`\times {\displaystyle \frac{[(\mathrm{Re}(Q))^2\xi +2\gamma _2q^2\mathrm{Re}(Q)\eta \gamma _2^2q^2\xi ]}{|Q|}}\},`$
$`R_2`$ $`=`$ $`{\displaystyle \frac{q^2B}{4|Q|^2}}\{(B+1)|Q|\xi ({\displaystyle \frac{1}{t_2^{B+2}}}{\displaystyle \frac{1}{t_1^{B+2}}})`$
$`[(B+1)({\displaystyle \frac{1}{t_2^{B+2}}}+{\displaystyle \frac{1}{t_1^{B+2}}})+{\displaystyle \frac{1}{|Q|}}({\displaystyle \frac{1}{t_1^{B+1}}}{\displaystyle \frac{1}{t_2^{B+1}}})][\mathrm{Re}(Q)\xi +\gamma _2q^2\eta ]\},`$
$`I_1`$ $`=`$ $`{\displaystyle \frac{q^3B}{4|Q|^3}}\{\{{\displaystyle \frac{1}{t_1^{B+1}}}[1{\displaystyle \frac{(B+1)|Q|}{t_1}}]{\displaystyle \frac{1}{t_2^{B+1}}}[1+{\displaystyle \frac{(B+1)|Q|}{t_2}}]\}[\mathrm{Re}(Q)\eta +\gamma _2\xi ]`$
$`+\left\{{\displaystyle \frac{1}{t_1^{B+1}}}\left[2+{\displaystyle \frac{(B+1)|Q|}{t_1}}\right]+{\displaystyle \frac{1}{t_2^{B+1}}}\left[2{\displaystyle \frac{(B+1)|Q|}{t_2}}\right]\right\}`$
$`{\displaystyle \frac{[(\mathrm{Re}(Q))^2\eta 2\gamma _2\mathrm{Re}(Q)\xi \gamma _2^2q^2\eta ]}{|Q|}}\},`$
$`I_2`$ $`=`$ $`{\displaystyle \frac{q^3B}{4|Q|^2}}\{(B+1)|Q|\eta ({\displaystyle \frac{1}{t_2^{B+2}}}{\displaystyle \frac{1}{t_1^{B+2}}})`$ (B5)
$`+[(B+1)({\displaystyle \frac{1}{t_2^{B+2}}}+{\displaystyle \frac{1}{t_1^{B+2}}})+{\displaystyle \frac{1}{|Q|}}({\displaystyle \frac{1}{t_1^{B+1}}}{\displaystyle \frac{1}{t_2^{B+1}}})][\mathrm{Re}(Q)\eta \gamma _2\xi ]\}.`$
## Appendix C Formulae for Hypergeometric Functions

The hypergeometric functions that occur in this paper are all of the form $`{}_{2}{}^{}F_{1}^{}(M,jB;M+N;1z)`$, where $`M`$ and $`N`$ are integers and both are positive. The following formula expresses them as two finite sums, one with $`N`$ terms and the other with $`M`$ terms. It is obtained by applying first Abramowitz and Stegun (1965, hereafter AS) equation (15.3.9) and then their equation (15.3.5). The result is

$`{\displaystyle \frac{(1z)^{M+N1}}{(M+N1)!}}{}_{2}{}^{}F_{1}^{}(M,jB;M+N;1z)`$ $`=`$ $`{\displaystyle \frac{(1)^M}{(N1)!}}{\displaystyle \underset{k=0}{\overset{N1}{}}}{\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{N1}{k}\right)(z)^k}{_{s=Nk}^{M+Nk1}(jBs)}}`$ (C1)
$`+{\displaystyle \frac{z^{NjB}}{(M1)!}}{\displaystyle \underset{k=0}{\overset{M1}{}}}{\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{M1}{k}\right)(z)^k}{_{s=1+k}^{N+k}(jBs)}}.`$
It is significant that the products of terms in $`B`$ in the denominators are all subparts of a longer product of terms in $`B`$ which multiplies the hypergeometric functions in the definition (4.2) of the coefficients $`I_{j,\mathrm{},N}`$ which appear in our series for the lensing invariants. Hence, these coefficients never have poles in $`B`$ and are all analytic functions of $`B`$. For example, the two $`N=0`$ integrals needed for the evaluation of the extra term in equation (69) are:

$`I_{1,2}`$ $`=`$ $`{\displaystyle \frac{t_1^{B1}}{16|Q|^3}}\left[(B2)(B3)+2(B1)(B3)\left({\displaystyle \frac{t_2}{t_1}}\right)(B1)(B2)\left({\displaystyle \frac{t_2}{t_1}}\right)^2+2\left({\displaystyle \frac{t_2}{t_1}}\right)^{B1}\right],`$
$`I_{1,0}`$ $`=`$ $`{\displaystyle \frac{t_1^2t_2^{B3}}{16|Q|^3}}\left[(B1)(B2)2(B1)(B3)\left({\displaystyle \frac{t_2}{t_1}}\right)+(B2)(B3)\left({\displaystyle \frac{t_2}{t_1}}\right)^22\left({\displaystyle \frac{t_2}{t_1}}\right)^{3B}\right].`$
The similarity of these two expressions is an instance of the general symmetry relation:

$$I_{j,\mathrm{},N}(t_1,t_2)=I_{j,2j\mathrm{},N}(t_2,t_1).$$
(C2)
Hence, only two distinct $`I`$ terms are needed for the first order corrections for the first moments

$`{\displaystyle \underset{4\mathrm{images}}{}}\mu _ip_ix_i`$ $`=`$ $`{\displaystyle \frac{B((1\gamma _1)\xi \gamma _2\eta )}{(1\gamma _1^2\gamma _2^2)^2}}+K_{30}^{10}\xi ^3+K_{21}^{10}\xi ^2\eta +K_{12}^{10}\xi \eta ^2+K_{03}^{10}\eta ^3+O(|\zeta |^5),`$ (C3)
$`{\displaystyle \underset{4\mathrm{images}}{}}\mu _ip_iy_i`$ $`=`$ $`{\displaystyle \frac{B((1+\gamma _1)\eta \gamma _2\xi )}{(1\gamma _1^2\gamma _2^2)^2}}+K_{30}^{01}\xi ^3+K_{21}^{01}\xi ^2\eta +K_{12}^{01}\xi \eta ^2+K_{03}^{01}\eta ^3+O(|\zeta |^5),`$ (C4)
where

$`K_{30}^{10}`$ $`=`$ $`b^2I_{1,0,1}+abI_{1,1,1}+abI_{1,2,1}+a^2I_{1,3,1},`$
$`K_{21}^{10}`$ $`=`$ $`cq[3bI_{1,0,1}+(b2a)I_{1,1,1}+(2ba)I_{1,2,1}+3aI_{1,3,1}],`$
$`K_{12}^{10}`$ $`=`$ $`q^2[(ab+2c^2)I_{1,0,1}+(a^22c^2)I_{1,1,1}+(b^22c^2)I_{1,2,1}+(ab+2c^2)I_{1,3,1}],`$
$`K_{03}^{10}`$ $`=`$ $`cq^3[aI_{1,0,1}+aI_{1,1,1}bI_{1,2,1}+bI_{1,3,1}],`$
$`K_{03}^{01}`$ $`=`$ $`cq[bI_{1,0,1}+bI_{1,1,1}aI_{1,2,1}+aI_{1,3,1}],`$
$`K_{21}^{01}`$ $`=`$ $`q^2[(ab+2c^2)I_{1,0,1}+(b^22c^2)I_{1,1,1}+(a^22c^2)I_{1,2,1}+(ab+2c^2)I_{1,3,1}],`$
$`K_{12}^{10}`$ $`=`$ $`cq^3[3aI_{1,0,1}+(a2b)I_{1,1,1}+(2ab)I_{1,2,1}+3bI_{1,3,1}],`$
$`K_{30}^{10}`$ $`=`$ $`q^4[a^2I_{1,0,1}+abI_{1,1,1}+abI_{1,2,1}+b^2I_{1,3,1}],`$ (C5)
and the constants $`a`$, $`b`$, and $`c`$ are

$$a=\frac{Bq^2}{2}\left[1+\frac{\mathrm{Re}Q}{2|Q|}\right],b=\frac{Bq^2}{2}\left[1\frac{\mathrm{Re}Q}{2|Q|}\right],c=\frac{\gamma _2Bq^3}{2|Q|}.$$
(C6)
Finally, let us caution against using eq. (C1) for numerical calculations for all but small values of $`M`$ and $`N`$, and for $`z`$ close to 1, because it can become a difference between two nearly equal quantities. The hypergeometric series, which converges rapidly for $`z`$ close to 1, should be used instead.

## Appendix D Convergence of Series and Analytic Continuation

To prove that the series of equation (68) converges for sufficiently small $`|\zeta |`$, we need a bound on the coefficients $`I_{j,2m}`$. The cases $`B<2`$ and $`B>2`$ require different treatments. For $`B<2`$, we use the integral (66), replace the powers of $`(w+t_1)`$ in the denominator with the smaller $`(w+t_2)`$, and derive the bound

$$\left|I_{j,2m}\right|<\frac{1}{\pi }_0^{\mathrm{}}\frac{w^{jB1}dw}{(w+t_2)^{2j+2}}=t_2^{(B2)j2}_0^{\mathrm{}}\frac{v^{jB1}dv}{(v+1)^{2j+2}}.$$
(D1)
The simpler integral can be evaluated in terms of Gamma functions (AS, eq. 6.2.1). As this bound is independent of $`m`$, we can take it out of the inner summation over $`m`$ when we use it in the infinite series in equation (68). The inner sum can then be carried out to give simply $`(|\lambda |^2+|\nu |^2)^j=|\zeta |^{2j}`$. The series for the magnification is therefore majorized by the series

$$\frac{Bq^2}{\pi }\underset{j=1}{\overset{\mathrm{}}{}}\frac{\mathrm{\Gamma }[jB]\mathrm{\Gamma }[2+(2B)j]}{(2j+1)!}t_2^{j(B2)2}|\zeta |^{2j}.$$
(D2)
Applying the ratio test and using Stirling’s formula (AS, eq. 6.1.39) to approximate the Gamma functions, we find that the series (68) converges if

$$|\zeta |<2\left(\frac{t_2}{2B}\right)^{(2B)/2}B^{B/2}.$$
(D3)
This condition is met uniformly for all $`B`$ in $`[1,2]`$ when

$$|\zeta |<\sqrt{t_2},$$
(D4)
and then the infinite series (68) for $`(B)`$ is absolutely and uniformly convergent.

The integral (66) is not valid for $`B>2`$, and we must then use equation (67) and find a bound for the hypergeometric function. We can use the hypergeometric series (AS, eq. 15.1.1) for that function because its fourth argument lies in $`(0,1)`$. As all the other arguments are positive, all the terms of the infinite series are positive. Moreover, each of them is less than the corresponding term of the infinite series for $`{}_{2}{}^{}F_{1}^{}(2j+2,jB;2j+2;1t_2/t_1)`$ because $`\mathrm{}2j`$. This majorizing series is geometric (AS, eq. 15.1.8) and sums to $`(t_2/t_1)^{jB}`$. The terms $`(jBs)`$ in the product are now all positive for large enough $`j`$ because $`s2j+1`$, so we ignore the $`s`$ terms and obtain, for large $`j`$, the bound

$$\left|I_{j,2m}\right|<\frac{(jB)^{2j+1}}{(2j+1)!}t_2^{2m2j1}t_1^{jB2m1}.$$
(D5)
The series for the magnification is now majorized by the series

$$\frac{Bq^2}{t_1t_2}\underset{j=1}{\overset{\mathrm{}}{}}\frac{(jB)^{2j+1}t_1^{jB}}{(2j+1)!}\underset{m=0}{\overset{\mathrm{}}{}}\left(\genfrac{}{}{0pt}{}{j}{m}\right)\left(\frac{|\lambda |}{t_1}\right)^{2m}\left(\frac{|\nu |}{t_2}\right)^{2m2j}.$$
(D6)
The inner sum is

$$\left[\frac{|\lambda |^2}{t_1^2}+\frac{|\nu |^2}{t_2^2}\right]^j,$$
(D7)
and its maximum possible value is $`|\zeta |^{2j}/t_2^{2j}`$ because $`|\lambda |^2+|\nu |^2=|\zeta |^2`$. Hence, the infinite series (68) for the magnification invariant is majorized for $`B>2`$ by the series

$$\frac{Bq^2}{t_1t_2}\underset{j=1}{\overset{\mathrm{}}{}}\frac{(jB)^{2j+1}t_1^{jB}}{(2j+1)!}\left(\frac{|\zeta |^2}{t_2^2}\right)^j.$$
(D8)
The ratio test show that this series is convergent, and that the infinite series for the magnification is therefore absolutely convergent, if

$$|\zeta |<\frac{2t_2}{eBt_1^{B/2}}.$$
(D9)
Thus, however large $`B`$ is, the series for the magnification invariant converges for sufficiently small offset. Actually, if $`t_1`$, the larger root of the equation $`P^2=|Q|^2`$, is less than unity (as it is if the on-axis shear $`\gamma _1`$ is negative) and

$$\gamma _1>\frac{(\gamma _1^2+\gamma _2^2)q^2}{1q^2},$$
(D10)
then

$$Bt_1^{B/2}\left[\frac{e}{2}\mathrm{ln}\left(\frac{1}{t_1}\right)\right]^1$$
(D11)
for all $`B`$, and the magnification series (68) has a finite radius of convergence that is independent of $`B`$. When $`t_1=1`$ (as it is when there is no shear), then the bound (D9) decreases as $`1/B`$.

Regardless, if we need the magnification invariant $`(B)`$ for a model for which $`B=B_{\mathrm{max}}>2`$, then we can find a $`\delta `$, independent of $`B`$, such that the series for $`(B)`$ converges absolutely for $`|\zeta |<\delta `$ for $`1BB_{\mathrm{max}}`$. That is, our series is uniformly convergent as a function of $`B`$. The explicit expressions that we have for each term in the series (68), from equation (67) and equation (C1) of Appendix C, shows that each is composed of polynomials and exponential functions of $`B`$. Hence, each term of the uniformly convergent series is an analytic function of $`B`$. It follows, by a well-known theorem of complex analysis (e.g Levinson & Redheffer 1970), that the sum $`(B)`$ of the series is also an analytic function of $`B`$ for $`|\zeta |<\delta `$. It can therefore be continued beyond the range $`1<B<2`$ for which it was originally derived by the analysis of §4.1 to apply to models with larger $`B`$.

With only minor changes, the series for the configuration invariants can be handled in a similar way with a similar outcome; the infinite series (76) for $`S_{m,n}(B)`$ is absolutely convergent when the same inequalities (D4) and (D9) as for $`(B)`$ hold. Hence, the series for the configuration invariants of the four bright images also remain valid for $`B>2`$ for sufficiently small offset $`|\zeta |`$.

Formulas for one of the four bright images become singular at the radial caustic where it merges with the weak fifth image. Therefore we can not expect our series to converge there, and inequality (D9) does not hold at the known analytical $`|\zeta |`$ values given in Appendix A for special points on the radial caustic. However, the tangential caustic need not limit the convergence, even though series lose physical significance there, because the singularities of the two merging images can cancel. One can verify that the $`B<2`$ inequality (D3) may indeed hold at the cusps of the tangential caustic.


format: text

text:
# A first order phase transition induced by magnetic field and temperature

## Abstract

Taking the pseudobinary C15-Laves phase compound Ce(Fe0.96Al0.04)2 as a paradigm for studying a ferromagnetic(FM) to antiferromagnetic(AFM) phase transition, we present interesting thermomagnetic history effects in magnetotransport measurements across this FM-AFM transition. We argue that these distinctive hysteretic features can be used to identify the exact nature -first order or second order - of this kind of transition in magnetic systems where electrical transport is strongly correlated with the underlying magnetic order. A comparison is made with the similar FM-AFM transitions observed in Nd and Pr-based manganese compounds with perovskite-type structure.

The nature of ferromagnetic (FM) to antiferromagentic (AFM) transition in the perovskite-type manganese oxide compounds Nd1/2Sr1/2MnO3 and Pr1/2Sr1/2MnO3 has been the subject of close scrutiny in recent years . The FM-AFM transition observed in these compounds is taken as a sort of a prototype of a first order transition, and certain thermomagnetic features have been highlighted which are thought to be generic of a first order phase transition . In the same spirit we have undertaken a study of FM-AFM phase transition in the pseudobinary C15-Laves phase compound Ce(Fe0.96Al0.04)2. We find striking thermomagnetic history effects in magnetotransport measurements across the FM-AFM transition in this interesting system.

CeFe2, with its relatively low Curie temperature ( T$`{}_{C}{}^{}`$230K) and reduced magnetic moment ($``$ 2.3$`\mu _B/f.u.`$) , is on the verge of a magnetic instability . Neutron measurement has shown the presence of antiferromagnetic fluctuations in the FM ordered state of CeFe2 below 100K . With small but suitable change in electronic structure caused by doping with elements like Co, Al, Ru, Ir, Os and Re at the Fe-site of CeFe2 , these antiferromagnetic fluctuations get stabilized into a low temperature AFM state, and after certain concentration of dopants (usually 5 to 10%) this AFM phase replaces the FM phase altogether .

While most recent experimental efforts are mainly focussed on understanding the cause of this magnetic instability  in CeFe2, there exists one other aspect of the observed magnetic properties which needs proper attention, viz. the exact nature of the FM-AFM transition. We have recently addressed this second question in Ru and Ir-doped CeFe2 alloys . In this paper we shall focus on the Al-doped CeFe2 alloys. In contrast to the Ru, Co and Ir doped alloys, a distinct co-existence of FM and AFM phase has been reported around the FM-AFM transition temperature in Al-doped CeFe2 alloys . We report here interesting thermomagnetic history dependence of magnetotransport in a Ce(Fe0.96Al0.04)2 alloy. We argue that these thermomagnetic history effects arise due to the first order nature of the FM-AFM transition. These effects are broader manifestations of the behavior reported earlier in manganese compounds , and can be used to identify a first order FM-AFM transition in a new system, especially in those where the electrical transport is strongly correlated with the underlying magnetic order.

The details of the preparation and characterization of the sample can be found in Ref.10. The samples from the same batch have been used earlier in the study of bulk magnetic and transport properties , and neutron measurements . We have used a superconducting magnet and cryostat system (Oxford Instruments, UK) for magnetotransport measurements as a function of temperature (T) and applied magnetic field (H). The resistivity is measured using a standard dc-four probe technique.

The inset of Fig. 1(a) shows the magnetization (M) vs T plot for the Ce(Fe0.96Al0.04)2 sample measured with an applied field of 2 mT. The sample undergoes a paramagnetic (PM) to FM transition at around 195K, followed by a lower temperature FM-AFM transition around 90K. These results are in consonance with the earlier bulk properties  and neutron measurements . We shall now study the field dependence of magnetoresistance in various temperature regimes. In the main panel of Fig.1 we present resistivity ($`\rho `$) as a function of H at T=3K, 5K, 20K, obtained after initial zero-field-cooling (ZFC) the sample to the temperature concerned. The $`\rho `$ vs H plot at T$``$120K (not shown here) is that of a typical ferromagnet, showing clear negative magnetoresistance. In the antiferromagnetic regime (see Fig. 1(a)-(c)), we see the clear signature of a field induced ferromagnetic transition at a field HM, where the resistivity decreases sharply with the increase in H. (The slight increase in $`\rho `$ in the field regime H$``$HM, indicating the positive magnetoresistance of the AFM state is not quite visible in the same scale). Although the change in resistivity due to this field induced transition is not as drastic as in Nd1/2Sr1/2MnO3 , this is within an order of magnitude of those obtained with the similar applied fields in Pr1/2Sr1/2MnO3 . It is to be noted here that the FM-AFM transition in the present sample is not accompanied by a metal-insulator transition and both the FM and AFM states remain metallic. On reducing the field from well above HM, a distinct hysteresis is observed in the $`\rho `$ vs H plot (see Fig. 1(a)-(c)). We attribute this to the first order nature of the field induced AFM-FM transition. While reducing H from well inside the FM state i.e. H$`>>H_M`$, the FM state continues to exist as supercooled metastable state below HM up to a certain metastability field H . Between HM and H fluctuations will help in the formation of droplets of the stable AFM state, and at H an infinitisimal fluctuation will drive the whole system to the stable AFM state. Similar hysteresis in the $`\rho `$-H plots of Nd1/2Sr1/2MnO3 and Pr1/2Sr1/2MnO3 has also been attributed to the first order nature of the phase transition . The role of thermal fluctuations is expected to be reduced in the very low temperature regime, and this is clearly seen in the $`\rho `$-H plots at T=3K and 5K (see Fig. 1(b) and (c)). In this temperature regime on reduction of the applied H to zero the $`\rho `$(H=0) lies distinctly below the intial ZFC-$`\rho `$(H=0), thus giving rise to an open hysteresis loop. This kind of open hysteresis loop has earlier been reported for Nd1/2Sr1/2MnO3  but not for Pr1/2Sr1/2MnO3 . We attribute this behaviour to the existence of a residual metastable FM state even when the applied H is reduced to zero. On increasing H on the negative side, the $`\rho `$(H) curve is clearly not symmetric to the virgin $`\rho `$(H) curve on the positive H side (see Fig. 1(b)-(c)). The field induced AFM-FM transition, however, takes place at the same $`|H_M|`$ (see Fig. 1(b)-(c)). A distinct hysteresis is observed on reducing H from the same $`|H_{max}|`$ as in the positive side. The envelope hysteresis loop now closes at H=0, i.e. it merges with the starting value where from the field excursion on the negative direction had started. Further on increasing H to the positive side, the $`\rho `$(H) curve now follows a path which is distinctly below the virgin $`\rho `$(H) curve but very symmetric to the $`\rho `$(H) curve on the negative side in the increasing H cycle. This $`\rho `$(H) curve (henceforth will be termed as forward envelope curve) merges with the virgin curve in the H-regime well beyond HM. We have checked this distinct difference between the virgin curve and the forward envelope curve in the negative H side also, by drawing a virgin curve in the negative H direction after zero field cooling the sample (see dashed line in Fig. 1(b)-(c)). This anomalous behaviour of virgin curve lying distinctly outside the envelope curve is observed  in the field dependence of magnetization as well (see inset of Fig. 1(b)).

From a closer inspection of the published results on the Nd1/2Sr1/2MnO3 sample , we expect similar behaviour to take place in that sample in the temperature regime T$``$20K. Specifically the $`\rho `$-H curve was shown to have an open hysteresis-loop (see figs. 2 B to 2D of ref.1). We expect that if H was reduced to -12T in that sample, and then raised back to +12T, this second leg would be closed similar to our Fig. 1(a). Similarly an isothermal M-H measurement in that compound should show a virgin curve lying outside the envelope curve as in the inset of Fig. 1(b).

Another interesting aspect worth noting in our present sample is that the slope of the $`\rho `$(H) return envelope curves change sharply on crossing H=0 in either direction (see Fig. 1(b)-(c)). It almost flattens in the low field regime H$`<|H_M|`$ on either side of H=0, as if the domains of the residual FM state still retain their previous memory. With the striking similarity of the field dependence of resistivity in Nd1/2Sr1/2MnO3 (see Fig. 2(b)-(d) of Ref.1) with those in the positive H cycle in the present sample (Fig. 1(b)-(c)), it is quite tempting to predict the similar sharp change in slope of R(H), as H changes sign of H in Nd1/2Sr1/2MnO3 as well.

Supercooling/superheating and metastability have been identified as key elements to explain the thermomagnetic history effects associated with the first order FM-AFM transition in perovskite-type manganese oxide systems . To explain the strong temperature dependence of this thermomagnetic irreversibility in our present system, we now invoke in addition the concept of the limit of the metastability H(T) ( or T(H)). The implicit assumption involved here is that the difference between the phase transition line HM(T) ( or TN(H)) and the limit of metastability H(T) ( or T(H) ) widens with the decrease (increase) in temperature (magnetic field). To support this assumption we shall now study the temperature dependence of resistivity in the presence of various applied magnetic fields.

In Fig. 2 we present $`\rho `$ vs T plots with H=0, 0.5T, 2T and 3T. Appearance of magnetic superzones  at the FM-AFM transition (TN) give rise to a distinct structure in the form of a local minimum in the $`\rho `$(T) (see Fig.2). There is a marked hysteresis associated with this transition, emphasizing the first order nature of the transition. We argue that the FM (AFM) phase exists as supercooled (superheated) metastable phase in the cooling (heating) cycle in this hysteretic temperature regime. In the high (low) temperature reversible regime the only magnetic phase is the stable FM (AFM) phase. In an earlier zero field neutron measurement  a clear co-existence of FM-AFM phase was observed for a substantial temperature regime below the onset of phase transtion. This alongwith the gradual development of the cubic to rhombohedral structual distortion  as the temperature is lowered, fit naturally in our present picture of first order phase transition. In the presence of finite H  the TN is suppressed and the hysteresis is enhanced substantially, so much so that with H=2T no reversible regime is observed below TN down to 5K. This in turn implies that the limit of metastability T(H) (below which one should see the reversible response of the stable AFM phase) gets suppressed even faster in comparison to TN(H), and the metastable regime (encompassed between TN(H) and T(H)) widens with the decrease in T or with increase in H. This is again in perfect consonance with our observation in the isothermal field dependence of resistivity.

Summarizing our results, the important findings of the present study are the following:

1. The butterfly $`\rho `$(H) hysteresis loop with anomalous virgin curve, which does not close at H=0 in the low temperature regime. This is complimented with magnetization study as well.
2. The distinct hysteresis in the $`\rho `$(T) curve at AFM-FM transition, which gets enhanced in the presence of applied magnetic field.

These observations along with the earlier neutron measurements  help to establish the first order nature of the AFM-FM transition in the CeFe2-based pseudobinary systems. This information in turn will be important for developing a theoretical model to explain the interesting electro-magnetic properties of CeFe2, which does not exist so far. Most importantly, thermomagnetic history effects observed here can be used as generic features to identify a first order FM-AFM phase transition. In fact the double hysteresis or butterfly loop in polarization measurements is regularly used to identify a first order transition in the ferroelectric/antiferroelectric materials . (However, we are unaware of any report regarding the anomalous nature of the virgin curve associated with the butterfly loops in the ferroelectric materials). Such a field induced first order transition can be explained in terms of free energy curves obtained by expanding in a power series in polarization and retaining only terms with even powers up to sixth order . Phase coexistence and metastability across the first order phase transition, and hence hysteresis arise naturally out of such free energy curves. Although, sharp rise in field induced magnetization (and associated hysteresis)  and sharp drop in the field induced resistivity  in magnetic systems are regularly taken as signatures of AFM-FM transition, no clearcut inference can be made from such observations regarding the exact nature of the phase transition. Our present study along with those on perovskite-type manganese oxide systems  attempt to fill-up this gap and provide some relatively easy means to identify the nature of such magnetic transitions from standard bulk properties measurements. It should be noted here that the resistivity/magnetoresistance is not a thermodynamic quantity, whereas magnetization/staggered magnetization is, and can be used as an order parameter in the Landau type free energy expansion. However, intricate correlation between the resistivity/magnetoresistance and magnetization/staggered magnetization across the FM-AFM transition, makes these magnetotransport properties suitable observables in the study of the nature of the phase transition. Further, these samples with well characterized electro-magnetic properties can be used as paradigms to study the various interesting aspects of a first order transition, namely nucleation and growth, and path dependence of the transition, in a relatively easy and reproducible manner. Such studies are not that easy and reproducible in the more common cases of first order transition, like melting (solidification) of solids (liquids) or vaporization (condensation) of liquids (gases).