format: text

text:
# References

Is it possible to clone using an *arbitrary* blank state?

Anirban Roya111res9708@isical.ac.in , Aditi Sen(De)b222aditisendein@yahoo.co.in  and Ujjwal Senb333ujjwalsen@yahoo.co.in

a*Physics and Applied Mathematics Unit, Indian Statistical Institute, 203 BT Road, Kolkata 700035, India*

b*Department of Physics, Bose Institute, 93/1 APC Road, Kolkata 700009, India*

## Abstract

We show that in a cloning process, whether deterministic inexact or probabilistic exact, one can take an arbitrary blank state while still using a fixed cloning machine.

Quantum information cannot be cloned. There cannot exist a machine which can produce two (or more) exact copies of an arbitrary state in a deterministic manner . Replication is not allowed even when the input state is taken at random from a given set of two non-orthogonal states . It has been further shown that probabilistic exact cloning is not possible if the input state is from a given linearly dependent set .

Although exact cloning is not possible, one can approximately clone an arbitrary input state . In this scheme a *fixed* state is taken as the blank state, depending on which and the initial machine state, the cloning operation is constructed. We extend this operation such that any *arbitrary* (pure or mixed) state, taken as the blank copy, can do the job.

The optimal universal $`12`$ inexact qubit cloner of Bruß*et al.* takes an arbitrary qubit $`|\psi \psi |=\frac{1}{2}(I+\stackrel{}{s}.\stackrel{}{\sigma })`$ along with a fixed blank state $`|b`$ and a machine state $`|M`$ as input. An entangled state of the three qubits is produced as the output such that the reduced density matrices of the first two qubits are two similar approximate copies $`\rho =\frac{1}{2}(I+\eta \stackrel{}{s}.\stackrel{}{\sigma })`$ of $`|\psi \psi |`$ with $`\eta =\frac{2}{3}`$. The unitary operator realizing this process is defined on the combined Hilbert space of the input qubit, the blank qubit and machine by

$$U^{}|0|b|M=\sqrt{\frac{2}{3}}|00|m+\sqrt{\frac{1}{6}}(|01+|10)|m_{}$$

$$U^{}|1|b|M=\sqrt{\frac{2}{3}}|11|m_{}+\sqrt{\frac{1}{6}}(|01+|10)|m$$
(1)
where $`|b`$ is a *fixed* blank state (in a two-dimensional Hilbert space), $`|M`$ is the initial state of the machine, $`|m`$ and $`|m_{}`$ being two mutually orthonormal states of the machine Hilbert space. The two clones are to surface at the first and second qubits. Note that the machine has turned out to be a qubit.

As it stands, the unitary operator $`U^{}`$ depends on the blank state $`|b`$ and the machine state $`|M`$. And the quality of the clones could be badly affected if $`|b`$ gets changed to an unknown state, say by some environment-induced decoherence. We show that by suitably constraining the unitary operator it is possible to keep the clones intact, even in this changed scenario. After dealing with the $`12`$ qubit cloner, we show that the same is true for the most general cloner, the $`NM`$ *qudit* (elements of a $`d`$-dimensional Hilbert space) cloner. We carry over these considerations to the case of probabilistic exact cloning.

Let us suppose that for the $`12`$ qubit cloner, the machine state $`|M`$ belongs to a four-dimensional Hilbert space. And let the unitary operator $`U`$ be defined on the combined Hilbert space of the input qubit, blank qubit and the four-dimensional Hilbert space of the machine by

$$U|0|b|M=\sqrt{\frac{2}{3}}|00|M_0+\sqrt{\frac{1}{6}}(|01+|10)|M_1$$

$$U|1|b|M=\sqrt{\frac{2}{3}}|11|M_1+\sqrt{\frac{1}{6}}(|01+|10)|M_0$$

$$U|0|b_{}|M=\sqrt{\frac{2}{3}}|00|M_2+\sqrt{\frac{1}{6}}(|01+|10)|M_3$$

$$U|1|b_{}|M=\sqrt{\frac{2}{3}}|11|M_3+\sqrt{\frac{1}{6}}(|01+|10)|M_2$$

where $`M_i|M_j=\delta _{ij}`$ ($`i,j=0,1,2,3)`$ and $`b|b_{}=0`$.

Let $`|B=c|b+d|b_{}`$ be an arbitrary pure state of the Hilbert space of the blank qubit. Then

$$U|0|B|M=\sqrt{\frac{2}{3}}|00|X+\sqrt{\frac{1}{6}}(|01+|10)|X^{}$$

$$U|1|B|M=\sqrt{\frac{2}{3}}|11|X^{}+\sqrt{\frac{1}{6}}(|01+|10)|X$$
(2)
where

$$|X=c|M_0+d|M_2$$

$$|X^{}=c|M_1+d|M_3$$

are orthogonal. This form is exactly the same as in equation (1). Thus an arbitrary input qubit $`|\psi `$ would be just as well cloned by equation (2) as it would be through equation (1).

We now consider the most general cloning machine, the one that produces $`M`$ approximate copies of the input, when $`N(<M)`$ $`d`$-dimensional inputs are provided . The corresponding unitary operator need only be defined on the symmetric subspace444The symmetric subspace is defined by the linear span of the set of all tensor product states $`|\psi |\psi \mathrm{}N`$ times, $`|\psi `$ being any $`d`$-dimensional state.  of the $`d^N`$-dimensional Hilbert space of the $`N`$ input qudits. It is defined by

$$U_{NM}^{}|\stackrel{}{n}|R=\underset{\stackrel{}{j}=0}{\overset{MN}{}}\alpha _{\stackrel{}{n}\stackrel{}{j}}|\stackrel{}{n}+\stackrel{}{j}|M_\stackrel{}{j}$$
(3)
where $`\stackrel{}{n}=(n_1,n_2,\mathrm{}.,n_d)`$, $`|\stackrel{}{n}`$ is a completely symmetric and normalised state with $`n_i`$ systems in $`|i`$ with $`_{i=1}^dn_i=N`$, $`\stackrel{}{j}=(j_1,j_2,\mathrm{}.,j_d)`$ with $`_{k=1}^dj_k=MN`$, $`|R`$ denoting the combined state of the $`MN`$ fixed pure $`d`$-dimensional blank states $`|b_d`$ and the initial state $`|M`$ of the cloning machine and $`|M_\stackrel{}{j}`$ denoting the orthonormal states of the cloning machine. And

$$\alpha _{\stackrel{}{n}\stackrel{}{j}}=\sqrt{\frac{(MN)!(N+d1)!}{(M+d1)!}}\sqrt{\underset{k=1}{\overset{d}{}}\frac{(n_k+j_k)!}{n_k!j_k!}}$$

There are $`s`$ equations required to define $`U_{NM}^{}`$, where $`s`$ is the dimension of the symmetric subspace . The required dimension of the machine is $`D=\frac{(MN+d1)!}{(MN)!(d1)!}`$.

Now we proceed as we had done for the $`12`$ qubit cloner. If we want to allow an arbitrary pure state (possibly entangled) of the $`d^{MN}`$-dimensional Hilbert space of the blank states to act as the new blank state and still produce the same outputs, we have to use a $`(D\times d^{MN})`$-dimensional machine. The new unitary operator $`U_{NM}`$ satisfies, along with the $`s`$ equations in (3), $`(d^{MN}1)s`$ more equations corresponding to the $`d^{MN}1`$ more blank states on which the new operator is to be defined. Then by linearity, for an arbitrary pure blank state $`|B_1`$, we would have

$$U_{NM}|\stackrel{}{n}|B|M=\underset{\stackrel{}{j}=0}{\overset{MN}{}}\alpha _{\stackrel{}{n}\stackrel{}{j}}|\stackrel{}{n}+\stackrel{}{j}|X_\stackrel{}{j}$$
(4)
where $`X_\stackrel{}{j}`$ are orthonormal states of the machine. This has the same form as eq. (3) and thus would be equally efficient in producing the requisite approximate copies.

Similar considerations carry over to the case of probabilistic exact cloning. Although we only consider the case qubits, the considerations carry over to higher dimensions. In the case of qubits, instead of the $`U_1^{}`$ defined by

$$U_1^{}|\psi _0|b|M=\sqrt{\gamma }|\psi _0|\psi _0|m+\sqrt{1\gamma }|\mathrm{\Phi }$$

$$U_1^{}|\psi _1|b|M=\sqrt{\gamma }|\psi _1|\psi _1|m+\sqrt{1\gamma }|\mathrm{\Phi }$$

with $`\gamma =1/\left(1+\left|\psi _0|\psi _1\right|\right)`$ ($`|m`$ and $`|\mathrm{\Phi }`$ are orthogonal), $`|\psi _0`$, $`|\psi _1`$ being two non-orthogonal states which are to be probabilistically cloned, we define $`U_1`$ as

$$U_1|\psi _0|b|M=\sqrt{\gamma }|\psi _0|\psi _0|M_0+\sqrt{1\gamma }|\mathrm{\Phi }$$

$$U_1|\psi _1|b|M=\sqrt{\gamma }|\psi _1|\psi _1|M_0+\sqrt{1\gamma }|\mathrm{\Phi }$$

$$U_1|\psi _0|b_{}|M=\sqrt{\gamma }|\psi _0|\psi _0|M_1+\sqrt{1\gamma }|\mathrm{\Phi }^{}$$

$$U_1|\psi _1|b_{}|M=\sqrt{\gamma }|\psi _1|\psi _1|M_1+\sqrt{1\gamma }|\mathrm{\Phi }^{}$$

where $`|M_0`$, $`|M_1`$, $`|\mathrm{\Phi }`$, $`|\mathrm{\Phi }^{}`$ are mutually orthogonal. Then

$$U_1|\psi _0|B|M=\sqrt{\gamma }|\psi _0|\psi _0|m^{}+\sqrt{1\gamma }|\mathrm{\Phi }^{\prime \prime }$$

$$U_1|\psi _1|B|M=\sqrt{\gamma }|\psi _1|\psi _1|m^{}+\sqrt{1\gamma }|\mathrm{\Phi }^{\prime \prime }$$

for an arbitrary blank qubit $`|B=c|b+d|b_{}`$ so that $`|m^{}=c|M_0+d|M_1`$ and $`|\mathrm{\Phi }^{\prime \prime }=c|\mathrm{\Phi }+d|\mathrm{\Phi }^{}`$ are orthogonal states. Consequently, the probabilistic cloning goes through with the same optimal efficiency even if we use an arbitrary blank pure qubit.

As we have mentioned, the main motivation behind consideration of an arbitrary blank state was decoherence. Decoherence, however, usually produces *mixed* states. But in our discussion we have only considered pure states. We now show that the input blank state could as well be an arbitrary mixed state. For definiteness, let us consider only mixed qubits. Any such mixed state can be written as $`\rho =a_1|a_1a_1\left|+a_2\right|a_2a_2|`$ where $`a_1,a_20`$, $`a_1+a_2=1`$ and $`a_1|a_2=0`$. Since $`|a_1`$ and $`|a_2`$ are pure states, by linearity it follows from (2) that the same unitary operator that allowed an arbitrary pure blank state, would just as well clone even when the blank state is a mixed qubit.

To summarize, we have shown that an unknown blank state can be used for cloning, whether it is deterministic inexact or probabilistic exact.

We thank the anonymous referee for making useful suggestions to revise our earlier manuscript. We acknowledge Guruprasad Kar, Sibasish Ghosh, Debasis Sarkar, Pinaki Pal and Mridula Kanoria for helpful discussions. A.S. and U.S. thanks Dipankar Home for encouragement and U.S. acknowledges partial support by the Council of Scientific and Industrial Research, Government of India, New Delhi.


text:
# The ordering temperature and critical exponents of the binomial Ising spin glass in dimension three

## I Introduction

The Edwards-Anderson Ising spin glass (ISG) has been studied intensively for over twenty five years. In this model Ising spins on a regular lattice in dimension three are coupled by random near neighbour interactions EA . Despite its deceptively simple definition, the basic physics of this canonical complex system still eludes consensus. Recently interest has focused on the correct description of the ground state and the low energy excited states ground , but finite temperature aspects are also highly important, in particular the critical behaviour at the freezing temperature $`T_g`$. Over the years a number of numerical efforts have been consecrated to accurate measurements near $`T_g`$, with the aim first of demonstrating that $`T_g`$ is non-zero, and then of determining the critical exponents precisely BY ; og ; KY ; pala ; pom ; Bal . The long relaxation times and the need to average over large numbers of independent samples renders the task numerically laborious. In addition, it has not always been appreciated that corrections to finite size scaling can be singularly vicious in this system.

Here we will concentrate on the problem of the accurate and reliable determination of $`T_g`$ and the critical exponents for the ISG with binomial ($`\pm J`$) near neighbour interactions in a dimension three simple cubic lattice ($`L\times L\times L`$),

$$H=J_{i,j}S_iS_j$$
(1)
with $`J_{i,j}`$ chosen at random as $`\pm 1`$. As this system is the dropsophilia of spin glass studies, it is important to establish as precisely and reliably as possible the basic parameters of its critical behaviour. We will show that a parameter set can be obtained consistent with most published data from different complementary approaches, if corrections to finite size scaling are properly taken into account.

We will principally refer to the simulations of references og ; KY ; pala ; pom ; Bal ; the estimates drawn from the analyses of the numerical data are shown in Table I. We will comment on these data as we go along. The main numerical contribution of the present work is to present non-equilibrium data of higher statistical accuracy than before, and to allow for corrections to finite size scaling in the analysis of these results also.

As a general rule, the apparent exponent values are strongly correlated with the value assumed for the ordering temperature, so once $`T_g`$ is well determined it is relatively easy to obtain good values for the exponents. On the contrary, if the wrong $`T_g`$ is used, the estimates of the exponents will be biased. Finite size scaling methods have been widely used, but corrections to finite size scaling have often been overlooked, which is dangerous. In Table I the estimates of the ordering temperatures and of the exponents from different publications should not be considered as independent of each other, but each parameter set should be taken as a whole.

The autocorrelation function at time $`t`$ is defined by

$$q(t)=<S_i(0).S_i(t)>$$
(2)
where the average ($`<..>`$)is over all spins. In practice a large number of further averages (denoted by $`[..]`$)are taken over many samples all of size $`L^3`$ spins with microscopically different arrangements of the bonds. We will use $`q`$ to mean $`q(t)`$ at times long enough for the system to have entered the equilibrium fluctuation regime. Periodic (toroidal) boundary conditions are generally but not always used (see Bal ). The equilibrium spin glass susceptibility is given by

$$\chi _L(T)=L^3[<q^2>]$$
(3)
where the average is taken over time once the sample has been annealed to equilibrium either by carrying out enough annealing steps or by using the more sophisticated and efficient parallel tempering method KY ; Bal . The non-equilibrium susceptibility $`\chi (t)`$ is the spin glass susceptibility at time $`t`$

$$\chi (t)=L^3[<q^2(t)>]$$
(4)
after a simple anneal which starts from a totally uncorrelated (infinite temperature) state and lasts a time $`t`$ huse .

## II Binder parameter

The first large scale simulations on this problem were those of Ogielski og  who did massive dynamic and static measurements on samples as large as $`L=64`$. His values were considered authoritative until finite size scaling techniques were later intensively exploited BY ; KY ; pala ; Bal . The standard scaling method to estimate $`T_g`$ from finite size data in ISGs is through the use of the dimensionless Binder parameter,

$$g_L(T)=[3<q^4>/<q^2>^2]/2$$
(5)
where $`<q^2>`$ and $`<q^4>`$ are moments of the fluctuations of the autocorrelation parameter $`q(t)`$ in thermal equilibrium at temperature $`T`$. The averages are taken first over time and then over samples. In absence of corrections to finite size scaling, the family of curves $`g_L(T)`$ for different sizes $`L`$ must all intersect precisely at $`g_0`$ for $`T=T_g`$. For the 3d binomial ISG, the $`g_L(T)`$ curves come together near $`T=1.20`$ BY  but very good statistics are needed to see intersections KY ; Bal . Precisely because the $`g_L(T)`$ curves lie very close together, even small corrections to finite size scaling can have a drastic effect on the crossing point temperatures. As corrections to finite size scaling are already visible for the spin glass susceptibility pala ; pom , corrections can be expected a fortiori for the Binder parameter. If $`\omega `$ is the leading correction to scaling exponent, at $`T_g`$ the successive $`g_L`$ values for increasing $`L`$ can be expected to behave as

$$g_L=g_0[1+constL^\omega ]$$
(6)
It has been pointed out Bal  that in these circumstances it is helpful to use the ”quotient method”. An elementary version of this method is to use ”crossing points”. Successive intersection temperatures $`T^{}(L)`$, corresponding to intersections of curves $`g_L(T)`$ and $`g_{2L}(T)`$, are plotted. If the leading correction exponent $`\omega `$ dominates, then $`T^{}(L)`$ values should extrapolate to the true infinite $`L`$ ordering temperature $`T_g`$ as

$$T^{}(L)=T_g+constL^{(\omega +1/\nu )}$$
(7)
Extrapolation to infinite $`L`$ provides a precise estimate for $`T_g`$, if the corrections to finite size scaling are fully under control. In the case of the 3d binomial ISG, accurate $`g_L(T)`$ data for sizes $`L=5,10,20`$ Bal , for $`L=6,8,12,16,24`$ KY  and for $`L=3,4`$ pom  can be used to obtain intersections. It should be noted that for technical reasons the boundary conditions used by Ballesteros et al Bal  on a specially built dedicated computer are helicoidal while for the other calculations the boundary conditions are periodic. A direct comparison of the raw data from cruz  and KY  indicates that systematic differences in the results because of the different boundary conditions are negligible for the spin glass susceptibility at $`L=20`$. The sequence of Binder parameter curves for $`L=3,4,5,6,8`$ pom ; Bal ; KY  behaves smoothly, although the $`L=5`$ curve is the only one in the set to have helicoidal boundary conditions. A subset of the data from the different sources is shown for illustration in figure 1.

It can be seen that the behaviour of the $`g_L(T)`$ curves and hence of the crossing temperatures varies in a non-monotonic manner with $`L`$. (We suggest below that this may be because $`\omega `$ is high and analytic corrections may be non-negligible). The $`\omega `$ and $`\nu `$ values of Ballesteros et al Bal  correspond to an effective exponent $`(\omega +1/\nu )=1.3(4)`$. (Estimates of $`\omega `$ vary considerably from author to author and we will return to this point later). In figure 2 we have plotted $`T^{}(L)`$ in the form of equation (2) using the central value of the exponent. The pairs of sizes $`[2L,L]`$ we have used are  and  where the boundary conditions are helicoidal, ,, and  where the boundary conditions are periodic. There is no very obvious systematic difference between the two sets of points with different boundary conditions. Even for the larger sizes the crossing temperature is still increasing fairly rapidly with increasing $`L`$ ; the very high quality data of Bal  show a $`L=10,20`$ crossing already at $`T^{}=1.118(10)`$, so the central estimate given by KY , $`T_g=1.11`$, is definitely too low. The curve in figure 2 (a two paramereter spline fit) is just a schematic indication of one possible extrapolation to to infinite $`L`$ which would appear to lead to a significantly higher estimate for $`T_g`$ than that of Bal , $`T_g=1.138(10)`$. Alternative curves could be drawn through the data points, but the Binder parameter data (including the high quality results of Bal ) point to a higher $`T_g`$ than their central estimate. The non-monotonic behaviour of the crossing points means that not only the leading correction term is operative. Thus in the case of this particular ISG system the Binder parameter method can only give a rough estimate of $`T_g`$, because the correction terms are badly under control.

It should be noted that because the behaviour of the crossing temperature $`T^{}(L)`$ is not monotonic, if only measurements of $`g_L(T)`$ for $`L=4,8,`$ and $`16`$ were available one would see an excellent unique intersection of the three $`g_L(T)`$ curves, but this intersection point is NOT at the true infinite size limit temperature $`T_g`$. It is an artefact of the non-monotonic behaviour of $`T^{}(L)`$.

Binder parameter data can also be directly used to estimate the exponent $`\nu `$, as

$$d(g_L)/d(1/T)L^{1/\nu }$$
(8)
at $`T_g`$. Using a single $`g_L(T)`$ data set, the value of $`\nu `$ estimated in this way depends strongly on the $`T_g`$ estimate. The variation of the $`g_L(T)`$ slopes with $`L`$ in the $`g_L`$ data of Ballesteros et al provide an apparent $`\nu `$ value which is about $`2.15`$ if $`T_g=1.138`$ (their preferred value); if $`T_g`$ is higher then the associated estimate for $`\nu `$ from the same data set is lower. Once again estimates can be affected also by correction terms. (Similarly, the apparent $`\eta `$ from spin glass susceptibility data changes rapidly with the assumed $`T_g`$, Figure 7).

## III Correlation length

Ballesteros et al Bal  do not in fact base their $`T_g`$ estimate on their Binder parameter data, but use an independent criterion : the intersection of curves for the ratio of the size dependent correlation length to $`L`$. However they show data for three sizes only ($`L=5,10,20`$) and do not discuss possible corrections to finite size scaling for this parameter, even for $`L`$ as low as $`L=5`$. Even though their intersection point for these three sizes appears to be unique and very clear, the $`T_g`$ estimate from this criterion ($`T_g=1.138(10)`$) is lower than the prefered $`T_g`$ value obtained from other measurements which we will discuss below, $`T_g=1.195(15)`$. Although both series of measurements show clear evidence for a well established transition, and though the difference between the estimated transition temperatures is only about $`4\%`$, it would be more satisfactory if complete agreement could be reached.

## IV Series

It is of interest to note that quite independent information from series calculations can throw further light on the values of the ordering temperature and the exponents. Long series results on the 3d binomial system singh  provide sets of parameters $`[T_g,\gamma ]_a`$ corresponding to different approximants, Figure 3. The parameter pair $`[T_g,\gamma ]`$ representing the true phase transition can be expected to lie somewhere within the cloud of approximant points. In Figure 3 we represent the approximant points from singh  together with the equivalent points from the different numerical simulation estimates from Table I (using $`\gamma =\nu (2\eta )`$). It can be seen that while the approximant cloud has a relatively wide spread in temperature and so is not very discriminatory for $`T_g`$, the cloud has a fairly narrow distribution in $`\gamma `$. The $`\gamma `$ values corresponding to the parameter sets of references Bal , KY and pala  in Table I lie systematically above the series data, while the $`[T_g,\gamma ]`$ pairs of og  and pom  fall neatly in the centre of the approximant cloud. In dimension 4 the agreement between series singh  and simulations d4  is excellent, so there seems no a priori reason to expect the series results to be strongly biased in dimension 3.

## V Non-equilibrium methods

We have proposed a novel technique including non-equilibrium data for identifying $`T_g`$ numerically through the combination of three simulations lwb ; pom . We will use it again in the present context in a slightly modified form.

First, the non-equilibrium spin glass susceptibility for samples quenched to $`T_g`$ increases with time as $`t^h`$ where

$$h=(2\eta )/z$$
(9)
$`\eta `$ and $`z`$ being the standard equilibrium static and dynamic exponents huse .

Secondly the autocorrelation function decay $`q(t)`$ at $`T_g`$ is of the form $`t^x`$ , if a sample is initially annealed at $`T_g`$ during a long waiting time $`t_w`$, under the condition that the measuring time for the decay $`t`$ is much shorter than $`t_w`$ og ; rieger . The exponent $`x`$ at $`T_g`$ is related to the equilibrium exponents through $`x=(d2+\eta )/2z`$ og .

Thirdly, the equilibrium spin glass susceptibility at $`T_g`$ increases with sample size as BY

$$\chi _L/L^2L^\eta $$
(10)
The effective exponent $`h(T)`$ can be measured directly as a function of trial temperatures $`T`$. It can also be evaluated indirectly and independently through combining the effective exponent $`\eta (T)`$ from the finite size scaling of the equilibrium susceptibility, with the effective decay exponent $`x(T)`$. We will refer to this derived value as $`h^{}(T)`$ :

$$h^{}(T)=2x(T)(2\eta (T))/(d2+\eta (T))$$
(11)
Measurements of each of the three effective parameters ($`h(T),x(T),\eta (T)`$) can be made at a series of trial temperatures $`T`$; consistency dictates that at $`T_g`$ we must have $`h=h^{}`$ so when $`h(T)=h^{}(T)`$, then $`T=T_g`$.

Once again, care must be taken to avoid errors arising from corrections to finite size scaling. For the measurements of $`h(T)`$, large samples can readily be used, so it is relatively easy to insure that for the time range used the system is far from the saturation limit for that particular size, see zheng . At short times there can however be finite time corrections, the analogue of finite size corrections, as the clusters of correlated spins that are gradually building up with annealing time $`t`$ have a $`t`$ dependent finite size, even if the sample size $`L`$ can treated as effectively infinite. Including the finite time correction, we should write parisi

$$\chi (t)=At^h[1Bt^{w/z}]$$
(12)
with $`h(T)`$ as above, and $`w`$ the leading dynamic correction to scaling exponent (which does not have to be equal to the static exponent $`\omega `$ma ). The factors $`A`$ and $`B`$ may vary with $`T`$. An example of data is shown for two alternative plots in figures 4 and 5. (Rather than $`\chi (t)`$ we plot $`\chi (t)1`$ because at infinite $`T`$, $`\chi (t)=1`$. Measurements were averaged over runs on 6766 samples of size $`L=28`$ at temperature $`T=1.2`$. For this size saturation to the finite size limited equilibrium susceptibility would only occur at $`t10^8`$). For this particular temperature the best fit gives for the important exponent $`h(T)`$ a value of $`0.390(5)`$. Our results for $`h(T)`$ at different temperatures, figure 7, are in good agreement with those of huse ; bray , but have higher statistical accuracy and allow for the correction to scaling which was not considered in the earlier work. We have checked that our estimates for $`h(T)`$ do not change when we use sample sizes $`L`$ from $`20`$ to $`28`$; this means that we should be in the limit such that the $`h(T)`$ values represent the infinite size behaviour.

Estimates of $`w/z`$ from the analyses of the non-equilibrium susceptibility data have statistical fluctuations but do not show any systematic trend with $`T`$ over the fairly narrow range of $`T`$ values near $`T_g`$ that we have explored, so our estimate of $`w/z`$ is not sensitive to the exact choice for $`T_g`$. The mean value is $`w/z=0.48(10)`$. As a number of independent measurements have shown that the dynamic exponent $`z`$ is near $`6`$ (og ; bray ; pom  and our estimate given below at the end of the analysis), we can deduce from $`w/z`$ that the dynamic correction to scaling exponent $`w=2.9(6)`$. This turns out to be very close to the static exponent value $`\omega 2.8`$ estimated from an analysis of the size dependence of the spin glass susceptibility $`\chi (L)`$ pom . It is distinctly higher than the estimate $`\omega =0.8(4)`$ given by Bal . (An attempt to fit the $`q^2(t)`$ data with a correction exponent $`w/z`$ of the order of $`1/6`$ fails completely).

There is an independent check on the expected size of $`\omega `$. Estimates of $`\omega `$ can be obtained from series expansion data through $`\omega =\mathrm{\Delta }_1/\nu `$ klein  where $`\mathrm{\Delta }_1`$ the appropriate correction to scaling parameter for the series work. In dimension 4 a reliable value $`\omega 3.0`$ was obtained klein ; in dimension 3 the estimate by the same authors is $`3.0`$ with rather lower accuracy. In the sequence of Ising ferromagnets pelissetto  and for the sequence of percolation models perco , $`\omega `$ increases regularly with $`d_ud`$, where $`d_u`$ is the upper critical dimension ($`4`$ for Ising ferromagnets, $`6`$ for the percolation model as for ISGs). The high value of $`\omega `$, near $`3`$ in dimension 3 from the ISG series work, is consistent with the value from the present simulation analysis and is also consistent with there also being a general upward trend for $`\omega `$ with increasing $`d_ud`$ in ISGs. A much lower estimate for $`\omega `$ Bal  seems inconsistent with the series work.

A large correction exponent means that the correction factors tend towards $`1`$ quite rapidly as $`L`$ or $`t`$ increase, so the large size or long time limiting behaviour can be reached fairly readily. However for large $`\omega `$ it may not be safe to ignore analytic corrections to scaling Bal , which could perhaps explain the non-monotonic behaviour of $`T_L^{}`$ from the Binder parameter data seen above, figure 2.

For the estimates of $`x(T)`$, the measurement runs for $`q(t)`$ must be preceded by long annealing runs. It turns out that after an appropriate anneal the algebraic form with constant $`x(T)`$ sets in from rather short times (a few Monte Carlo steps), and the measurements are self averaging og ; rieger . Deviations from the asymptotic behaviour are so small that no meaningful analysis in terms of corrections to scaling could be carried out (see for instance the $`q(t)`$ data shown in og ). As a consequence it is relatively easy to obtain robust and accurate values of $`x(T)`$ on large samples. The data of references og  and rieger  are accurate and in excellent agreement with each other, figure 6. Our own data are in full agreement with these values so we have relied on the published measurements. In ISGs the relaxation function continues to be algebraic for temperatures below $`T_g`$, so there is no problem in defining $`x(T)`$, while above $`T_g`$ an additional relaxation term with a cutoff function form appears. Ogielski introduced phenomenologically the stretched exponential decay in this context to parametrise the cutoff function, and it has been frequently used since. For present purposes, at temperatures where the pure algebraic decay no longer holds, one is clearly above $`T_g`$.

Again, as for $`h(T)`$, $`x(T)`$ can be measured on almost arbitrarily large samples. This is an important advantage of methods where achieving strict thermal equilibrium is not obligatory.

Numerically the most demanding of the three measurements is that of the size dependent equilibrium spin glass susceptibility $`\chi (L,T)`$. We have relied on the excellent data of KY , which extend up to sizes $`L`$ from $`16`$ to $`24`$ depending on the temperature. The method of allowing for corrections to finite size scaling is outlined in pom . Essentially, a reliable estimate of the effective temperature dependent exponent $`\eta (T)`$ can be obtained from a fit to $`\chi (L,T)`$ including the correction factor

$$\chi (L,T)/L^2L^{\eta (T)}(1+constL^\omega )$$
(13)
with an effective exponent $`\eta (T)`$ at and below $`T_g`$ and a temperature independent $`\omega `$. The correction term turns out to be more important than the statistical error only for $`L<8`$. (We have checked that the raw equilibrium $`\chi (L,T)`$ data exhibited by cruz  are in excellent agreement with those of KY ). The $`\chi (L,T)`$ results can be taken to provide good infinite size limit estimates for $`\eta (T)`$, figure 7, with individual errors of about $`\pm 0.03`$ including statistical and systematic contributions. Above $`T_g`$ the $`\chi (L,T)/L^2`$ curves start to bend down with large $`L`$, as would be expected from the general scaling form.

The results for $`h(T)`$ and $`h^{}(T)`$ as defined above are shown in figure 8. There is a clean intersection at $`T=1.195(15)`$ which we can take to be a reliable estimate for $`T_g`$, essentially free of corrections to finite size scaling because in each data set care has been taken to use as large samples as possible, and to allow for corrections to finite size scaling whenever appropriate. The clean intersection should set at rest any lingering suspicion marinari  that there is not a true transition in this spin glass at a critical temperature $`T_g`$ close to that estimated by Ogielski og  many years ago.

We obtain the corresponding exponents from an analysis of the various data, with $`T_g`$ in hand, Table I. Consistent and precise values are obtained using finite size scaling for $`\chi (L,T)`$ KY  with correction terms pom , and from the large sample data shown by Ogielski og . The value estimated for $`z`$ directly from the critical values of $`h`$ and $`x`$ is $`5.65(15)`$ which can be compared with $`6.0(8)`$ from og . The major source of error for all the exponents (as reflected in the error bars) is the estimated uncertainty in $`T_g`$.

Finally, if we plot $`\eta (d)`$ as a function of dimension $`d`$ we obtain the curve shown in Figure 9, where we have included the present result (Table 1) for dimension $`3`$, and values from jerome , endo ; ian , klein ; ian  and the mean field value($`\eta =0`$) for dimensions $`2,4,5`$ and $`6`$ respectively. (We can note that $`\eta `$ remains well defined even in dimension 2 where $`T_g`$ tends to zero jerome ). The curve is a guide for the eye. The overall shape of $`\eta (d)`$ is strikingly similar to that for the analagous series of $`\eta (d)`$ values in the percolation model adler . The present estimate for $`\eta (d=3)`$ intrapolates smoothly into the sequence of values obtained independently for the other dimensions.

## VI Conclusion

In conclusion, we have once again considered the problem of obtaining reliable estimates for the ordering temperature of the binomial interaction ISG in dimension three and for the associated critical exponents. It is very important to allow for corrections to finite size scaling for this particular system. It is clear from figures 1 and 2 that Binder parameter crossing points in this ISG are strongly affected by substantial corrections to finite size scaling even up to and beyond sizes as large as $`L=20`$. For the spin glass susceptibility $`\chi `$, corrections to scaling are also present but their influence virtually dies out by $`L8`$. From non-equilibrium susceptibility data we obtain an estimate for the value of the dynamic correction to scaling exponent which is large, $`w=2.9(6)`$, and very similar to a previous estimate of the static correction to scaling exponent $`\omega `$ pom .

From an analysis of non-equilibrium data on large samples together with equilibrium susceptibilites pom ; lwb , our estimates for $`T_g=1.195(15)`$, and the associated static critical exponents, are in excellent agreement with those of Ogielski og  who made dynamic and static measurements on even larger samples (up to $`L=64`$). They are close to those of Palassini and Caracciolo pala , who carefully extrapolated to infinite size, explicitly allowing for corrections to finite size scaling. They are consistent with high quality series results singh , which come from a method complementary to and independent of the Monte Carlo numerical simulations. When Binder parameter data are extrapolated to large $`L`$, the $`T_g`$ estimate obtained is imprecise because of uncertainties in the corrections to finite size scaling, but within the errors it is consistent with the estimates using the other techniques. There is however a puzzling residual disagreement with the work of Ballesteros et al Bal  concerning the precise value of the critical temperature $`T_g`$.

We gratefully acknowledge the use of the computing facilities of IDRIS.


text:
# Can A Gravitational Quadruple Lens Produce 17 images?

## 1 Introduction

Gravitational lensing phenomenon is ubiquitous throughout the universe and is poised to bear a substantial weight of the explorations of diverse objects from planetary systems (Bennett and Rhie, 2000) to the universe itself (Kaiser, 1998; Tyson et al., 1998; Schechter, 2000). The utility of gravitational lensing dates back to the inception of Einstein’s general relativity (Einstein 1911, 1916; see Weinberg 1972; SEF 1992) whose first confirmation as the theory of gravity was offered by the measurements of the astrometric effect of the gravitational lensing by the sun of the background stars (Eddington 1919; see Weinberg 1972). The more dramatic lensing effect is the multiplicity of the images which arises due to the fact that there can be multiple null geodesics (photon paths) that connect a distance radiation emission source and an observer because of the focusing effect of the intervening lensing masses.

The first gravitationally lensed multiple image object to be discovered was double quasar Q0957+561A,B (Walsh et al., 1979) whose two images separated by about 6 arcseconds plank the elliptical lensing (cD) galaxy located at about 5 arcseconds from the brighter image A. A topological argument showed that a smooth extended mass distribution generates an odd number of images (Burke, 1981)111This paper has a special feature: It has only one page! In the case of $`n`$-point lenses, Poincaré-Hopf index of Burke’s vector field is $`(1n)`$, and so the total number of images is even when $`n`$ is odd and vice versa. Also, see section 2, and note the statement on logarithmic singularities in SEF, pp. 175., and so with the first double quasar surfaced the problem of missing images which may be too dim to be identified in the middle of the stars of the lensing galaxy or which may have been altered due to the hard nucleus effect of the stars of the lensing galaxy in which the assumption of continuum mass distribution has to be modified.

The effective lensing range of a point mass is given by the Einstein ring radius which is roughly the geometric mean of the Schwarzschild radius and the (reduced) distance of the lensing mass. So, the lensing by a point mass is short-ranged and its singular nature can dominate the lensing behavior within its lensing range. The Schwarschild radius of a star is $`10^5`$ (light) sec, the size of the universe is $`4\times 10^{17}`$ sec, and so the Einstein ring radius of a star at a cosmological distance is $`2\times 10^6`$ sec. Since the typical distance between two stars in a galaxy is $`10^8`$ sec, the effective number of stars that affect the local lensing behavior can be a few to many depending on the local surface mass density of the galaxy where the image under consideration is; the particle effect of the stars in a galaxy as a continuuum is ususally discussed under the subject title of quasar microlensing because the stellar granularity 222The effect of galactic granularity in a cluster lensing is of mesolensing where the astrometric shifts due to the galaxy nuggets can not be ignored. CL0024+1654 is a well-knwon example. As is common to any physical system, this intermediate scale phenomenon requires careful studies(Broadhurst et al., 2000). mainly affects the photometric behavior rather than the global astrometric behavior or the time delay (Young , 1981). Careful studies of the optical light curves of Q0957+561A,B have shown the absence of microlensing of time scales 40 - 120 days (Gil-Merino et al., 2000). The radio time delays of the double quasar at $`4`$cm and $`6`$cm are somewhat inconsistent with the optical time delay $`420`$ days (Haarsma et al., 1999), and the residual systematics may be an indicator of radio activities to be learned.

Quasar microlensing may be considered gravitational scintillation by a random distribution of an unknown effective number of point lenses (Young , 1981). In a microlensing of a star, the lensing range of a microlensing star is order of $`1`$AU (the reduced distance of the lens is $`10^6`$ times smaller than in quasar microlensing), and so the gravitational lensing involves only a single scattering and the lens is a low multiplicity $`n`$-point lens whose multiplicity is determined by the number of gravitationally bound objects within the lensing range of the total lensing mass. These low multiplicity point lenses are particularly interesting because the simplicity allows discoveries and measurements of the lensing systems which may be compact stars, planetary systems, or free-floating planets that are difficult to probe with other means. Observations of protoplanetary disks (Robberto et al. , 1999), density wake fields, and dust rings (Ozernoy et al., 2000) are being actively pursued, and microlensing will offer extensive statistics on planets as the final products of the planet formations and evolutions. Microlensing will also offer extensive statistics on stellar remnant black holes in the Galaxy adding our knowledge on the end point of the stellar evolution theories, the IMF, and the early history of the Galaxy.

Microlensing is a relatively recent experimental phenomenon (Alcock, 2000) advanced with large format CCDs, fast affordable computers, high density data storage devices, and our total ignorance on dark matter species (Paczynski, 1986). The experimental necessity to interpret each lensing event exactly, as opposed to statistical nature of quasar microlensing (Wyithe and Turner, 2000), affords examining of the low multiplicity point lens systems to exhausting details. We describe the simple behavior of the simplest lenses that is shared by arbitrary $`n`$-point lens systems: the $`n`$-point lens equation is almost analytic and that imposes topological constraints on the caustic curves in equations (18) and (12). In the case of binary lenses, the contraints are strong enough to specify all the topological configurations of the caustic curves (Rhie, 2001). With a bit of systematic toiling, one can also comb through the behavior of the class of triple lenses which is particularly useful for multiple planetary systems (Bennett and Rhie, 2001) or circumbinary planets (Bennett et al., 2000).

In the following sections, we summarize the properties of the $`n`$-point lens equations, construct quadruple lenses that produce up to 15 images and argue that 15 is the maximum possible number of images of any quadruple lenses. We discuss the notion of image domains and illustrate the correlation between the limit points as the “markers” of the positive image domains and the number of positive images of a point source. We conclude without a rigorous mathematical proof that the number of positive images of an $`n`$-point lens is bounded by the number of finite limit points $`2(n1):n>1`$ and so the total number of images of a point source by an $`n`$-point lens can not exceeed $`5(n1):n>1`$. We leave it as a conjecture, having no shortage of margin, for an algebraic topologist or a combinatorian (perhaps with a finite Erdös number) who may find it interesting to complete the counting. We construct “necklace” point mass configurations that have positive images at $`(n+1)`$ distinct limit points to show that the lower bound of the maximum number of images is $`3(n+1)`$. We comment on (thick) Einstein ring images in relation to the image domains of the interior of caustic curves.

## 2  The $`n`$-point Lens Equation

A gravitational lens of total mass $`M`$ consisting of $`n`$ gravitationally bound point masses whose interdistances are negligibly small compared to the distances of the center of mass to the positions of the observer ($`D_1`$) and the radiation emission source ($`D_2`$) is described by the following (two-dimensional linear) lens equation where $`\omega `$ is the position of the unlensed radiation emission source and $`z`$ is the position of an image of the radiation emission source generated by the point lens masses $`ϵ_jM`$ located at $`x_j:j=1,\mathrm{},n`$.

$$\omega =z\underset{j=1}{\overset{n}{}}\frac{ϵ_j}{\overline{z}\overline{x}_j};\underset{j=1}{\overset{n}{}}ϵ_j=1.$$
(1)
Here the lens plane where the two-dimensional position variables reside and perpendicular to the line of sight of the (unlensed) source star passes through the center of mass of the lens elements; this choice is useful when the focus is on the physics of the lensing system such as a microlensing planetary system. The lens equation has been normalized to be dimensionless such that the unit distance scale of the lens plane parameterized by $`\omega `$ or $`z`$ is set by the Einstein ring radius of the total mass $`M`$ located at the center of mass.

$$1=R_E=\sqrt{4GMD};\frac{1}{D}=\frac{1}{D_1}+\frac{1}{D_2}$$
(2)
where $`D`$ is the reduced distance. (As is familiar from the reduced mass in mechanics, the reduced distance is smaller than the smaller of the two distances $`D_1`$ and $`D_2`$.)

The $`n`$-point lens equation (1) can be embedded in an ($`n^2+1`$)-th order analytic polynomial equation (Witt, 1990), and so the number of images of a point source can not exceed ($`n^2+1`$) because an ($`n^2+1`$)-th order analytic polynomial (a complex function of only $`z`$ or $`\overline{z}`$) always has ($`n^2+1`$) solutions. A single lens produces two images, a binary lens produces three or five images, a triple lens produces four, six, eight, or ten images, and it has been a puzzle whether a quadruple lens can produce up to $`n^2+1=17`$ images (Mao et al., 1997). The mininum number of images of an $`n`$-point lens is $`n+1`$ and one can count them easily by looking at the lens equation for $`\omega =\mathrm{}`$: there is one positive image at $`\mathrm{}`$ and one negative image at each of the $`n`$ lens positions. The smooth continuum mass in Burke’s odd number theorem has one (unlensed) positive image near the source for large $`\omega `$.

### 2.1 Linear Differential Properties of the Lens Equation

The linear differential behavior of the lens equation plays an important role in lensing and is described by the Jacobian matrix of the lens equation.

$$𝒥=\left(\begin{array}{c}1\overline{\kappa }\\ \kappa 1\end{array}\right);\kappa \overline{\omega }\frac{\overline{\omega }}{z}=\underset{j=1}{\overset{n}{}}\frac{ϵ_j}{(zx_j)^2}$$
(3)
If an image of an infinitesimally small source is formed at $`z`$, the parity of the image is given by the sign of the determinant $`J`$ of the Jacobian matrix and the magnification of the image is given by $`A=|J|^1`$.

$$J(z)=1|\kappa (z)|^2;J(z)1$$
(4)
The set of solutions to $`J(z)=0`$ occupies a central position in lensing because that is where the magnification $`A(z)`$ becomes large and the lensing signals are most obvious. The lens equation is stationary where its Jacobian determinant $`J`$ vanishes, and so the curve $`J=0`$ is called the critical curve. Images flip the parity across the critical curve because $`J`$ changes its sign. In fact, the images with opposite parities form pairs that are mapped to the same source positions under the lens equation.

Let’s consider a real function $`y=x^2`$ for an easy illustration of the criticality: it is stationary or the linear derivative vanishes at $`x=0`$; two points $`x=\pm \delta x`$ are mapped to one point $`y=\delta x^2`$; the Jacobian determinant $`J_1`$ of the $`1\times 1`$ Jacobian matrix $`𝒥_1`$ is $`J_1=𝒥_1=dy/dx`$, and $`J_1>0`$ for $`x=\delta x`$ ($`>0`$) and $`J_1<0`$ for $`x=\delta x`$; $`y=0`$ on which the critical point $`x=0`$ is mapped under the real function divides the $`y`$-axis into two regions: $`\{y|y>0\}`$ where each point $`y`$ has two solutions $`x=\pm \sqrt{y}`$ and $`\{y|y>0\}`$ where each point $`y`$ has no solutions. The set of points onto which the critical curve is mapped under the lens equation is called the caustic curve and each caustic point defines a direction (transverse to the caustic curve) that is locally analogous to the $`y`$-axis where the caustic point is at $`y=0`$.

In the neighborhood of a caustic point, point sources on one side of the caustic point (say $`y>0`$) produce pairs of images with opposite parities across the corresponding critical point (analogous to $`x=0`$) while point sources on the other side of the caustic point (say $`y<0`$) do not. Thus, the caustic curve defines the boundary across which the number of images changes by two (or a multiple of two at the intersection points). In other words, in the neighborhood of the critical curve the lens equation is two-to-one, and two images appear from or disappear into the critical curve as the source crosses the caustic curve. The italicized statement is well-known, and the focus of this paper is to discuss what is underlying the behavior of caustic curves as the boundaries of the caustic domains and to convey that caustic curves are easily comprehensible benign curves despite the hostile appearance of the spikiness near the cusps. It all starts from a simple observation that the Jacobian matrix is determined by one analytic function $`\kappa (z)`$. The lens equation is a function of both variables $`z`$ and $`\overline{z}`$, but the derivatives are (anti)analytic. So, complex coordinates are useful for the $`n`$-point lens equation as are spherical coordinates for spherically symmetric systems.

The Analytic Function $`\kappa (z)`$ Defines a Potential Flow on the Lens Plane. An analytic function describes a two-dimensional harmonic flow of an inviscid fluid whose streamlines are conserved except at the “sources” and “sinks” given by the poles and zeroes of the analytic function (Landau and Lifshitz, 1982).

$$\overline{}\kappa =0^2\kappa =0$$
(5)
The $`n`$ lens positions are the poles of $`\kappa `$ and are double poles. In the neighborhood of a lens position, $`zx_j`$,

$$\kappa \frac{ϵ_j}{(zx_j)^2}.$$
(6)
The zeros of $`\kappa `$ are referred to as (finite) limit points because $`J`$ takes the maximum value $`1`$ where $`\kappa =0`$. Equation (3) shows that there are $`2(n1)`$ limit points and they are simple zeros unless degenerate. If $`z_{}`$ is a limit point, near the limit point $`zz_{}`$,

$$\kappa \kappa (z_{})(zz_{}).$$
(7)
The infinity behaves as a double zero. If $`z\mathrm{}`$,

$$\kappa \underset{1}{\overset{n}{}}\frac{ϵ_j}{z^2}=\frac{1}{z^2}.$$
(8)
If the phase angle of $`\kappa `$ is $`2\phi `$,

$$\kappa =|\kappa |e^{2i\phi },$$
(9)
then $`\phi `$ changes by $`\pi `$ around each (finite) limit point, $`2\pi `$ around the infinity, and $`2\pi `$ around each lens position where the positive orientation is counterclockwise as usual. So, the phase angle of $`\kappa `$ defines the stream lines or field lines that flow from the lens positions to the finite limit points and infinity. See figure 1 for the field lines of a binary lens. The magnitude $`|\kappa |`$ behaves as the potential and decreases monotonically along the field lines from $`|\kappa |=\mathrm{}`$ at the lens positions to $`|\kappa |=0`$ at the limit points. The equipotential curves ($`|\kappa |=`$ constant) are smooth because the lens equation is smooth everywhere except at the poles ($`\kappa =\mathrm{}`$) and are orthogonal to the field lines ($`\phi =`$ constant). The total number of field lines that cross an arbitrary equipotential curve ($`|\kappa |=`$ constant) which may consist of multiple loops is $`\mathrm{\Delta }\phi =2\pi n`$.

The Critical Curve is but an Equipotential Curve given by $`|\kappa |=1`$. The critical curve is an equipotential curve (with the potential value $`1`$) since $`|\kappa |=1`$ where $`J=0`$, and it is easy to visualize the configurations because of the orderly behavior of the family of equipotential curves. If we consider the lens plane as a large two-sphere (with one point $`\mathrm{}`$ at infinity – the well-known notion of one-point compactification), the equipotential curves and field lines can be considered latitudes and longitudes of a multiple-pole sphere with the north poles at the lens positions (“sources”) and the south poles at the limit points and infinity (“sinks”). Check the critical curves in this paper for the predictable configurations. The normals to the equipotential curves and field lines form a right-handed basis vector field.

$$\{|\kappa |,\phi \}2\{\overline{}|\kappa |,\overline{}\phi \}$$
(10)
The multiplicity of the “sources” and “sinks” implies that there are “forks” in the streamlines where the equipotential curves bifurcate. The bifurcation points (or saddle points: $`\kappa =0`$) are four-prong vertices (or $`\times `$-points), and the basis vectors are not uniquely defined at the bifurcation points. There are three bifurcation points in every binary lens (Rhie, 2001). For higher $`n`$ ($`>2`$), the number varies with the lens parameters.

A critical loop (a connected closed curve with $`J=0`$) can be assigned an integer or a half-integer because of the analyticity of $`\kappa `$. If the closed curve encloses $`N_\times `$ lens positions and $`N_{}`$ limits points, we define the associated topological charge $`e`$ as follows. (It is valid for an arbitrary closed curve that is analytically connected to the critical curve: the residue theorem.)

$$e=N_\times \frac{N_{}}{2}$$
(11)
Then the net flux of the field lines that flow out through the critical loop is $`\mathrm{\Delta }\phi =2\pi e`$. The charge of a limit point is negative ($`e=1/2`$) because the field lines flow in while the charge of a lens position is positive ($`e=1`$). The relative sign reflects the opposite orientations of the critical loops around a limit point and a lens position, which derives from the analytic behavior of $`\kappa `$ around a zero and a pole in equations (7) and (6). The critical curve of a lens can consist of many loops, and the topological charges associated with the critical loops are subject to a constraint because the total number of field lines is determined by the number of lens elements $`n`$: $`\mathrm{\Delta }\phi =2\pi n`$.

$$\underset{j}{}|e_j|=n$$
(12)
If $`z`$ is a Critical Point, $`\omega (z)`$ is called a Caustic Point. Since the lens equation is smooth (and so continuous) in the neighborhood of the critical curve, the caustic curve $`\{w(z)||\kappa (z)|=1\}`$ shares the connectivity of the critical curve and is smooth except at the stationary points where cusps form. The caustic curve is made of the same number of loops the critical curve is made of and bifurcates where the critical curve bifurcates. So a caustic loop is associated with the same topological charge of the corresponding critical loop, and we will see that the charge $`|e|`$ measures the smooth rotation of the tangent of the caustic loop. At the cusps, the tangents change the directions abruptly by $`\pi `$. The critical points that are mapped to cusps under the lens equation are called precusps.

In the case of a single lens (Schwarzschild lens), every point on the critical curve (Einstein ring) is a precusp, and the point caustic may be considered a degenerate cusp. If we look at the caustics in this paper, the cusps are just a finite number of punctuations on the otherwise smooth curves. So our guiding thought for an intuitive or impressionistic understanding (or handwaving) would be: How badly are the critical curves of $`n`$-point lenses deviated from the circular Einstein ring of a single lens? If we look at the critical curves in this paper, we do find some comforting similarities in the segments around lens positions. That is where the critical curve is more or less tangential to the critical direction ($`E_{}`$ defined below). The segments around limit points include parts where the critical curve is more or less normal to the critical direction. If we look at the critical curve in figure LABEL:fig\_idomain, its (extrinsic) curvature sign changes around the precusps corresponding to the off-the-axis cusps, but $`E_{}=ie^{i\phi }`$ smoothly rotates with increasing $`\phi `$ along the critical curve. So, it should not be surprising that the connected critical curve behaves differently around the limit points from around the lens positions.

The critical condition $`J=0`$ implies that (at least) one of the eigenvalues of the Jacobian matrix $`𝒥`$ vanishes on the critical curve. The $`2\times 2`$ matrix $`𝒥`$ has two eigenvalues 333The Jacobian determinant is the product of the two eigenvalues, and $`J=0`$ result in two possible types of critical curves: “tangential critical curve” where the critical curve is tangent to the critical direction ($`E_{}`$; see section 2) and “radial critical curve” where the critical curve is normal to $`E_{}`$ (SEF pp231). A single (Schwarzschild) lens has only a “tangential critical curve” which is the Einstein ring (the second factor of $`J`$ is positive and $`0`$) and each point of the “tangential critical curve” is a precusp. In general, the critical curve of an $`n`$-point lens is a mixture of the parts that are more or less tangential to $`E_{}`$ and parts that are more or less normal to $`E_{}`$., and $`\lambda _{}`$ vanishes on the critical curve.

$$\lambda _\pm =1\pm |\kappa |$$
(13)
The eigenvectors ($`\pm e_\pm `$: blind to the senses) are conveniently described by the basis vectors ($`e_\pm `$) (with definite orientations).

$$e_+\left(\genfrac{}{}{0pt}{}{e^{i\phi }}{e^{i\phi }}\right),e_{}\left(\genfrac{}{}{0pt}{}{ie^{i\phi }}{ie^{i\phi }}\right);e_\pm =\sqrt{2}$$
(14)
An arbitrary vector $`dz`$ can be decomposed in terms of the basis vectors: $`dz=dz_+E_++dz_{}E_{}`$ where $`dz_\pm `$ are real and $`E_\pm `$ denote the upper components of $`e_\pm `$.

The Critical Direction is $`\pm E_{}`$. The Caustic Curve is Tangent to $`\pm E_+`$. Let $`dz`$ be an arbitrary displacement from a critical point $`z_c`$, then $`\lambda _{}(z_c)=0`$, and $`d\omega =dz_+\lambda _+E_+=2dz_+E_+`$. So the lens equation is critical or stationary in the eigendirection $`\pm E_{}`$: if $`dz_+=0`$, then $`d\omega =0`$. If $`dz=dz_{}E_{}`$, the lens equation restricted to the critical direction is quadratic (Rhie and Bennett, 2000).

$$\delta \omega _{2c}=\frac{1}{2}\overline{}Jdz_{}^2=\frac{1}{4}dz_{}^2\left(_+JE_++_{}JE_{}\right),$$
(15)
where $`dJ=dz_+_+J+dz_{}_{}J`$. In an exact analogy to the real function $`y=x^2`$ we discussed above, the equation (15) has two solutions ($`\pm dz_{}`$) only if $`\delta \omega _{2c}/\overline{}J>0`$. Then, the remaining question is: Where is $`\overline{}J`$ pointing? In the $`z`$-plane, it is a normal vector to the critical curve (or any equipotential curve): it is normal outward to the critical curve with positive topological charge $`e>0`$, and normal inward to the one with $`e<0`$. In order to determine the relative direction of $`\overline{}J`$ in the neighborhood of the caustic curve (in the $`\omega `$-plane), we need to find a reference direction that is well defined intrinsically. We first note that the caustic curve is always tangent to the eigendirection $`\pm E_+`$. So the normal component $`(_{}JE_{})`$ determines on which side of the caustic curve a point source at $`\omega (z_c)+\delta \omega _{2c}`$ produces two images at $`z_c\pm dz_{}`$ and we can use the (extrinsic) curvature vector of the caustic curve as the reference direction.

If $`dz=dt`$ is tangent to the critical curve ($`d\phi >0`$), $`d\omega (z)`$ is tangent to the causic curve and is in the direction of $`(_{}JE_+)`$.

$$\delta \omega _1=2dt_+E_+=\frac{|dt|}{|\overline{}J|}_{}JE_+$$
(16)
The second order term shows that the caustic curve bends in the direction of $`(_{}JE_{})`$.

$$\delta \omega _{2t}=\frac{1}{2}\overline{}\overline{\kappa }dt^2=\frac{1}{4}|dt|^2(_+JE_+_{}JE_{})$$
(17)
So, if $`\pm E_+`$ define the $`x`$-axis such that the caustic curve is convex upward, then the quadratic equation (15) has two solutions where $`y\delta \omega _{2c}>0`$ and none where $`y<0`$. We also observe from (16) and (17) that the tangent to the caustic curve rotates clockwise ($`d\phi >0`$) unless $`_{}J=0`$. See figure 2 for an illustration. So, in order to form a caustic domain that produces many images, the caustic curve will have to wind around many times (here always counterclockiwise) intersecting itself and nesting. The caustic curve in figure 2 has total winding number 5.

The lens equation restricted to the critical curve is stationary ($`\delta \omega _1=0`$) where $`_{}J=0`$, and the caustic curve develops a cusp. If $`m`$ is the number of the cusps of a caustic loop, the tangent of the caustic loop rotates by $`2\pi |e|`$ (clockwise) in total on the smooth segments and flips by $`m\pi `$ at the cusps resulting in the net rotation of $`2\pi q=m\pi 2\pi |e|`$ where $`q`$ is the winding number.

$$q=\frac{m}{2}|e|$$
(18)
A caustic loop with $`q2`$ intersects itself and can form a caustic domain where a point source generates $`2q`$ more images than outside the caustic loop. $`|e|`$ is easily determined from the corresponding critical loop and subject to the constraint $`_j|e_j|=n`$ in equation (12). Caustic loops intersect and nest each other, and the number of images jumps by two at each crossing of the caustic loop. The case $`q=1`$ is like an $`m`$-gon (polygon with $`m`$ vertices) and $`|e|`$ amounts to the total internal angle in units of $`2\pi `$. A critical loop enclosing a limit point has $`|e|=1/2`$, and so the corresponding caustic loop has “internal angle” $`\pi `$ and is triangular: $`(q,m)=(1,3)`$. The caustic curve of a quadruple lens in figure 2 consists of two triangular caustics and one with 12 cusps: $`(q,m)=(3,12)`$, and it satisfies the constraint $`_j|e_j|=1/2+1/2+3=4`$.

A corollary will be the theorem of even number of cusps (SEF, pp213). From equations (18) and (12),

$$\underset{j}{}\frac{m_j}{2}=n+\underset{j}{}q_j.$$
(19)
Since the RHS is an integer, $`_jm_j`$ is an even number. Thus in an $`n`$-point lens, the number of caustic loops with odd number of cusps is always even. In a binary, trioids: $`(q,m)=(1,3)`$ appear in a pair that are reflection symmetric with respect to the lens axis. The triple lens in figure 1 in Rhie (1997) consists of two caustic loops with odd number of cusps: $`(q,m)=(1,3)`$ and $`(2,9)`$. The quadruple lens in figure 2 includes two trioids.

## 3 Can a Quadruple Lens Produce 17 Images?

The fact that there are triple lenses that can produce up to 10 images was discovered while we were investigating the location of a point source that produces images at the limit points (Rhie, 1997). There are four (finite) limit points in a triple lens, and the number of images of the point source that produces images at the four limit points must be at least ten because there are always two more negative images than positive images. Since a triple lens can produce no more than 10 images, a correlation between the maximum possible number of images and the limit points was suspected. We found that a point source can produce an image at each of $`2(n1)`$ limit points only in the cases of $`n=2`$ and $`n=3`$.

If a quadruple lens were to produce 17 images with all the positive images at the limit points, there would have to be 7 limit points because there are always three more negative images in a quadruple lens. A quadruple lens has only 6 (finite) limit points. We find quadruple lenses with 6 distinct (or non-degenerate) limit points that can have caustic domains $`𝒟^5`$ where a source produces 15 images (Fig. 2, Fig. 7), and the limit points ($`z_{}`$) are not mapped onto one point in the $`\omega `$-plane (Fig. 3). We can see from figure 2 and figure 3 that one extra limit point and associated $`|e|=1/2`$ caustic curve could make a domain $`𝒟^6`$ with 17 images, but that would violate the constraint $`_j|e_j|=4`$ (and the theorem of even number of cusps).

The quadruple lens in figure 4 is triple-like and allows a source position that has an image at every limit point. The limit points are (doubly) degenerate due to the high degree of symmetry of the lens configuration, and the caustic curve is relatively simple accommodating only up to $`𝒟^2`$ where a source generates 9 images (3 positive images and 6 negative images); a source at the center produces 3 images on each three-fold symmetry axis. If the three lens elements off the center are on a circle of radius $`a`$, the 3 positive images of the source at the center are formed at the limit points when $`a=2^{\frac{1}{6}}=0.89`$.

The quadruple lens in figure 5 has five distinct limit points because the one at the center is doubly degenerate, and the highest degree caustic domain is $`𝒟^4`$ where a source generates 13 images. As the two limit points at the center in figure 3 merges, the two intersecting triangular caustics in figure 2 “merge and unwrinkle” into a quadratic caustic in figure 5. An intermediate step after the “merger” is shown for two cases in figure 7: the domains $`𝒟^5`$ are split and diminished in size. If we consider a family of lenses in which four equal mass lens elements are equally spaced on a circle of radius $`a`$ (such that $`x_1=a`$), the source at the center of mass generates 5 positive images at the limit points when $`a=3^{\frac{1}{4}}/2=0.658`$. The total number of images is 13, and the source is in $`𝒟^4`$ which is the highest degree caustic domain for the family of regular rhombus lens configurations.

What is implied is that the maximum number of images can be obatined when the configuration of the lens elements is largely symmetric to maximize the gravitational interference of all the lens elements as a whole but not so symmetric as to generate degenerate limit points. When the lens elements are dispersed, the gravitational interference is fractionized and the lens system effectively behaves like a linear sum of lower number of lens systems (Fig. 8). When the limit points are degenerate, each degeneracy costs 2 images. In fact, the limit points do not have to be degenerate to cost images: When they are close enough, the critical loops around two limit points merge to form one connected positive image domain. The residual highest degree domains $`𝒟^5`$ in figure 7 disappear quickly when the lens configurations become more regular. One trivial nontheless worthy point to note is that a degenerate limit point can host only one image of a point source (property of the lens equation) even though the degeneracy has its direct bearing on the $`e`$-charge of a loop (such as critical loop) enclosing it (linear differential behavior of the lens equation). Then, the significance of the limit points must lie in that they are the “markers” of positive image domains, which we will conjecture to be the case after briefly examining the notion of image domains.

The Images of the Caustic Curve Defines the Image Domains. We are convinced by now that the caustic curve of a lens defines the caustic domains, and it is clear what we must have meant by the fact that the caustic domain $`𝒟^5`$ of the quadruple lens in figure 2 generates 15 images. There are 15 image domains in the image plane whose boundaries are mapped onto the boundary of the $`𝒟^5`$ under the lens equation. The image plane with the domain boundaries of the quadruple lens is a bit of mess to use for an illustration even though we can say that it is rather pretty. So, let’s examine the class of binary lenses that offers a simple laboratory case for the notion of image domains; the caustic curves are simple curves (winding number $`q=1`$) and there are only two types of caustic domains: outside for 3 images and inside for 5 images. The fact that the caustic curve is sandwiched between the domains of 3 images and 5 images implies that the caustic curve itself generates 4 images. In other words, there are 4 curves in the image plane that are mapped onto the caustic curve, one of which is obviously the critical curve. The critical curve is smooth as we are all familiar with by now, but the others do have cuspy points. The lens equation at the cuspy points is smooth without criticality, and so the corresponding cusps (kinky points) on the caustic curve generate kinky points in the $`z`$-plane.

Figure 9 shows the image domains of a binary lens with a connected caustic with six cusps that are defined by the 4 images curves of the caustic curve shown in the inset. Smooth one is the critical curve, and the other 3 curves are tangent to the critical curve at the precusps (each curve is attached to the critical curve at two precusps). The 3 curves are smooth at the precusps on the critical curve where the lens equation (restricted to the critical curve) is stationary, hence each of the 3 curves have four cuspy points that are the images of the cusps of the caustic curve where the lens equation is non-critical. The four curves that are the image curves of the caustic curve form 5 image domains around the critical curve. They are where the 5 images of an emission source inside the caustic curve form and consist of 3 negative image domains (labelled $`𝒟_{}^1`$) and 2 positive image domains ($`𝒟_+^1`$). The limit points are marked by $``$ as is the case throughout this paper, and they are inside the two positive domains $`𝒟_+^1`$. The entire lens plane is divided into 8 domains by the image curves of the caustic curve: $`5𝒟^1`$’s and $`3𝒟^0`$’s, and they account for 3 images for a source outside the caustic curve and 5 images for one inside. If we consider our Local Group as a binary lens made of the Galaxy and M31 at a cosmological distance, an elliptical source galaxy can fill the inside caustic domain to produce a large ring image with two holes which we may refer to as “thick Einstein ring” with two holes.

In the neighborhood of a precusp the image domains are divided into four sections, which is related to the fact that the lens equation restricted to the critical direction is cubic at a precusp and produces 3 solutions when the source is inside the caustic and 1 solution when the source is outside. The four domains around a precusp consist either of 3 positive image domains and 1 negative image domain (say, type $`+`$), or of 1 positive image domain and 3 negative image domains (type $``$). The two cusps on the lens axis are of type $`+`$, and the four cusps off the lens axis are of type $``$. When a source crosses the caustic through a cusp of type $`+`$ ($``$), two positive (negative) images are created and move in a tangential direction to the critical curve, and one positive (negative) image crosses the critical curve changing its parity. The net parity remains the same as it does when the caustic crossing is away from a cusp. It should be worth reemphasizing that the well-known statement that two images of opposite parities appear or disappear at a caustic crossing which is strictly valid for line caustic crossings has to be relaxed for cusp crossings because the pair of images that appear or disappear have the same parities (rh97).

These microscopic discussions which may seem irrelevant or annoyingly nitpicking fit together rather naturally once we draw the image domains, and simple intuitive interpretations fall out: A source near a cusp inside the caustic has 3 highly magnified images because they are all near the critical curve; if the source moves along a line caustic, two images remain large because they are near the critical curve and the 3rd one shrinks because it is away from the critical curve; the 3rd image did not come into the picture as the other two small images when we considered the quadratic equation for two large images in equation (15) because of its distance from the critical curve. The mysterious behavior of the parities of the two images that appear or disapper across the caustic curve we discussed above (and once was a source of caustic comments from referees) is just a partner swapping process among the three highly magnified images whose process one can comprehend easily by considering a source moving along the caustic curve past a cusp; the caustic curve near a cusp is reflection-symmetric with respect to its tangent vector (eigenvector $`E_+`$) at the cusp, and a source approaching the cusp along the symmetry axis from inside the caustic generates two images that are reflection-symmetric with respect to the symmetry axis and one image on the symmetry axis; the two images off the symmetry axis ($`E_+`$) are on the same side of the critical curve and so have the parity.

Figure 10 shows the image domains of a case where the caustic curve with six cusps is split into two quadroids. The long arrow indicates 3 large image curves of the 4-cupsed caustic centered around $`(0,1)`$. If a quasar with its host galaxy covers this caustic loop, there will be four images of the quasar and the finite size galaxy will form a ring image with one hole. Two of them are postive images and the other two are negative images. The fourth image curve of the caustic is inside the the image ring of the other caustic loop and is marked by the shorter arrow. This small image domain generates the fifth image (negative image) and has a high probability to become a missing image because of its smallness. What we note is that when the critical curve splits, the number of positive image domains for each caustic loop is never more than the number of positive image domains of the connected caustic curve which is the same as the number of the finite limit points.

In the case of lensing by an elliptical projected mass distribution (zeroth approximation of a galaxy lens), the caustic consists of one quadroid (or astroid), and the fifth image of a quasar lensed by a galaxy is more likely to be missing due to the brightness of the sizeable lensing galaxy. However, the color difference of the lensing galaxy and the lensed quasar may reveal the small fifth image in high resolution images444PG 1115+080, B1608+656, and B1938+666 (Fig. 2 of Kochanek et al. 2000) are caustic covering lensing events and one can infer from their ring images the characteristic assemblage of the image domains around critical curves. The image domain configurations depend on the lensing mass distributions, and the apparent (image) rings depend also on the source position with respect to the caustic. . The position of the fifth image in relation to the image ring will be one of the clinching constraints in the reconstruction of the image domains and so the lens characteristics. Draw a circle of radius $`\sqrt{0.3}=0.548`$ centered at the lens position (marked $`\times `$) of the smaller “thick Einstein ring” (arrowed) in figure 10 and note that the circle goes through just about the middle of the image ring donut. The circle is the Einstein ring of the lensing mass $`0.3`$ as a single lens. Thus, one can more or less read off the lensing mass from the image ring. The apprearance of the observed image ring of a host galaxy of a quasar depends on the position of the source with respect to the caustic, of course, but the distribution of the bright core images of the quasar offers information on the point source position. We will discuss the image domains, “thick Einstein rings”, and Einstein rings of extended mass distribution lenses elsewhere (Rhie, 2001).

## 4  A Conjecture on the Number of Images

We have meandered through the properties of the caustic curve, the image curves of the caustic curve, image domains, the critical curve and smoothly rotating Jacobian eigenvectors on the critical curve, and the limit points and lens positions as the “sinks” and “sources” of the potential defined by $`|\kappa (z)|`$. The number of images jumps across the caustic curve because of the pairs of positive and negative images across the critical curve. Thus, if we refer to the image curves of the caustic curve other than the critical curve as precaustic curves, the critical curve is “padded” from both sides by the positive and negative image domains defined by the precaustic curves. Precaustic curves are tangentially attached to the critical curve at designated precusps, and their $`J`$ values are either positive (positive precaustic curve) or negative (negative precaustic curve) except at the precusps where they are tangent to the critical curve. The number of positive images of a given caustic is given by the number of positive image domains that are mapped to the caustic domain by the lens equation, and figure 9 and 10 illustrate that each positive image domain defined by a positive precaustic curve is associated with a finite limit point. So, a conjecture follows.

1. Conjecture: The number of positive images of an $`n`$-point lens is bounded by the number of finite limit points $`2(n1)`$ for $`n>1`$. When $`n=1`$, there is one limit point at $`\mathrm{}`$ and a point source has one positive image as is well known.
2. Corollary: The total number of images of an $`n`$-point lens can not exceed $`5(n1)`$ when $`n>1`$. When $`n=1`$, there are two images.
3. The maximum number of images of an $`n`$-point lens can not be less than $`3(n+1):n3`$.

The corollary follows from the conjecture because there are always $`(n1)`$ more negative images in an $`n`$-point lens. In order to prove the third item, we consider a lens consisting of $`n`$ equal masses equally spaced on a circle of radius $`a`$ can produce an image at every distinct limit point. We choose a source at the center of circle to maximally utilize the symmetry and examine the cases of $`n=2^k`$. The lens equation (1) becomes as follows where $`z_jzx_j`$.

$$z=\underset{1}{\overset{n}{}}\frac{1}{n}\frac{1}{\overline{z}_j}\overline{F}_n$$
(20)
For $`n=2`$, we let $`x_{1,2}=\pm a`$, and

$$\overline{F}_2=\frac{\overline{z}}{\overline{z}^2a^2}.$$
(21)
For $`n=4`$, we add two equal masses at $`x_{3,4}=\pm ia`$, and

$$\overline{F}_4=\overline{F}_2(a^2)+\overline{F}_2(a^2)=\frac{\overline{z}^3}{\overline{z}^4a^4}.$$
(22)
The pattern is clear: $`\overline{F}_8=\overline{F}_4(a^4)+\overline{F}_4(a^4)`$, etc, and

$$\overline{F}_n==\frac{\overline{z}^{n1}}{\overline{z}^na^n}.$$
(23)
The limit points $`z_{}`$ are the zeros of $`F_n`$.

$$0=z^{n2}(z^2+(n1)a^n)$$
(24)
The center $`z_{}=0`$ is an (n-2)-th order zero and satifies the lens equation (and so an image). Other $`n`$ zeros are on the circle of radius $`r_{}=(n1)^{1/n}a:r_{}<a`$.

$$z_{}=r_{}e^{\frac{\pi }{n}},r_{}e^{\frac{3\pi }{n}},\mathrm{}r_{}e^{\frac{(n1)\pi }{n}}$$
(25)
They are images when $`a=a_n`$.

$$a_n=n^{\frac{1}{2}}(n1)^{\frac{n2}{n}}$$
(26)
As $`n`$ becomes large, $`a_n`$ converges to $`1`$. Since the $`n`$-point “necklace” lens has $`(n+1)`$ distinct limit points, the source at the center produces no less than $`(n+1)`$ positive images, or equivalently no less than $`3n+1`$ images in total. One image is at the center, and $`3n`$ images are on three circles of radius $`r_{}`$ ($`n`$ positive images), $`r_<<a`$ ($`n`$ negative images), and $`r_>>a`$ ($`n`$ negative images). In the case of $`n=4`$, $`a=0.658`$, $`r_{}=\sqrt{3}/2=0.866`$, $`r_<=0.5`$, and $`r_>=1/2+\sqrt{7}/4=1.16`$. We have analysed the cases where $`n`$ is power of 2, but the (discrete) axial symmetry of the “necklace” lenses gurantees that $`3n+1`$ is the lower bound of the maximum total number of images of arbitrary $`n`$-point lenses.

The formula $`3n+1`$ produces the correct maximum possible number of images $`10`$ for $`n=3`$ because it happens to be that the two formulae for the number of limit points $`2(n1)`$ and $`n+1`$ coincide when $`n=3`$. Mao et al. (1997) numerically examined some cases of the family of “necklace” lenses and observed that the maximum number of images was $`3n+1`$. The formula $`3n+1`$ predicts too many images for $`n=2`$ and $`n=1`$ because the number of the limit points is less than $`(n+1)`$ and 2 less images for $`n=4`$. Here we have shown that $`3n+1`$ is the lower bound of the maximum possible number of images of an $`n`$-point lens for $`n3`$ by investigating source positions that generate positive images at the limit points as we did in rh97.