text:
# Least squares fitting of circles and lines

## 1 Introduction

Fitting simple contours (primitives) to experimental data is one of the basic problems in pattern recognition and computer vision. The simplest contours are linear segments and circular arcs. The need of approximating scattered points by a circle or a circular arc arises in physics , biology and medicine , archeology , industry , etc. The problem was studied since at least early sixties , and attracted much more attention in recent years due to its importance in image processing .

We study the least squares fit (LSF) of circles and circular arcs. This method is based on minimizing the mean square distance from the circle to the data points. Given $`n`$ points $`(x_i,y_i)`$, $`1in`$, the objective function is defined by

$$=\underset{i=1}{\overset{n}{}}d_i^2$$
(1.1)
where $`d_i`$ is the Euclidean (geometric) distance from the point $`(x_i,y_i)`$ to the circle. If the circle satisfies the equation

$$(xa)^2+(yb)^2=R^2$$
(1.2)
where $`(a,b)`$ is its center and $`R`$ its radius, then

$$d_i=\sqrt{(x_ia)^2+(y_ib)^2}R$$
(1.3)
The minimization of (1.1) is a nonlinear problem that has no closed form solution. There is no direct algorithm for computing the minimum of $``$, all known algorithms are either iterative and costly or approximative by nature.

The purpose of this paper is a general study of the problem. It consists of three parts. In the first one we address fundamental issues, such as the existence and uniqueness of a solution, the right choice of parameters to work with, and the general behavior of the objective function $``$ on the parameter space. These issues are rarely discussed in the literature, but they are essential for understanding of advantages and disadvantages of practical algorithms. The second part is devoted to the geometric fit – the minimization of the sum of squares of geometric distances, which is given by $``$. Here we evaluate three most popular algorithms (Levenberg-Marquardt, Landau, and Späth) and develop a new approach. In the third part we discuss an algebraic fit based on approximation of $``$ by a simpler algebraic function that can be minimized by a direct, noniterative algorithm. We compare three such approximations and propose some more efficient numerical schemes than those published so far.

Additional information about this work can be found on our web site .

## 2 Theoretical aspects

The very first questions one has to deal with in many mathematical problems are the existence and uniqueness of a solution. In our context the questions are: Does the function $``$ attain its minimum? Assuming that it does, is the minimum unique? These questions are not as trivial as they may appear.

2.1 Existence of LSF. The function $``$ is obviously continuous in the circle parameters $`a,b,R`$. According to a general principle, a continuous function always attains a minimum (possibly, not unique) if it is defined on a closed and bounded (i.e., compact) domain. Our function $``$ is defined for all $`a,b`$ and all $`R0`$, so its domain is not compact. For this reason the function $``$ fails to attain its minimum in some cases.

For example, let $`n3`$ distinct points lie on a straight line. Then one can approximate the data by a circle arbitrarily well and make $``$ arbitrarily close to zero, but since no circle can interpolate $`n3`$ collinear points, we will always have $`>0`$. Hence, the least squares fit by circles is, technically, impossible. For any circle one can find another circle that fits the data better. The best fit here is given by the straight line trough the data points, which yields $`=0`$. If we want the LSF to exist, we have to allow lines, as well as circles, in our model, and from now on we do this.

One can prove now that the function $``$ defined on circles and lines always attains its minimum, and so the existence of the LSF is guaranteed. A detailed proof is available in .

We should note that if the data points are generated randomly with a continuous probability distribution (such as normal or uniform), then the probability that the LSF returns a line, rather than a circle, is zero. This is why lines are often ignored and one practically works with circles only. On the other hand, if the coordinates of the data points are discrete (such as pixels on a computer screen), then lines may appear with a positive probability and have to be reckoned with.

2.2 Uniqueness of LSF. This question is not trivial either. Even if $``$ takes its minimum on a circle, the best fitting circle may not be unique, as several other circles may minimize $``$ as well. We could not find such examples in the literature, so we provide our own here.

Examples of multiple minima. Let four data points $`(\pm 1,0)`$ and $`(0,\pm 1)`$ make a square centered at the origin. We place another $`k4`$ points identically at the origin $`(0,0)`$ to have a total of $`n=k+4`$ points. This configuration is invariant under a rotation through the right angle about the origin. Hence, if a circle minimizes $``$, then by rotating that circle through $`\pi /2`$, $`\pi `$, and $`3\pi /2`$ we get three other circles that minimize $``$ as well.

Thus, we get four distinct best fitting circles, unless either (i) the best circle is centered at the origin or (ii) the best fit is a line. So we need to show that neither is the case. This involves some simple calculations. If a circle has radius $`r`$ and center at $`(0,0)`$, then $`=4(1r^2)+kr^2`$. The minimum of this function is attained at $`r=4/(k+4)`$, and it equals $`_0=4k/(k+4)`$. Also, the best fitting line passes through the origin and gives $`_1=2`$. On the other hand, consider the circle passing through three points $`(0,0)`$, $`(0,1)`$, and $`(1,0)`$. It only misses two other points, $`(1,0)`$ and $`(0,1)`$, and it is easy to see that it gives $`<2`$, which is less than $`_1`$ and $`_0`$ whenever $`k4`$. Hence, the best fit is a circle not centered at the origin, and so our construction works.

Figure 1: A data set (left) for which the objective function (right) has four minima.

Figure 1 illustrates our example, it gives the plot of $``$ with four distinct minima (the plotting method is explained below). By changing the square to a rectangle one can make $``$ have exactly two minima. By replacing the square with a regular $`m`$-gon and increasing the number of identical points at $`(0,0)`$ one can construct $``$ with exactly $`m`$ minima for any $`m3`$, see details in .

Of course, if the data points are generated randomly with a continuous probability distribution, then the probability that the objective function $``$ has multiple minima is zero. In particular, small random perturbations of the data points in our example on Fig. 1 will slightly change the values of $``$ at its minima, so that one of them will become a global (absolute) minimum and three others – local (relative) minima.

We note, however, that while the cases of multiple global minima are indeed exotic, they demonstrate that the global minimum of $``$ may change discontinuously if one slowly moves data points.

2.3 Local minima of the objective function. Generally, local minima are undesirable, since they can trap iterative algorithms and lead to false solutions.

We have investigated how frequently the function $``$ has local minima (and how many). In our experiment, $`n`$ data points were generated randomly with a uniform distribution in the unit square $`0<x,y<1`$. Then we applied the Levenberg-Marquard algorithm (described below) starting at 1000 different, randomly selected initial conditions. Every point of convergence was recorded as a minimum (local or global) of $``$. If there were more than one point of convergence, then one of them was the global minimum and the others were local minima. Table 1 shows the probabilities that $``$ had 0,1,2 or more local minima for $`n=5,\mathrm{},100`$ data points.

|  | 5 | 10 | 15 | 25 | 50 | 100 |
| --- | --- | --- | --- | --- | --- | --- |
| 0 | 0.879 | 0.843 | 0.883 | 0.935 | 0.967 | 0.979 |
| 1 | 0.118 | 0.149 | 0.109 | 0.062 | 0.031 | 0.019 |
| $`2`$ | 0.003 | 0.008 | 0.008 | 0.003 | 0.002 | 0.002 |

Table 1. Probability of 0, 1, 2 or more local minima of $``$ when $`n=5,\mathrm{},100`$ points are randomly generated in a unit square.

We see, surprisingly, that local minima only occur with probability $`<15`$%. The more points are generated, the less frequently the function $``$ has any local minima. Multiple local minima turn up with probability even less than 1%. The maximal number of local minima we observed was four, it happened only a few times in almost a million random samples we tested.

Generating points with a uniform distribution in a square produces completely “chaotic” samples without any predefined pattern. This is, in a sense, the worst case scenario. If we generate points along a circle or a circular arc (with some noise), then local minima virtually never occur. For example, if $`n=10`$ points are sampled along a 90o circular arc of radius $`R=1`$ with a Gaussian noise at level $`\sigma =0.05`$, then the probability that $``$ has any local minima is as low as $`0.001`$.

Therefore, in typical applications the objective function $``$ is very likely to have a unique (global) minimum and no local minima. Does this mean that a standard iterative algorithm, such as the steepest descent or Gauss-Newton or Levenberg-Marquardt, would converge to the global minimum from any starting point? Unfortunately, this is not the case, as we demonstrate next.

2.4 Plateaus and valleys on the graph of the objective function. In order to examine the behavior of $``$ we found a way to visualize its graph. Even though $`(a,b,R)`$ is a function of three parameters, it happens to be just a quadratic polynomial in $`R`$ when the other two variables are fixed. So it has a unique global minimum in $`R`$ that can be easily found. If we denote

$$r_i=\sqrt{(x_ia)^2+(y_ib)^2}$$
(2.1)
then the minimum of $``$ with respect to $`R`$ is attained at

$$\widehat{R}=\overline{r}:=\frac{1}{n}\underset{i=1}{\overset{n}{}}r_i$$
(2.2)
This allows us to eliminate $`R`$ and express $``$ as a function of $`a,b`$:

$$=\underset{i=1}{\overset{n}{}}(r_i\overline{r})^2$$
(2.3)
A function of two variables can be easily plotted. This is exactly how we did it on Fig. 1.

Figure 2: A simulated data set of 50 points.

Now, Fig. 2 presents a typical random sample of $`n=50`$ points generated along a circular arc (the upper half of the unit circle $`x^2+y^2=1`$) with a Gaussian noise added at level $`\sigma =0.01`$. Fig. 3 shows the graph of $``$ plotted by MAPLE in two different scales. One can clearly see that $``$ has a global minimum around $`a=b=0`$ and no local minima. Fig. 4 presents a flat grey scale contour map, where darker colors correspond to deeper parts of the graph (smaller values of $``$).

Figure 3: The objective function $``$ for the data set shown on fig. 2: a large view (left) and the minimum (right).

Fig. 3 shows that the function $``$ does not grow as $`a,b\mathrm{}`$. In fact, it is bounded, i.e. $`(a,b)_{\mathrm{max}}<\mathrm{}`$. The boundedness of $``$ actually explains the appearance of large nearly flat plateaus and valleys on Fig. 3 that stretch out to infinity in some directions. If an iterative algorithm starts somewhere in the middle of such a plateau or a valley or gets there by chance, it will have hard time moving at all, since the gradient of $``$ almost vanishes there. We indeed observed conventional algorithms getting “stuck” on flat plateaus or valleys in our tests. The new algorithm we propose below does not have this drawback.

Figure 4: A grey-scale contour map of the objective function $``$. Darker colors correspond to smaller values of $``$. The minimum is marked by a cross.

Second, there are two particularly interesting valleys that stretch roughly along the line $`a=0`$ on Figs. 3 and 4. One of them, corresponding to $`b<0`$, has its bottom point at the minimum of $``$. The function $``$ slowly decreases along the valley as it approaches the minimum. Hence, any iterative algorithm starting in that valley or getting there by chance should, ideally, find its way downhill and arrive at the minimum of $``$.

The other valley corresponds to $`b>0`$, it is separated from the global minimum of $``$ by a ridge. The function $``$ slowly decreases along this valley as $`b`$ grows. Hence, any iterative algorithm starting in this valley or getting there “by accident” will be forced to move up along the $`b`$ axis, away from the minimum of $``$, all the way to $`b=\mathrm{}`$.

If an iterative algorithm starts at a randomly chosen point, it may go down into either valley, and there is a good chance that it descends into the second valley and then escapes to infinity. For the sample on Fig. 2, we applied the Levenberg-Marquardt algorithm starting at a randomly selected point in the square $`5\times 5`$ about the centroid $`x_c=\frac{1}{n}x_i`$, $`y_c=\frac{1}{n}y_i`$ of the data. We found that the algorithm escaped to infinity with probability about 50%.

Unfortunately, such “escape valleys” are inevitable. We have proved  that for almost every data sample there was a pair of valleys stretching out in opposite directions, so that one valley descends to the minimum of $``$, while the other valley descends toward infinity.

We now summarize our observations. There are three major ways in which conventional iterative algorithms may fail to find the minimum of $``$:

* when they converge to a local minimum;
* when they slow down and stall on a nearly flat plateau or in a valley;
* when they diverge to infinity along a decreasing valley.

We will not attempt here to deal with local minima of $``$ causing the failure of type (a), since local minima occur quite rarely. But the other two types of failures (b) and (c) can be drastically reduced, if not eliminated altogether, by adopting a new set of parameters, which we introduce next.

2.5 Choice of parametrization. The trouble with the natural circle parameters $`a,b,R`$ is that they become arbitrarily large when the data are approximated by a circular arc with low curvature. Not only does this lead to a catastrophic loss of accuracy when two large nearly equal quantities are subtracted in (1.3), but this is also ultimately responsible for the appearance of flat plateaus and descending valleys that cause failures (b) and (c) in iterative algorithms.

We now adopt a parametrization used in , in which the equation of a circle is

$$A(x^2+y^2)+Bx+Cy+D=0$$
(2.4)
Note that this gives a circle when $`A0`$ and a line when $`A=0`$, so it conveniently combines both types of contours in our model. The parameters $`A,B,C,D`$ must satisfy the inequality $`B^2+C^24AD>0`$ in order to define a nonempty circle or a line. Since the parameters only need to be defined up to a scalar multiple, we can impose a constraint

$$B^2+C^24AD=1$$
(2.5)
If we additionally require that $`A0`$, then every circle or line will correspond to a unique quadruple $`(A,B,C,D)`$ and vice versa. For technical reasons, though, we do not restrict $`A`$, so that circles have a duplicate parametrization, see (2.7) below. The conversion formulas between the natural parameters and the new ones are

$$a=\frac{B}{2A},b=\frac{C}{2A},R=\frac{1}{2|A|}$$
(2.6)
and

$$A=\pm \frac{1}{2R},B=2Aa,C=2Ab,D=\frac{B^2+C^21}{4A}$$
(2.7)
The distance from a point $`(x_i,y_i)`$ to the circle can be expressed, after some algebraic manipulations, as

$$d_i=2\frac{P_i}{1+\sqrt{1+4AP_i}}$$
(2.8)
where

$$P_i=A(x_i^2+y_i^2)+Bx_i+Cy_i+D$$
(2.9)
One can check that $`1+4AP_i0`$ for all $`i`$, see below, so that the function (2.8) is well defined.

The formula (2.8) is somewhat more complicated than (1.3), but the following advantages of the new parameters can outweigh the higher cost of computations.

First, the is no need to work with arbitrarily large values of parameters anymore. One can effectively reduce the parameter space to a box

$$\{|A|<A_{\mathrm{max}},|B|<B_{\mathrm{max}},|C|<C_{\mathrm{max}},|D|<D_{\mathrm{max}}\}$$
(2.10)
where $`A_{\mathrm{max}}`$, $`B_{\mathrm{max}}`$, $`C_{\mathrm{max}}`$, $`D_{\mathrm{max}}`$ can be determined explicitly . This conclusion is based on some technical analysis, and we only outline the main steps. Given a sample $`(x_i,y_i)`$, $`1in`$, let $`d_{\mathrm{max}}=\mathrm{max}_{i,j}\sqrt{(x_ix_j)^2+(y_iy_j)^2}`$ denote the maximal distance between the data points. Then we observe that (i) the distance from the best fitting line or circle to the centroid of the data $`(x_c,y_c)`$ does not exceed $`d_{\mathrm{max}}`$. It also follows from (2.2) that (ii) the best fitting circle has radius $`Rd_{\mathrm{max}}/n`$, hence $`|A|n/2d_{\mathrm{max}}`$. Thus the model can be restricted to circles and lines satisfying the two conditions (i) and (ii). Under these conditions, the parameters $`A,B,C,D`$ are bounded by some constants $`A_{\mathrm{max}}`$, $`B_{\mathrm{max}}`$, $`C_{\mathrm{max}}`$, $`D_{\mathrm{max}}`$. A detailed proof of this fact is contained in .

The effective boundedness of the parameters $`A,B,C,D`$ renders them a preferred choice for numerical calculations. In particular, when the search for a minimum is restricted to a closed box (2.10), it can hardly run into vast nearly flat plateaus that we have seen on Fig. 3.

Second, the objective function is now smooth in $`A,B,C,D`$ on the parameter space. In particular, no singularities occur as a circular arc approaches a line. Circular arcs with low curvature correspond to small values of $`A`$, lines correspond to $`A=0`$, and the objective function and all its derivatives are well defined at $`A=0`$.

Third, recall the two valleys shown on Figs. 3-4, which we have proved to exist for almost every data sample . In the new parameter space they become two halves of one valley stretching continuously across the hyperplane $`A=0`$ and descending to the minimum of $``$. Therefore, any iterative algorithm starting anywhere in that (now unique) valley would converge to the minimum of $``$ (maybe crossing the plane $`A=0`$ on its way). There is no escape to infinity anymore.

As a result, the new parametrization can effectively reduce, if not eliminate completely, the failures of types (b) and (c). This will be confirmed experimentally in the next section.

## 3 Geometric fit

3.1 Three popular algorithms. The minimization of the nonlinear function $``$ cannot be accomplished by a finite algorithm. Various iterative algorithms have been applied to this end. The most successful and popular are

* the Levenberg-Marquardt method.
* Landau algorithm .
* Späth algorithm

Here (a) is a short name for the classical Gauss-Newton method with the Levenberg-Marquardt correction . It can effectively solve any least squares problem of type (1.1) provided the first derivatives of $`d_i`$’s with respect to the parameters can be computed. The Levenberg-Marquardt algorithm is quite stable and reliable, and it usually converges rapidly (if the data points are close to the fitting contour, the convergence is nearly quadratic). Fitting circles with the Levenberg-Marquardt method is described in many papers, see, e.g. .

The other two methods are circle-specific. The Landau algorithm employs a simple fixed-point iterative scheme, nonetheless it shows a remarkable stability and is widely used in practice. Its convergence, though, is linear . The Späth algorithm makes a clever use of an additional set of (dummy) parameters and is based on alternating minimization with respect to the dummy parameter set and the real parameter set $`(a,b,R)`$. At each step, one set of parameters is fixed, and a global minimum of the objective function $``$ with respect to the other parameter set is found, which guarantees (at least theoretically) that $``$ decreases at every iteration. The convergence of Späth’s algorithm is, however, known to be slow .

The cost per iteration is about the same for all the three algorithms: the Levenberg-Marquardt requires $`12n+41`$ flops per iteration, Landau takes $`11n+5`$ flops per iteration and for Späth the flop count is $`11n+13`$ per iteration (in all cases, a prior centering of the data is assumed, i.e. $`x_c=y_c=0`$).

We have tested the performance of these three algorithms experimentally, and the results are reported in Section 3.3.

3.2 A new algorithm for circle fitting. In Section 2.5 we introduced parameters $`A,B,C,D`$ subject to the constraint (2.5). Here we show how the Levenberg-Marquardt scheme can be applied to minimize the function (1.1) in the parameter space $`(A,B,C,D)`$.

First, we introduce a new parameter – an angular coordinate $`\theta `$ defined by

$$B=\sqrt{1+4AD}\mathrm{cos}\theta ,C=\sqrt{1+4AD}\mathrm{sin}\theta ,$$

so that $`\theta `$ replaces $`B`$ and $`C`$. Now one can perform an unconstrained minimization of $`=d_i^2`$ in the three-dimensional parameter space $`(A,D,\theta )`$. The distance $`d_i`$ is expressed by (2.8) with

$`P_i`$ $`=`$ $`A(x_i^2+y_i^2)+\sqrt{1+AD}(x_i\mathrm{cos}\theta +y_i\mathrm{sin}\theta )+D`$
$`=`$ $`Az_i+Eu_i+D`$
where we denote, for brevity, $`z_i=x_i^2+y_i^2`$, $`E=\sqrt{1+4AD}`$, and $`u_i=x_i\mathrm{cos}\theta +y_i\mathrm{sin}\theta `$. The first derivatives of $`d_i`$ with respect to the parameters are

$$d_i/A=(z_i+2Du_i/E)R_id_i^2/Q_i$$

$$d_i/D=(2Au_i/E+1)R_i$$

$$d_i/\theta =(x_i\mathrm{sin}\theta +y_i\mathrm{cos}\theta )ER_i$$

where $`Q_i=\sqrt{1+4AP_i}`$ and

$$R_i=\frac{2(1Ad_i/Q_i)}{Q_i+1}$$

Then one can apply the standard Levenberg-Maquardt scheme, see details in . The resulting algorithm is more complicated and costly than the methods described in 3.1, it requires $`39n+40`$ flops and one trigonometric function call per iteration. But it converges in a fewer iterations than other methods, so that its overall performance is rather good, see the next section.

This approach has some pitfalls – the function $``$ is singular (its derivatives are discontinuous) when $`1+4AP_i=0`$ for some $`i`$ or when $`1+4AD=0`$. We see that

$$1+4AP_i=\frac{(x_ia)^2+(y_ib)^2}{R^2}$$

This quantity vanishes if $`x_i=a`$ and $`y_i=b`$, i.e. when a data point coincides with the circle’s center, and this is extremely unlikely to occur. In fact, it has never happened in our tests, so that we did not use any security measures against the singularity $`1+4AP_i=0`$.

On the other hand, we see that

$$1+4AD=\frac{a^2+b^2}{R^2}$$

which vanishes whenever $`a=b=0`$. This singularity turns out to be more serious – when the circle’s center computed iteratively approaches the origin, the algorithm often gets stuck because it persistently tries to enter the forbidden area $`1+4AD<0`$. We found a simple remedy: whenever the algorithm attempts to make $`1+4AD`$ negative, we shift the origin by adding a vector $`(u,v)`$ to all the data coordinates $`(x_i,y_i)`$, and then recompute the parameters $`(A,D,\theta )`$. The vector $`(u,v)`$ should be of size comparable to the average distance between the data points, and its direction can be chosen randomly. The shift had to be applied only occasionally and its cost was negligible.

Remark. Another convenient parametrization of circles (and lines) was proposed by Karimäki . He uses three parameters: the signed curvature ($`\rho =\pm 1/R`$), the distance of closest approach ($`d`$) to the origin, and the direction of propagation ($`\phi `$) at the point of closest approach. Karimäki’s parameters $`\rho `$, $`d`$, $`\phi `$ are similar to ours $`A`$, $`D`$, $`\theta `$, and can be also used in the alternative Levenberg-Marquardt scheme.

3.3 Experimental tests. We have generated $`10000`$ samples of $`n`$ points randomly with a uniform distribution in the unit square $`0<x,y<1`$. For each sample, we generated $`1000`$ initial guesses by selecting a random circle center $`(a,b)`$ in a square $`5\times 5`$ around the centroid $`(x_c,y_c)`$ of the sample and then computing the radius $`R`$ by (2.2). Every triple $`(a,b,R)`$ was then used as an initial guess, and we ran all the four algorithms starting at it.

For each sample and for each algorithm, we have determined the number of runs, from 1000 random initial guesses, when the iterations converged to the global minimum of $``$. Dividing that number by the total number of generated initial guesses, 1000, we obtained the probability of convergence to the minimum of $``$. We also computed the average number of iterations in which the convergence took place. At this stage, the samples for which the function $``$ had any local minima, in addition to the global minimum, were eliminated. The fraction of such samples was less than 15%, as one can see in Table 1. The remaining majority of samples were used to compute the overall characteristics of each algorithm: the grand average probability of convergence to the minimum of $``$ and the grand mean number of iterations the convergence took. We note that, since samples with local minima are eliminated, the probability of convergence will be a true measure of the algorithm’s reliability. All failures to converge will occur when the algorithm falters and cannot find any minima, which we consider as the algorithm’s fault.

This experiment was repeated for $`n=5,10,\mathrm{},100`$. The results are presented on Figures 5 and 6, where the algorithms are marked as follows: LAN for Landau, SPA for Späth, LMC for the canonical Levenberg-Maquardt in the $`(a,b,R)`$ parameter space, and LMA for the alternative Levenberg-Maquardt in the $`(A,D,\theta )`$ parameter space.

Figure 5: The probability of convergence to the minimum of $``$ starting at a random initial guess.

Figure 5 shows the probability of convergence to the minimum of $``$, it clearly demonstrates the superiority of the LMA method. Figure 6 presents the average cost of convergence, in terms of flops per data point, for all the four methods. The fastest algorithm is the canonical Levenberg-Marquardt (LMC). The alternative Levenberg-Marquardt (LMA) is almost twice as slow. The Späth and Landau methods happen to be far more expensive, in terms of the computational cost, than both Levenberg-Marquardt schemes.

Figure 6: The average cost of computations, in flops per data point.

The poor performance of the Landau and Späth algorithms is illustrated on Fig. 7. It shows a typical path taken by each of the four procedures starting at the same initial point $`(1.35,1.34)`$ and converging to the same limit point $`(0.06,0.03)`$ (the global minimum of $``$). Each dot represents one iteration. For LMA and LMC, subsequent iterations are connected by grey lines. One can see that LMA (hollow dots) heads straight to the target and reaches it in about 5 iterations. LMC (solid dots) makes a few jumps back and forth but arrives at the target in about 15 steps. On the contrary, SPA (square dots) and LAN (round dots) advance very slowly and tend to make many short steps as they approach the target (in this example SPA took 60 steps and LAN more than 300). Note that LAN makes an inexplicable long detour around the point $`(1.5,0)`$. Such tendencies account for the overall high cost of computations for these two algorithms as reported on Fig. 6.

Figure 7: Paths taken by the four algorithms on the $`ab`$ plane, from the initial guess at $`(1.35,1.34)`$ marked by a large cross to the minimum of $``$ at $`(0.06,0.03)`$ (a small cross).

Next, we have repeated our experiment with a different rule of generating data samples. Instead of selecting points randomly in a square we now sample them along a circular arc of radius $`R=1`$ with a Gaussian noise at level $`\sigma =0.01`$. We set the number of points to 20 and vary the arc length from $`5^\mathrm{o}`$ to $`360^\mathrm{o}`$. Otherwise the experiment proceeded as before, including the random choice of initial guesses. The results are presented on Figures 8 and 9.

Figure 8: The probability of convergence to the minimum of $``$ starting at a random initial guess.

We see that the alternative Levenberg-Marquardt method (LMA) is very robust – its displays a remarkable 100% convergence across the entire range of the arc length. The reliability of the other three methods is high (up to 95%) on full circles ($`360^\mathrm{o}`$) but degrades to 50% on half-circles. Then the conventional Levenberg-Marquardt (LMC) stays on the 50% level for all smaller arcs down to $`5^\mathrm{o}`$. The Späth method breaks down on $`50^\mathrm{o}`$ arcs and the Landau method breaks down even earlier, on $`110^\mathrm{o}`$ arcs.

Figure 9 shows that the cost of computations for the LMA and LMC methods remains low for relatively large arcs, but it grows sharply for very small arcs (below $`20^\mathrm{o}`$). The LMC is generally cheaper than LMA, but, interestingly, becomes more expensive on arcs below $`30^\mathrm{o}`$. The cost of the Späth and Landau methods is, predictably, higher than that of the Levenberg-Marquardt schemes, and it skyrockets on arcs smaller than half-circles making these two algorithms prohibitively expensive.

Figure 9: The average cost of computations, in flops per data point.

We emphasize, though, that our results are obtained when the initial guess is just picked randomly from a large square. In practice one can always find more sensible ways to choose an initial guess, so that the subsequent iterative schemes would perform much better than they did in our tests. We devote the next section to various choices of the initial guess and the corresponding performance of the iterative methods.

## 4 Algebraic fit

Generally, iterative algorithms for minimizing nonlinear functions like (1.1) are quite sensitive to the choice of the initial guess. As a rule, one needs to provide an initial guess close enough to the minimum of $``$ in order to ensure a rapid convergence.

The selection of an initial guess requires some other, preferably fast and non-iterative, procedure. In mass data processing, where speed is a factor, one often cannot afford relatively slow iterative methods, hence a non-iterative fit is the only option.

A fast and non-iterative approximation to the LSF is provided by the so called algebraic fit, or AF for brevity. We will describe three different versions of the AF below.

4.1 “Pure” algebraic fit. The first one, we call it AF1, is a very simple and old method, it has been known since at least 1972, see , and then rediscovered and published independently by many people . In this method, instead of minimizing the sum of squares of the geometric distances (1.1)–(1.3), one minimizes the sum of squares of algebraic distances

$`_1(a,b,R)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}[(x_ia)^2+(y_ib)^2R^2]^2`$ (4.11)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}(z_i+Bx_i+Cy_i+D)^2`$
where $`z_i=x_i^2+y_i^2`$ (as before), $`B=2a`$, $`C=2b`$, and $`D=a^2+b^2R^2`$. Now, differentiating $`_1`$ with respect to $`B,C,D`$ yields a system of linear equations:

$`M_{xx}B+M_{xy}C+M_xD`$ $`=`$ $`M_{xz}`$
$`M_{xy}B+M_{yy}C+M_yD`$ $`=`$ $`M_{yz}`$
$`M_xB+M_yC+nD`$ $`=`$ $`M_z`$
where $`M_{xx},M_{xy}`$, etc. denote moments, for example $`M_{xx}=x_i^2`$, $`M_{xy}=x_iy_i`$. Solving this system (by Cholesky decomposition or another matrix method) gives $`B,C,D`$, and finally one computes $`a,b,R`$.

The AF1 algorithm is very fast, it requires $`13n+31`$ flops to compute $`a,b,R`$. However, it gives an estimate of $`(a,b,R)`$ that is not always statistically optimal in the sense that the corresponding covariance matrix may exceed the Rao-Cramer lower bound, see . This happens when data points are sampled along a circular arc, rather than a full circle. Moreover, when the data are sampled along a small circular arc, the AF1 is very biased and tends to return absurdly small circles . Despite these shortcomings, though, AF1 remains a very attractive and simple routine for supplying an initial guess for iterative algorithms.

4.2 Gradient-weighted algebraic fit (GWAF). In the next subsections we will show how the “pure” algebraic fit can be improved at a little extra computational cost. First, we use again the circle equation (2.4) and note that the minimization of (4.11) is equivalent to that of

$$_1(A,B,C,D)=\underset{i=1}{\overset{n}{}}(Az_i+Bx_i+Cy_i+D)^2$$
(4.12)
under the constraint $`A=1`$. We will show that some other constraints lead to more accurate estimates.

The best results are achieved with the so called gradient-weighted algebraic fit, or GWAF for brevity. In its general form it goes like this. Suppose one wants to approximate scattered data with a curve described by an implicit polynomial equation $`P(x,y)=0`$, the coefficients of the polynomial $`P(x,y)`$ playing the role of parameters. The “pure” algebraic fit is based on minimizing

$$_\mathrm{a}=\underset{i=1}{\overset{n}{}}[P(x_i,y_i)]^2$$

where one of the coefficients of $`P`$ must be set to one to avoid the trivial solution in which all the coefficients turn zero. Our AF1 is exactly such a scheme. On the other hand, the GWAF is based on minimizing

$$_\mathrm{g}=\underset{i=1}{\overset{n}{}}\frac{[P(x_i,y_i)]^2}{P(x_i,y_i)^2}$$
(4.13)
here $`P(x,y)`$ is the gradient of the function $`P(x,y)`$. There is no need to set any coefficient of $`P`$ to one, since both numerator and denominator in (4.13) are homogeneous quadratic polynomials of parameters, hence the value of $`_\mathrm{g}`$ is invariant under the multiplication of all the parameters by a scalar. The reason why $`_\mathrm{g}`$ works better than $`_\mathrm{a}`$ is that we have, by the Taylor expansion,

$$\frac{|P(x_i,y_i)|}{P(x_i,y_i)}=d_i+𝒪(d_i^2)$$

where $`d_i`$ is the geometric distance from the point $`(x_i,y_i)`$ to the curve $`P(x,y)=0`$. Thus, the function $`_g`$ is simply the first order approximation to the classical objective function $``$ in (1.1).

The GWAF is known since at least 1974 . It was applied specifically to quadratic curves (ellipses and hyperbolas) by Sampson in 1982 , and recently became standard in computer vision industry . This method is well known to be statistically optimal, in the sense that the covariance matrix of the parameter estimates satisfies the Rao-Cramer lower bound . We plan to investigate the statistical properties of the GWAF for the circle fitting problem in a separate paper, here we focus on its numerical implementation.

In the case of circles, $`P(x,y)=A(x^2+y^2)+Bx+Cy+D`$ and $`P(x,y)=(2Ax+B,2Ay+C)`$, hence

$`P(x_i,y_i)^2`$ $`=`$ $`4Az_i^2+4ABx_i+4ACy_i+B^2+C^2`$ (4.14)
$`=`$ $`4A(Az_i+Bx_i+Cy_i+D)+B^2+C^24AD`$
and the GWAF reduces to the minimization of

$$_\mathrm{g}=\underset{i=1}{\overset{n}{}}\frac{[Az_i+Bx_i+Cy_i+D]^2}{[4A(Az_i+Bx_i+Cy_i+D)+B^2+C^24AD]^2}$$
(4.15)
This is a nonlinear problem that can only be solved iteratively, see some general schemes in . However, there are two approximations to (4.15) that lead to simpler and noniterative solutions.

4.3 Pratt’s approximation to GWAF. If data points $`(x_i,y_i)`$ lie close to the circle, then $`Az_i+Bx_i+Cy_i+D0`$, and we approximate (4.15) by

$$_2=\underset{i=1}{\overset{n}{}}\frac{[Az_i+Bx_i+Cy_i+D]^2}{B^2+C^24AD}$$
(4.16)
The objective function $`_2`$ was proposed by Pratt in 1987 , who clearly described its advantages over the “pure” algebraic fit. Pratt proposed to minimize $`_2`$ by using matrix methods, see below, which were computationally expensive. We describe a simpler yet absolutely reliable numerical algorithm below.

Converting the function $`_2`$ back to the original circle parameters $`(a,b,R)`$ gives the minimization problem

$$_2(a,b,R)=\frac{1}{4R^2}\underset{i=1}{\overset{n}{}}[x_i^2+y_i^22ax_i2by_i+a^2+b^2R^2]^2\mathrm{min}$$
(4.17)
In this form it was stated and solved by Chernov and Ososkov in 1984 . They found a stable and efficient noniterative algorithm that did not involve expensive matrix computations. Since 1984, this algorithm has been used in experimental nuclear physics. We propose an improvement of their algorithm below. The objective function (4.17) was also derived by Kanatani in 1998 , but he did not suggest any numerical method for minimizing it.

The minimization of (4.15) is equivalent to the minimization of the simpler function (4.12) subject to the constraint $`B^2+C^24AD=1`$. We write the function (4.12) in matrix form as $`_1=𝐀^T\mathrm{𝐌𝐀}`$, where $`𝐀=(A,B,C,D)^T`$ is the vector of parameters and $`𝐌`$ is the matrix of moments:

$$𝐌=\left(\begin{array}{cccc}M_{zz}& M_{xz}& M_{yz}& M_z\\ M_{xz}& M_{xx}& M_{xy}& M_x\\ M_{yz}& M_{xy}& M_{yy}& M_y\\ M_z& M_x& M_y& n\end{array}\right)$$
(4.18)
Note that $`𝐌`$ is symmetric and positive semidefinite (actually, $`𝐌`$ is positive definite unless the data points are interpolated by a circle or a line). The constraint $`B^2+C^24AD=1`$ can be written as $`𝐀^T\mathrm{𝐁𝐀}=1`$, where

$$𝐁=\left(\begin{array}{cccc}\hfill 0& \hfill 0& \hfill 0& \hfill 2\\ \hfill 0& \hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 1& \hfill 0\\ \hfill 2& \hfill 0& \hfill 0& \hfill 0\end{array}\right)$$
(4.19)
Now introducing a Lagrange multiplier $`\eta `$ we minimize the function

$$_{}=𝐀^T\mathrm{𝐌𝐀}\eta (𝐀^T\mathrm{𝐁𝐀}1)$$

Differentiating with respect to $`𝐀`$ gives $`\mathrm{𝐌𝐀}\eta \mathrm{𝐁𝐀}=0`$. Hence $`\eta `$ is a generalized eigenvalue for the matrix pair $`(𝐌,𝐁)`$. It can be found from the equation

$$\mathrm{det}(𝐌\eta 𝐁)=0$$
(4.20)
Since $`Q_4(\eta ):=\mathrm{det}(𝐌\eta 𝐁)`$ is a polynomial of the fourth degree in $`\eta `$, we arrive at a quartic equation $`Q_4(\eta )=0`$ (note that the leading coefficient of $`Q_4`$ is negative). The matrix $`𝐁`$ is symmetric and has four real eigenvalues $`\{1,1,2,2\}`$. In the generic case, when the matrix $`𝐌`$ is positive definite, by Sylvester’s law of inertia the generalized eigenvalues of the matrix pair $`(𝐌,𝐁)`$ are all real and exactly three of them are positive. In the special case when $`𝐌`$ is singular, $`\eta =0`$ is a root of (4.20). To determine which root of $`Q_4`$ corresponds to the minimum of $`_2`$ we observe that $`_2=𝐀^T\mathrm{𝐌𝐀}=\eta 𝐀^T\mathrm{𝐁𝐀}=\eta `$, hence the minimum of $`_2`$ corresponds to the smallest nonnegative root $`\eta _{}=\mathrm{min}\{\eta 0:Q_4(\eta )=0\}`$.

The above analysis uniquely identifies the desired root $`\eta _{}`$ of (4.20), but we also need a practical algorithm to compute it. Pratt  proposed matrix methods to extract all eigenvalues and eigenvectors of $`(𝐌,𝐁)`$, but admitted that those make his method a costly alternative to the “pure algebraic fit”. A much simpler way to solve (4.20) is to apply the Newton method to the corresponding polynomial equation $`Q_4(\eta )=0`$ starting at $`\eta =0`$. This method is guaranteed to converge to $`\eta _{}`$ by the following theoretical result:

###### Theorem 1

The polynomial $`Q_4(\eta )`$ is decreasing and concave up between $`\eta =0`$ and the first nonnegative root $`\eta _{}`$ of $`Q_4`$. Therefore, the Newton algorithm starting at $`\eta =0`$ will always converge to $`\eta _{}`$.

A full proof of this theorem is provided in . The resulting numerical scheme is more stable and efficient than the original Chernov-Ososkov solution . We denote the above algorithm by AF2.

The cost of AF2 is $`16n+16m+80`$ flops, here $`m`$ is the number of steps the Newton method takes to find the root of (4.20). In our tests, 5 steps were enough, on the average, and never more than 12 steps were necessary. Hence the average cost of AF2 is $`16n+160`$.

4.4 Taubin’s approximation to GWAF. Another way to simplify (4.15) is to average the variables in the denominator:

$$_3=\underset{i=1}{\overset{n}{}}\frac{[Az_i+Bx_i+Cy_i+D]^2}{[4A^2z+4ABx+4ACy+B^2+C^2]^2}$$
(4.21)
where

$$z=\frac{1}{n}\underset{i=1}{\overset{n}{}}z_i=\frac{1}{n}M_z,x=\frac{1}{n}M_x,y=\frac{1}{n}M_y$$

This idea was first proposed by Agin  but became popular after a publication by Taubin , and it is known now as Taubin method.

The minimization of (4.21) is equivalent to the minimization of $`_1`$ defined by (4.12) subject to the constraint

$$4A^2M_z+4ABM_x+4ACM_y+B^2n+C^2n=1$$
(4.22)
This problem can be expressed in matrix form as $`_1=𝐀^T\mathrm{𝐌𝐀}`$, see 4.3, with the constraint equation $`𝐀^T\mathrm{𝐂𝐀}=1`$, where

$$𝐂=\left(\begin{array}{cccc}4M_z& 2M_x& 2M_y& 0\\ 2M_x& n& 0& 0\\ 2M_y& 0& n& 0\\ 0& 0& 0& 0\end{array}\right)$$
(4.23)
is a symmetric and positive semidefinite matrix. Introducing a Lagrange multiplier, $`\eta `$ as in 4.3, we arrive at the equation

$$\mathrm{det}(𝐌\eta 𝐂)=0$$
(4.24)
Here is an advantage of Taubin’s method over AF2: unlike (4.20), (4.24) is a cubic equation for $`\eta `$, we write it as $`Q_3(\eta )=0`$. It is easy to derive from Sylvester’s law of inertia that all the roots of $`Q_3`$ are real and positive, unless the data points belong to a line or a circle, in which case one root is zero. As in 4.3, the minimum of $`_3`$ corresponds to the smallest nonnegative root $`\eta _{}=\mathrm{min}\{\eta 0:Q_3(\eta )=0\}`$.

Taubin  used matrix methods to extract eigenvalues and eigenvectors of the matrix pair $`(𝐌,𝐂)`$. A simpler way is to apply the Newton method to the corresponding polynomial equation $`Q_3(\eta )=0`$ starting at $`\eta =0`$. This method is guaranteed to converge to the desired root $`\eta _{}`$ since $`Q_3`$ is obviously decreasing and concave up between $`\eta =0`$ and $`\eta _{}`$ (note that the leading coefficient of $`Q_3`$ is negative). We denote the resulting algorithm by AF3.

The cost of AF3 is $`16n+14m+40`$ flops, here $`m`$ is the number of steps the Newton method takes to find the root of (4.24). In our tests, 5 steps were enough, on the average, and never more than 13 steps were necessary. Hence the average cost of AF3 is $`16n+110`$. This is 50 flops less than the cost of AF2.

4.5 Nonalgebraic (heuristic) fits. Some experimenters also use various simplistic procedures to initialize an iterative scheme. For example, some pick three data points that are sufficiently far apart and find the interpolating circle . Others place the initial center of the circle at the centroid of the data . Even though such “quick and dirty” methods are generally inferior to the algebraic fits, we will include two of them in our experimental tests for comparison. We call them TRI and CEN:

* TRI: Find three data points that make the triangle of maximum area and construct the interpolating circle.
* CEN: Put the center of the circle at the centroid of the data and then compute the radius by (2.2).

We note that our TRI actually requires $`𝒪(n^3)`$ flops and hence is far more expensive than any algebraic fit. In practice, though, one can often make a faster selection of three points based on the same principle .

4.5 Experimental tests. Here we combine the fitting algorithms in pairs: first, an algebraic (or heuristic) algorithm prefits a circle to the data, and then an iterative algorithm uses it as the initial guess and proceeds to minimize the objective function $``$. Our goal here is to evaluate the performance of the iterative methods described in Section 3 when they are initialized by various algebraic (or heuristic) prefits, thus ultimately we determine the quality of those prefits. We test all 5 initialization methods – AF1, AF2, AF3, TRI, and CEN – and all 4 iterative schemes – LMA, LMC, SPA, and LAN (in the notation of Section 3) – a total of $`5\times 4=20`$ pairs.

We conduct two major experiments, as we did in Sect. 3.3. First, we generate $`10000`$ samples of $`n`$ points randomly with a uniform distribution in the unit square $`0<x,y<1`$. For each sample, we determine the global minimum of the objective function $``$ by running the most reliable iterative scheme, LMA, starting at 1000 random initial guesses, as we did in Section 3.3. This is a credible (though, expensive) way to locate the global minimum of $``$. Then we apply all 20 pairs of algorithms to each sample. Note that no pair needs an initial guess, since the first algorithm in each pair is just designed to provide one.

After running all $`N=10000`$ samples, we find, for each pair $`[ij]`$ of algorithms ($`1i5`$ and $`1j4`$), the number of samples, $`N_{ij}`$, on which that pair successfully converged to the global minimum of $``$ (recall that the minimum was predetermined at an earlier stage, see above). The ratio $`N_{ij}/N`$ then represents the probability of convergence to the minimum of $``$ for the pair $`[ij]`$. We also find the average number of iterations the convergence takes, for each pair $`[ij]`$ separately.

Figure 10: The probability of convergence to the minimum of $``$ for 20 pairs of algorithms. The bar on the right explains color codes.

This experiment was repeated for $`n=10,20,50`$, and $`100`$ data points. The results are presented on Figures 10 and 11 by grey-scale diagrams. The bars of the right explain our color codes. For brevity, we numbered the algebraic/heuristic methods:

* 1 = AF1, 2 = AF2, 3 = AF3, 4 = TRI, and 5 = CEN

and labelled the iterative schemes by letters:

* a = LMA, b = LMC, c = SPA, and d = LAN

Each small square (cell) represents a pair of algorithms, and its color corresponds to a numerical value according to our code.

Fig. 10 shows the probability of convergence to the global minimum of $``$. We see that all the cells are white or light grey, meaning the reliability remains close to 100% (in fact, it is 97-98% for $`n=10`$, 98-99% for $`n=20`$ and almost 100% for $`n50`$.)

We conclude that for completely random samples filling the unit square uniformly, all five prefits are sufficiently accurate, so that any subsequent iterative method has little trouble converging to the minimum of $``$.

Figure 11: The cost in flops per data point. The bar on the right explains color codes.

Fig. 11 shows the computational cost for all pairs of algorithms, in the case of convergence. Colors represent the number of flops per data point, as coded by the bar on the far right. We see that the cost remains relatively low, in fact it never exceeds 700 flops per point (compare this to Figs. 6 and 9). The highest cost here is 684 flops per point for the pair TRI+LAN and $`n=10`$ points, marked by a cross.

The most economical iterative scheme is the conventional Levenberg-Marquardt (LMC, column b), which works well in conjunction with any algebraic (heuristic) prefit. The alternative Levenberg-Marquardt (LMA), Späth (SPA) and Landau (LAN) methods are slower, with LMA leading this group for $`n=10`$ and trailing it for larger $`n`$. There is almost no visible difference between the algebraic (heuristic) prefits in this experiment, except the TRI method (row 4) performs slightly worse than others, especially for $`n=100`$ (not surprisingly, since TRI is only based on 3 selected points).

One should not be deceived by the overall good performance in the above experiment. Random samples with a uniform distribution are, in a sense, “easy to fit”. Indeed, when the data points are scattered chaotically, the objective function $``$ has no pronounced minima or valleys, and so it changes slowly in the vicinity of its minimum. Hence, even not very accurate initial guesses allow the iterative schemes to converge to the minimum rapidly.

Figure 12: The probability of convergence to the minimum of $``$ for 20 pairs of algorithms. The bar on the right explains color codes.

We now turn to the second experiment, where data points are sampled, as in Sect. 3.3, along a circular arc of radius $`R=1`$ with a Gaussian noise at level $`\sigma =0.01`$. We set the number of points to $`n=20`$ and vary the arc length from $`20^\mathrm{o}`$ to $`180^\mathrm{o}`$. In this case the objective function has a narrow sharp minimum or a deep narrow valley, and so the iterative schemes depend on an accurate initial guess.

The results of this second experiment are presented on Figures 12 and 13, in the same fashion as those on Figs. 10 and 11. We see that the performance is quite diverse and generally deteriorates as the arc length decreases. The probability of convergence to the minimum of $``$ sometimes drops to zero (black cells on Fig. 12), and the cost per data point exceeds our maximum of 2000 flops (black cells on Fig. 13).

First, let us compare the iterative schemes. The Späth and Landau methods (columns c and d) become unreliable and too expensive on small arcs. Interestingly, though, Landau somewhat outperforms Späth here, while in earlier experiments reported in Sect. 3.3, Späth fared better. Both Levenberg-Marquartd schemes (columns a and b) are quite reliable and fast across the entire range of the arc length from $`180^\mathrm{o}`$ to $`20^\mathrm{o}`$.

Figure 13: The cost in flops per data point. The bar on the right explains color codes.

This experiment clearly demonstrates the superiority of the Levenberg-Marquardt algorithm(s) over fixed-point iterative schemes such as Späth or Landau. The latter two, even if supplied with the best possible prefits, tend to fail or become prohibitively expensive on small circular arcs.

Lastly, we extend this test to even smaller circular arcs (of 5 to 15 degrees) keeping only the two Levenberg-Marquardt schemes in our race. The results are presented on Figure 14. We see that, as the arc gets smaller, the conventional Levenberg-Marquardt (LMC) gradually loses its reliability but remains quite efficient, while the alternative scheme (LMA) gradually gives in speed but remains very reliable. Interestingly, both schemes take about the same number of iterations to converge, for example, on $`5^\mathrm{o}`$ arcs the pair AF2+LMA converged in 19 iterations, on average, while the pair AF2+LMC converged in 20 iterations. The higher cost of the LMA seen on Fig. 14 is entirely due to its complexity – one iteration of LMA requires $`39n+40`$ flops compared to $`12n+41`$ for LMC, see Sect. 3. Perhaps, the new LMA method can be optimized for speed, but we did not pursue this goal here. In any case, the cost of LMA per data point remains moderate, it is nowhere close to our maximum of 2000 flops (in fact, it always stays below 1000 flops).

Figure 14: The probability of convergence (left) and the cost in flops per data point (right) for the two Levenberg-Marquardt schemes.

We finally compare the algebraic (heuristic) methods. Clearly, AF1, TRI, and CEN are not very reliable, with CEN (the top row) showing the worst performance of all. Our winners are AF2 and AF3 (rows 2 and 3) whose characteristics seem to be almost identical, in terms of both reliability and efficiency.

The best pairs of algorithms in all our tests are these:

* AF2+LMA and AF3+LMA, which surpass all the others in reliability.
* AF2+LMC and AF3+LMC, which beat the rest in efficiency.

For small circular arcs, the pairs listed in (a) are somewhat slower than the pairs listed in (b), but the pairs listed in (b) are slightly less robust than the pairs listed in (a).

In any case, our experiments clearly demonstrate the superiority of the AF2 and AF3 prefits over other algebraic and heuristic algorithms. The slightly higher cost of these methods themselves (compared to AF1, for example) should not be a factor here, since this difference is well compensated for by the faster convergence of the subsequent iterative schemes.

Our experiments were done on Pentium IV personal computers and a Dell Power Edge workstation with 32 nodes of dual 733MHz processors at the University of Alabama at Birmingham. The C++ code is available on our web page .


format: text

text:
# Evolution of a Metal to Insulator Transition in Ca2-xNaxCuO2Cl2, as seen by ARPES

## I Introduction

The close proximity of an antiferromagnetic Mott insulating phase suggests that this insulator may be the natural starting point for understanding the unusually high superconducting transition temperatures as well as the unconventional normal state properties in the cuprates. The insulator to metal transition in strongly correlated materials is a well studied fundamental problem in solid state physics in its own right. Even by limiting ourselves to a one band Hubbard model, which may be an appropriate starting point for the half-filled cuprates, there remain several ways one could imagine the electronic structure evolving upon doping (see figure 1). One scenario is the chemical potential shifting to the top of the valence band similar to the case of a band insulator. An alternative scenario is states being created inside the Mott gap upon doping. These are the two most commonly discussed pictures with regard to the cuprates.

ARPES studies of the La2-xSrxCuO4 (LSCO) system show that spectral weight is transfered to low energy states, which lie inside the Mott gap, upon hole doping. In contrast, in the Bi2Sr2Ca1-xYxCu2O8+δ (Bi2212) and Nd2-xCexCuO4 (NCCO) systems there is work which shows that the chemical potential shifts in addition to a transfer of spectral weight. A potential problem for the Bi2212 system is the difficulty to obtain deeply underdoped samples. The oxychlorides, Sr2CuO2Cl2 and Ca2CuO2Cl2, are single layer cuprates that have provided the best photoemission data on the undoped insulating state to date. Thus the oxyhalides are ideal systems to study the metal to insulator transition with photoemission, as high quality single crystals, which cleave as easily as the Bi2212 system, can now be grown through the metal to insulator transition.

Recently, we have reported that the chemical potential shifts upon doping the insulator, Ca2CuO2Cl2(CCOC), with Na. Here we present the full body of evidence which supports arguments for such a shift of the chemical potential in the Ca2-xNaxCuO2Cl2 system. We examine the valence band and photon energy dependence of Na-doped CCOC and compare it with that found in the insulator. By making a detailed comparison of the low energy spectra we unambiguously show that the high energy pseudogap is in fact a result of the d-wave-like dispersion seen in the insulator, as had been suggested previously.

As we will show, residual effects from the insulator strongly distort the ”band structure” in these deeply underdoped samples resulting in a gap structure along the LDA-predicted Fermi surface that strongly deviates from the canonical d$`_{x^2y^2}`$ form. In order to reconcile this result with d-wave superconductivity, we conclude that the measured gap structure is a consequence of probing two different regimes. One regime is concerned with a small unresolved gap along a hole pocket about $`(\pi /2,\pi /2)`$. The second regime is in the vicinity of $`(\pi ,0)`$ where a large gap, which we identify as the pseudogap, is identified with the undoped insulator. This picture is consistent with the observed chemical potential shift.

The evolution of the electronic structure, however, is not solely a simple shift of the chemical potential since the spectral lineshape also evolves with doping. We further investigate here the two component structure found in the nodal direction of Ca2-xNaxCuO2Cl2. We show that this structure is not directly related to Tc, and from a MDC (Momentum Distribution Curve) analysis we find a “kink” in the dispersion which has the same energy as the ”dip” which divides the two components in the EDCs (Energy Distribution Curves). The observation of a kink near 50 meV in this system supports the observed universality of this type of behaviour in the cuprates. Finally, we discuss the similarities and differences between photoemission results from other cuprates and Na-doped CCOC.

## II Experimental

Na-doped CCOC single crystals were grown by a flux method under high pressure, by using a cubic anvil type pressure apparatus. A powder mixture of Ca2CuO2Cl2, NaCl, NaClO4, and CuO (1:0.2:0.2:0.1 molar ratio), sealed in a gold tube, was heated to 1250C and then slowly cooled to 1050C under high pressure. Uniform samples at various Na concentrations, x, were achieved by using different pressures during the synthesis. 4 GPa and 5 GPa gave a Na concentration of 10$`\%`$ and 12$`\%`$, respectively. Crystals with a maximum dimension of 1.5 $``$ 1.5 $``$ 0.1 mm3 were obtained. The Na content was estimated by comparing the c-axis lattice constant with powdered ceramic samples of known content. Magnetization measurements of Tc gave results consistent with the assigned doping concentrations. The Na concentrations used in this study were x=0.10(Tc=13) and x=0.12(Tc=22 K). Since superconductivity does not occur until x$``$0.08 in the Na-doped CCOC system, we consider our samples to lie in the heavily underdoped regime. As a reference, the optimal Tc in powder samples has been found to be 28 K. ARPES measurements were performed on beamline 5-4 of SSRL. The samples were easily cleaved in situ and measured at the indicated temperatures in the figure captions. The angular resolution was 0.29, while the energy resolution was always $``$30 meV. When examining the detailed low energy excitations (ie after figure 7) the energy resolution was improved to $``$15 meV.

## III Data and Results

### A Valence Band

We begin with the valence band spectra of Na-doped CCOC in figure 2, which illustrates a peculiarity of this system. Namely, from 8 to 13 eV binding energy there exists several pronounced features which rapidly vanish with time. The spectra from the cleave where these features were most prominent is shown in the top panel. Features are visible at 9, 10, and 12 eV which are not observed in the parent compound, CCOC. Figure 2b reveals that within 90 minutes following the initial cleave, the additional features from 8 to 11 eV have vanished. The feature at 12 eV also vanishes (not shown). For binding energies less than 7 eV no change in the spectra was observed over these time scales where particularly close attention was paid to the near $`E_F`$ features, less than 1 eV binding energy. The photon energy was adjusted to ensure the high binding energy features were not derived from a contribution of second order light from the monochromator. Comparing the spectra in panels a) and b) suggests a dispersive nature of the additional peaks not seen in the insulating compound. However, it was also observed that these features shift to lower binding energy as they age. We tentatively assign the origin of these features from 8 to 13 eV binding energy to flux inclusions, such as NaCl and NaClO4 which are used in excess during the synthesis; subsequently the surface of these flux inclusions are passivated by the photon following the cleave. Despite the uncertain nature of these high energy features, the fact that the low binding energy region of the spectra have no variation over the short time scale in which the high energy features vanish, gives us confidence that our low energy results on the Na-doped CCOC crystals are intrinsic.

Figure 3 compares the valence band spectra along the nodal direction for x=0 and x=0.1. The data on the insulator are consistent with previous reports on a similar material (Sr2CuO2Cl2), and the data for the Na-doped Ca2CuO2Cl2 were reproduced on multiple cleaves. The position of the valence band in the insulator with respect to the chemical potential can vary as much as 1 eV which is believed to be a result of pinning the chemical potential at different impurity levels that exist for the different cleaved surfaces. It is important to note that the valence band data of the insulator shown in figures 3 and 4 are at the minimum possible binding energy, since if they are shifted further towards the chemical potential, spectral weight from the tail of the valence band peaks would be pushed above the chemical potential falsely indicative of a metallic nature. The low energy excitations within 1 eV of the chemical potential, which the remainder of this paper will focus on, are barely visible on this broad scale. Generally, the two data sets are quite similar, although differing relative intensities for a few of the bands about (0,0) can be seen. The other difference is that the insulator data are sharper than those of the metal. Perhaps this is an indication that the surface of Na-doped Ca2CuO2Cl2 is not as well ordered as the pure sample grown under atmospheric conditions. Figure 4a plots the second derivative of the above spectra, which can be used to more clearly identify the dispersion of various features. This method agrees with the dispersion extracted by the peak positions of the various bands in the insulator and the metal as shown in panel 4b. The resolved bands appear to have shifted to lower binding energy with doping by roughly 300 meV, as would be expected if this were a simple band material. One apparent exception exists from the band at 5.6 eV binding energy in the metallic sample near (0,0) in figure 4b. However, due to the broadness of the features in the metal, it is likely that this band is a superposition of several peaks, resulting in this apparent discrepancy. Note that in the insulator several bands are resolved in this same region. Thus, the overall picture from the valence band data is highly suggestive of a chemical potential shift with doping. The most convincing demonstration for this is at $`(\pi ,\pi )`$ where the peaks for both concentrations are most clearly resolved. Here, the peak positions indicate that the insulator is shifted to higher energy by at least 300 meV relative to the metal. The specific value will depend on which insulating data set is used, since the chemical potential in the insulator can vary substantially as mentioned above. We again emphasize that the current binding energy of the insulator could not be found at smaller binding energies without spectral weight being shifted unrealistically above the chemical potential. Thus, neither charging of the insulator nor the variability of the chemical potential in the insulator can be responsible for the observed shift upon doping

### B Near EF Comparison

To further support this idea, figure 5 presents spectra near EF along the entire (0,0) to $`(\pi ,0)`$ and (0,0) to $`(\pi ,\pi )`$ lines, where all the spectra from this insulating data set have been shifted by 750 meV so that the centroid of the feature at $`(\pi /2,\pi /2)`$ is at the chemical potential. The value of 750 meV is specific to this particular data set, but is 300 meV more than allowed by the uncertainty in the position of the chemical potential. Along (0,0) to $`(\pi ,0)`$ the spectra are remarkably similar. In many instances fine details of the line shape match extremely well. For example, a second broad component in the electronic structure can be observed at roughly 600 meV higher binding energy from the first feature creating a heavily asymmetric lineshape. This indicates that upon doping the chemical potential simply drops to intersect the top of the valence band similar to the case of an ordinary band material. This would naturally explain the large high energy pseudogap seen at $`(\pi ,0)`$ as a remnant property of the insulator as first suggested by Laughlin.

Along (0,0) to $`(\pi ,\pi )`$ the situation is less clear. As the feature moves towards lower binding energy a sharp Fermi cutoff appears in the metallic samples, as anticipated by the large peak at the Fermi level from the energetically shifted spectra of the insulating sample. However, the match between the spectra becomes increasingly worse as $`(\pi /2,\pi /2)`$ is approached. Perhaps this simply indicates that a Fermi crossing has occurred before $`(\pi /2,\pi /2)`$, and thus the weight at $`(\pi /2,\pi /2)`$ is suppressed relative to the insulator or that some spectral weight transfer in addition to the shift of the chemical potential has occured. To see if the suppression was due to matrix element effects we also changed the polarization to maximize the nodal direction cross section, but it remained weaker. This suppression possibly indicates that increased scattering which could effectively smear out the $`k`$ resolution, and would have the greatest effect where the dispersion is the steepest (near $`(\pi /2,\pi /2)`$), has occured. However, this would require a nontrivial small angle scattering mechanism, which does not simply produce an angle independent background. Also, the cleavability of the Na-doped compounds is the same as the parent insulator, and laser reflections from the sample indicate flat surfaces, which argue against angle averaging as the cause for the observed suppression of weight. Typically, highly angle dependent valence band spectra as shown in figure 3 are indicative of a good surface, although the Na-doped valence band features are somewhat broader than in the case of the insulator. Pothuizen et al. noted that in SCOC the lineshape of the 2 eV binding energy peak at $`(\pi ,\pi )`$ is identical to the lineshape of the Zhang-Rice singlet at $`(\pi /2,\pi /2)`$. Given such a correlation a broadened low-energy spectral lineshape at $`(\pi /2,\pi /2)`$ in Na-doped CCOC is perhaps expected when the valence band at 2 eV has been broadened, as observed here (see figure 3). This suggests that the small differences between the doped system and the half filled insulator may indeed result from a difference in sample surface quality.

Considering that the features of hole doped CCOC track the dispersion of the features in the half filled insulator, one might wonder exactly how similar are the electronic states for the two doping levels. A good way to test this is by analyzing the photoemission matrix elements. Under the dipole approximation the photoemission intensity can be written as $`I|\mathrm{\Psi }_f|Ap|\mathrm{\Psi }_i|^2`$ where $`\mathrm{\Psi }_{i,f}`$ are the initial and final state wavefunctions, and $`Ap`$ is the perturbing Hamiltonian. Practically, this means that the cross section of the photoemission process will depend on the photon energy and experimental geometry as well as the wavefunctions themselves. Thus, if two wavefunctions are to be considered as the same, then they should have the same dependence on the experimental conditions. We have found that CCOC experiences a fairly dramatic change in the modulation of intensity along the nodal direction when changing the photon energy from 16.5 to 17.5 eV. In figure 6, we compare the EDCs along the nodal direction of a 10$`\%`$ Na-doped CCOC sample and an undoped CCOC sample using 16.5 eV and 17.5 eV photons. We notice that the change in the modulating intensity on going from $`E_\gamma `$=16.5 to 17.5 eV is similar for both samples. Namely, for 16.5 eV photons there is much more intensity at $`(\pi /2,\pi /2)`$ and higher $`k`$ values, which vanishes on going to 17.5 eV photons. Even the existence of the second feature at 600 meV higher binding energy from the lowest energy feature can be seen in both samples along the nodal direction. These results show that the gross features of the wavefunction of the electronic states within 2 eV of the chemical potential discussed thus far are remarkably similar. The peculiar photon energy dependence shown in figure 6 will be discussed separately with a much larger body of data. Along with the fairly rigid shift in dispersion of the low energy states and the valence band, the lineshape comparison and photon energy dependence show conclusively that the chemical potential indeed shifts upon doping the half filled insulator, Ca2CuO2Cl2, with Na.

### C Metallic Nature of Na-doped Ca2CuO2Cl2

A naive expectation of a simple chemical potential shift upon doping is the creation of a small hole pocket centered at $`(\pi /2,\pi /2)`$. However this contrasts with the majority of ARPES results from the optimal and overdoped regime which find a large Fermi surface centered about $`(\pi ,\pi )`$ as predicted by band theory. We shall thus take a closer look at the near EF electronic structure to see if the simple picture of a chemical potential shift as illustrated in figure 1b is correct. In other words, is the electronic structure unchanged with the exception of the position of the chemical potential? Figure 7 presents a sampling of the EDCs at selected $`k`$ points. The cuts shown in (a) through (d) are taken parallel to the nodal direction. For several of the cuts close to the nodal direction a broad feature starts at high binding energy and approaches the Fermi energy with increasing $`k`$, where a sharp Fermi cutoff of the spectra is indicative of a Fermi surface crossing. We will return to the presence of a peak-dip-hump like structure near the nodal Fermi crossing later. The Fermi crossings, which have been determined by the spectra with the minimum leading edge midpoint (the binding energy where the spectral intensity is half of the intensity at the peak position), have been highlighted in the figure. We note that while a minimum value of the leading edge midpoint identifies the spectra which corresponds to $`k_F`$, the value itself will be pushed back to higher binding energy if this portion of the Fermi surface is gapped. As the cuts approach $`(\pi ,0)`$ the spectral features become less pronounced and it quickly becomes impossible to identify a Fermi crossing as one can no longer clearly identify a peak position to meaningfully extract a leading edge. This is best illustrated in panel (e) where spectra along the ostensible Fermi surface predicted by band theory are presented. Near the node, there is a clear Fermi cutoff at EF, but this feature disappears quickly as one approaches the antinodal region.

One can also visualize the data with an intensity map at EF, obtained by integrating each EDC by a $`\pm `$4 meV window about EF. The highest intensity contour, shown in figure 8a, forms what appears as a small arc. The open circles indicate where a clear Fermi crossing can be observed as in the cuts shown in figure 7(a) through (c). Note this ends as the Fermi surface crosses the antiferromagnetic zone boundary. By assuming that the Fermi surface can not abruptly end in the middle of the Brillouin zone, one is left with two possibilities. Either the Fermi surface still lies along the one predicted by band theory and has been heavily gapped as one approaches $`(\pi ,0)`$ or the Fermi surface is a small hole pocket centered on or near $`(\pi /2,\pi /2)`$. We investigate these two possibilities below.

Let us first attempt to push the assumption that the large LDA-like Fermi surface holds for this heavily underdoped sample. We separate the large Fermi surface into two regions, as illustrated in figure 8b. In Region I the broad spectral function has a sharp Fermi cutoff, which clearly defines $`k_F`$, while in region II one must extrapolate the Fermi surface from the Fermi crossings found in region I assuming a large Fermi surface that does not terminate in the middle of the Brillouin zone. This is shown by the solid dots in figure 8. Note that in the region near $`(\pi ,0)`$, variations in the precise shape of the Fermi surface will not affect the following analysis as the spectra have little variation which is apparent from the EDCs in figure 7d.

Figure 9 plots the leading edge midpoint for spectra along the large Fermi surface just discussed relative to the value at $`\varphi `$=45 (the nodal position for a pure d$`_{x^2y^2}`$ functional form). In this way one can determine the gap structure along the Fermi surface. The open symbols indicate the leading edge midpoint for the Fermi crossing spectra found in region I. However, in region II the lack of a clear leading edge midpoint prevented us from identifying a Fermi crossing in the first place. Hence, we need an alternative method to describe the dispersion. Near $`(\pi ,0)`$ an energy scale (albeit a broad one) still exists in the spectra as indicated by a dot in figures 7 (d) and (e). Thus we can blindly normalize the maximum of the spectra below 300 meV to 1, and extract a leading edge midpoint via this simplified approach. This is what is shown by the filled symbols in figure 9. For a simple spectral function this procedure would give the correct gap value. In our case, one can see that close to the node this procedure roughly matches our more careful analysis, but it also allows us to extend our gap analysis beyond the regime where a clear Fermi crossing has been identified, albeit in a crude way. The solid black line is the expectation for the canonical d-wave functional form. As mentioned above, in the optimal and overdoped Bi2212 samples a simple spectral function and a large Fermi surface is observed, and the d-wave functional form well describes the experimental gap structure. However, in the heavily underdoped regime studied here, there is a strong deviation from the d-wave functional form. In order to reconcile this gap structure with d-wave superconductivity our assumption of a large Fermi surface centered at $`(\pi ,\pi )`$ must break down. In this case, figure 9 is not physically significant as it maps the energy along two physically different regions as defined in figure 8b, each with its own gap structure. In other words, by assuming a d-wave gap structure we have shown that only a fraction of the spectra used for figure 9 actually lie on the Fermi surface.

Region I does indeed represent a portion of the Fermi surface about $`(\pi /2,\pi /2)`$ which may contain a very small gap, while region II does not contain a Fermi crossing, and the resulting large gap identifies a larger energy scale. We have just ruled out the possibility that the Fermi surface in this sample is a large hole pocket centered at $`(\pi ,\pi )`$. If we maintain that the Fermi surface can not abruptly end in the middle of the Brillouin zone we must conclude that the Fermi surface forms a small hole pocket, as naively expected from the observed chemical potential shift. The question then becomes, how does such a Fermi surface develop a gap consistent with d-wave superconductivity. Given the experimental resolution we can not comment on whether or not the flat region in figure 9 would itself fit a $`|cos(k_xa)cos(k_ya)|`$ functional form, which is the simplest expectation for d-wave pairing, although models exist with small hole pockets and d-wave pairing for which this is not the case. We note that for a superconductor with a transition temperature of 13 K, the maximum of the gap from BCS theory is expected to be less than 2 meV at the antinode and even smaller along the observed arc, which is too small for our experimental resolution and is not the objective of this study. Perhaps this gap could be measured with an improved spectrometer.

By symmeterizing the hole pocket about the antiferromagnetic zone boundary one finds an area of 5.7$`\%`$, which implies a doping level of 11.4$`\%`$ due to the electron spin. This is in rough agreement with the hole doping of 10$`\%`$ determined by comparing Tc and lattice constants with results on powdered samples. Again we note that the hole pocket is also not necessarily symmetric about the antiferromagnetic zone boundary, which could substantially modify these values. The idea of a small hole pocket about $`(\pi /2,\pi /2)`$ in Na-doped CCOC is also consistent with the observations of shadow bands at higher binding energy. Presumably, the coherence factors responsible for the shadow band are very weak making it difficult to distinguish a second Fermi crossing from the background at EF. We stress that while our findings are quite consistent with the small Fermi pocket picture, we did not see the shadow Fermi surface at zero energy even though we do observe a shadow band at higher energies.

Meanwhile, the observation of a chemical potential shift in this material from the insulating state and the photon energy dependence of the lineshape emphasize that the features near $`(\pi ,0)`$ (region II) in the Na-doped Ca2CuO2Cl2 sample are a direct result of the magnetism which dominates the antiferromagnetic insulator. To illustrate that the energy scale near $`(\pi ,0)`$ is the same as that of the insulator, recall figure 5 where the spectra of the insulator matches that of the Na-doped sample at $`(\pi ,0)`$.

We have also repeated the above gap structure analysis on a sample with x=0.12 (Tc=22 K), and some of the results are presented in blue in figure 8 and as circles in figure 9. The data reproduce the above analysis with a few small quantitative differences. The Fermi arc increases in size (shown in figure 8b), consistent with an increase in doping, while the binding energy near $`(\pi ,0)`$ reduces slightly which demonstrates that this energy scale is decreasing with increasing doping. This is also expected as we have identified this energy scale with the magnetism of the insulator whose strength should diminish together with the antiferromagnetic correlations as holes are added to the system. This has also been observed in previous studies on the high energy pseudogap. We should note that the hidden message of describing the Fermi surface in the underdoped regime as a small pocket about $`(\pi /2,\pi /2)`$ is that the system must at some point transform upon doping to the large Fermi surface centered about $`(\pi ,\pi )`$ observed in the overdoped regime. In measuring the volume of the pocket for the x=0.12 sample we find the doping level to be 21$`\%`$. The fact that this analysis begins to show a deviation from a doping level of x perhaps begins to identify a transition to a state where the doping level should be described as 1-x (the expected volume of the large Fermi surface centered at $`(\pi ,\pi )`$).

In comparison with other cuprates, the lack of a Fermi cutoff in the EDC’s would normally be identified as an extremely large pseudogap. In an attempt to quantify this suppression of weight, we compare 10$`\%`$ Na-doped CCOC to 10$`\%`$ Dy-doped Bi2212 in figure 10. 10$`\%`$ Dy-doped Bi2212 is underdoped with a Tc of 65 K giving it a Tc/Max(Tc) of 0.72, while 10$`\%`$ Na-doped CCOC has a ratio of 0.46. A comparison with 17.5$`\%`$ Dy-doped Bi2212 (not shown) with a Tc/Max(Tc) ratio of 0.28 yielded a similar comparison. One clearly observes the loss of spectral weight at the Fermi energy for the Fermi momentum of $`(\pi ,0.3\pi )`$ when compared to the respective Fermi crossings along the nodal direction near $`(0.4\pi ,0.4\pi )`$. There are two common methods for characterizing the pseudogap with ARPES. The first, is by the shift of the leading edge midpoint of a spectra to higher binding energy relative to the chemical potential as indicated by the dashed lines. This is typically referred to as the low energy pseudogap. For Dy-doped Bi2212 this is reasonably well defined, and gives a value of 30 meV. For the Na-doped CCOC sample, we find a value of roughly 50 meV at $`(\pi ,0.3\pi )`$. However, as previously noted, the lack of a Fermi cutoff near $`(\pi ,0)`$ implies that the low-energy pseudogap is not well defined in the Na-doped CCOC spectra. As a result, we are inclined to characterize the pseudogap by what is known as the high energy pseudogap which identifies a larger energy scale in the spectra about $`(\pi ,0)`$, and is sometimes referred to as a “hump”. Here, we find a high energy pseudogap of roughly 300 meV for Na-doped CCOC, compared with 120 meV for Dy-doped Bi2212 as indicated by the triangles in figure 10. Independent of the method chosen, the pseudogap for Na-doped CCOC appears significantly larger than that of Bi2212 at a comparable doping level.

The high energy pseudogap seen in Bi2212 was first conjectured by Laughlin to be a result of the $`d`$-wave-like modulation of the dispersion of the insulating oxy-halides. From the observed shift of the chemical potential it is clear that the high energy pseudogap indeed directly results from the dispersion seen in the parent compound, although it appears that the energy scales of the high energy pseudogap have a small system dependent variation. The implication for such a large pseudogap relative to the Bi2212 system at comparable doping should be investigated, and the bilayer effect in Bi2212 will need to be taken into consideration. If the condensation energy for superconductivity is acquired by the gapping of low energy states near $`(\pi ,0)`$, then this would result in superconductivity being less favorable as evidenced by the smaller optimal Tc(28 K and 90 K for Na-doped CCOC and Bi2212, respectively).

### D Nodal Peak-Dip-Hump Structure

A two component structure can be seen in the low energy spectra of Na-doped CCOC at the Fermi crossing in the nodal direction. This was shown in figure 7e, where we plotted EDCs along the Fermi arc and then towards $`(\pi ,0)`$ as indicated in the cartoon inset. The reason for the difference in clarity of this structure between cleaves is unknown, but is likely due to a variation of inhomogeneities and experimental conditions which make it impossible to resolve the two components that appear as one in the $`(0.43\pi ,0.43\pi )`$ spectra of figure 5. Note that in general, the presence of sharp low energy features is a tribute to a high sample quality. The loss of the two component structure as one approaches $`(\pi ,0)`$ is indicative of the observed large pseudogap, which precludes the identification of any similar structure in the low energy excitations near $`(\pi ,0)`$.

The observation of a two component structure in these superconducting samples begs the question of how it may be related to superconductivity. Figure 11 presents temperature dependence of this feature from 20 K($``$ Tc) up to 120 K. It is clear that the peak-dip-hump like structure is observed up to 75 K, and is finally smeared out at 120 K. By cycling the temperature back down to 20 K, we ensure that the broadening at high temperatures is not an aging effect. Since even optimally doped polycrystaline samples have only a maximum Tc of 28 K, it is clear that the two component structure observed here does not turn on near Tc. This is in contrast with the more familiar two component structure in ARPES on the Bi2212 system which has now also been seen in YBCO and Bi2223, and does turn on with Tc. It is important to note that, aside from the different temperature dependence, the two component structure observed in Na-doped CCOC occurs along the nodal direction as opposed to the $`(\pi ,0)`$ region as in the other cases. Furthermore, we note that CCOC is a single layer material and thus the peak-dip-hump seen here can not be a consequence of bilayer splitting as has been conjectured for the peak-dip-hump structure observed in Bi2212.

In order to gain some insight into the origin of this structure figure 12 presents a MDC analysis along the nodal direction of a x=0.10 sample. The dashed curves are the single Lorentzian fits to selected MDCs. Similar procedure in other cuprates has been used to extract self energies from the fit parameters. Here the fits serve as a method for parameterizing the observed dispersion. Figures 12b and 12c plot the peak positions and widths as a function of the binding energy of the corresponding MDCs. A change of slope in the dispersion of the peak positions is seen near 55 meV. This value is attained from the intersection of two linear fits ranging from 0 to 30 meV and 100 to 200 meV. We shall refer to this feature as a “kink”. A similar kink in other cuprate systems has received significant attention in recent ARPES literature. Certainly, the similarity of the kink energy position to the dip energy position seen in the EDCs suggests that these two results are related. The MDC width also displays an accelerated decrease in width as indicated by the dotted lines in figure 12c.

What is the possible origin of a two component structure in the EDC near $`k_F`$ along the nodal direction? Pure macroscopic phase separation of metallic and undoped insulating domains can be ruled out, since the feature seen in the insulator would not have shifted to the chemical potential in this case. One possibility is that a sharp, coherent-like peak is growing with doping, while the incoherent features at half filling slowly vanish. This idea is captured in the recent phase string calculations done by Muthukumar, Weng, and Sheng. In this sense, the two components would hint at a balance between the antiferromagnetic insulator and the drive for the system to become metallic. An alternative scenario is that the peak-dip-hump structure is the result of coupling to a collective mode. This could also naturally explain the presence of the kink. Although still a debated interpretation, Lanzara et al. propose that the universally observed kink is due to electron-phonon coupling. In this regard, it is interesting to note that the phonon breathing mode is expected to be roughly 10 meV less for CCOC on the basis of its a-axis lattice constant than in the case of insulating Bi2212 or LSCO as determined by optics, while the kink energy position determined by ARPES is less by roughly the same amount.

Distinguishing between these two scenarios will be difficult. A problem with the first picture would be that the kink would have to be explained as an artifact of the MDC analysis due to the multiple electronic features which are not individually resolved. This is not appealing considering the apparent universality of the kink feature among cuprates, and one would then favor a coupling to a collective mode. However, it is surprising that the kink is seen in both LSCO where the states were created inside the gap, and in Na-doped CCOC where the states appear to originate from the effective lower Hubbard band which has now been shifted to the chemical potential. This suggests that a third possibility, combining both ideas, may also be considered for the doping evolution. Namely, that the chemical potential has shifted upon doping, the precise location would be system dependent, while states are created at EF which then couple to some collective excitations creating the universally observed kink (see figure 1d). If this is the case, the Ca2-xNaxCuO2Cl2 and La2-xSrxCuO4 system would no longer appear so different. Unfortunately, the present data can not distinguish between these three possibilities, and we must leave this as an opportunity for future experiments with more available doping levels to investigate.

## IV Conclusions and Discussion

The data from the valence band, near EF spectra, and photon energy dependence in the Ca2-xNaxCuO2Cl2 system provide convincing evidence that the chemical potential shifts to the top of the valence band upon doping the insulator, Ca2CuO2Cl2. Surprisingly, these results appear to present a different evolution across the metal to insulator transition than in LSCO where the chemical potential is pinned and states are created inside the gap upon doping.

With the comparison of CCOC we have shown that the high energy pseudogap is indeed a remnant property of the insulator. Meanwhile, by examining the low lying excitations we found that an attempt to use a large Fermi surface centered at $`(\pi ,\pi )`$ results in a gap structure which deviates from the canonical d-wave functional form. Instead a small Fermi pocket about $`(\pi /2,\pi /2)`$ as naively expected from a rigid chemical potential shift is a better description of heavily underdoped Ca2-xNaxCuO2Cl2. This is the case when we insist that the Fermi surface can not terminate at a point in the middle of the Brillouin zone. A Fermi surface arc could also explain the data if this were allowable in a non-Fermi liquid picture. We anticipate this region to also possess a gap structure which could differentiate between competing models in the underdoped regime. Finally, a two component structure is also observed in the EDCs along the nodal direction. This appears to be related to a kink seen in the MDC derived dispersion, which is now a universally observed feature of the cuprates.

The presence of features seen in the insulator within the low energy spectra of the Na-doped superconductor suggests a strong interplay between antiferromagnetism and superconductivity in the underdoped regime. One way to think of these results is to start from the overdoped regime where superconductivity is the first instability encountered. Upon further underdoping, magnetism begins to play a role which opens a large pseudogap, and suppresses superconductivity. Alternatively, one can start with the insulator, which forms a small pocket as carriers are introduced into the system. Superconductivity then creates a new instability that opens a small gap along the Fermi surface, which we are unable to detect given our resolution. We look forward to future experiments which will clarify Ca2-xNaxCuO2Cl2’s evolution to the more ”conventional” behaviour of the overdoped regime, and its relation to the other cuprate superconductors.

## V Acknowledgments

F.R. would like to thank M.B. Walker and N. Schwanger for generously providing assitance to fascilitate completion of this work, as well as A.-M.S. Tremblay and S.C. Zhang for helpful discussions. This research was carried out at the Stanford Synchrotron Radiation Laboratory which is operated by the D.O.E. Office of Basic Energy Science, Division of Chemical Sciences. The Office’s Division of Material Science provided funding for this research. The Stanford work is also supported by NSF grant DMR-0071897.