# Random Matrices and the Convergence of Partition Function Zeros in Finite Density QCD

## I Introduction

In contrast with the numerous successes of lattice QCD, simulations at finite chemical potential  have oscillated between mildly promising and outright frustrating. The source of the trouble lies in the following. In the Euclidean formulation of QCD, the chemical potential spoils the anti-Hermiticity of the Dirac operator. As a result, the fermion determinant is no longer a real number. In general, it has a complex phase. Hence, the action cannot serve as a statistical weight in a Monte Carlo sampling of field configurations.

Quenched simulations, where the fermion determinant is not included in the statistical weight, may provide a fairly reliable approximation to selected observables of the true unquenched theory. However, in the presence of a chemical potential, quenched simulations have produced consistently unphysical results . The reason is that the quenched theory is the $`N_f0`$ limit of an unphysical theory where the fermion determinant is replaced by its absolute value . This is a theory with a second, “conjugate” set of anti-quark species together with the normal quarks. Because of Goldstone bosons consisting of a quark and a conjugate anti-quark, the critical chemical potential in quenched QCD is half the pion mass. This phenomenon was demonstrated in lattice simulations by Gocksch  using a $`U(1)`$ toy model and was understood analytically in  by using a random matrix model inspired by QCD.

One important conclusion of further studies of the same RMM is that the phase of the fermion determinant leads to very large cancellations in the ensemble averaging . A measure of this phenomenon is given by the fact that the partition function is proportional to $`\mathrm{exp}(\mu ^2N)`$, where $`N`$ is the size of the random matrix, corresponding to the number of sites in a lattice simulation. Cancellations of this magnitude would require prohibitive statistics in order for a brute force simulation<sup>*</sup><sup>*</sup>*where the determinant is included as an observable in a quenched ensemble to be successful. We wish to note that in some models the construction of clever algorithms makes it possible to deal with these cancellations in Monte Carlo simulations .

It would be hard to overrate the potential importance of a successful lattice approach to QCD at finite chemical potential. At this time, there is not even an estimate for the value of the critical chemical potential from lattice simulations. In the absence of the guidance provided by the lattice, it is difficult to assess the many semi-empirical descriptions of nuclear matter at high density . For instance, an extension of the RMM to include temperature via the first Matsubara frequency  gives reasonable predictions about the phase diagram of QCD, even while ignoring most of its dynamics. In the present paper we wish to exploit the qualitative similarity between our simple RMM and $`N_c=3`$ QCD at finite chemical potential in an attempt to understand certain lattice results on the problem of finite $`\mu `$.

We are interested in the analytic dependence of the QCD partition function on the chemical potential $`\mu `$. This can be obtained by computing the coefficients of the expansion of the partition function in powers of the chemical potential or the fugacity. The Glasgow method of lattice QCD  is designed to do this. The unquenched partition function can be seen as the quenched average of the fermion determinant, i.e. averaged only with the gauge action. In general, one may also use an unquenched ensemble at some fixed chemical potential $`\mu =\mu _0`$, and include the inverse fermion determinant at that same value,

$`Z(\mu )=\mathrm{det}\text{ /}D(\mu )_{\mathrm{gauge}}=C{\displaystyle \frac{\mathrm{det}\text{ /}D(\mu )}{\mathrm{det}\text{ /}D(\mu _0)}}_{\mathrm{gauge},\mu =\mu _0}.`$ (1)   
Here $`C`$ is an irrelevant constant. One expects the efficiency of the averaging process to depend on the overlap between the quantity being averaged and the distribution used to generate the ensemble. When these two functions have their largest values in vastly different places, this is known as an overlap problem.

In the lattice Glasgow method the fermion determinant is expanded in powers of the fugacity $`\xi =\mathrm{exp}(\mu )`$. The expansion is finite and exact, since the fermion determinant is just an $`N\times N`$ matrix ($`N`$ is the number of lattice points times the number of degrees of freedom per site). It is obtained by writing $`det\text{ /}D\left(\mu \right)=\xi ^Ndet\left(P+\xi \right)`$, where $`P`$ is called the propagator matrix. The expansion coefficients, written in terms of the eigenvalues of $`P`$, are then obtained by ensemble averaging. The zeros of the partition function in the complex $`\mu `$ plane map out the phase structure. In particular, the ones close to the real axis define the critical value(s) of $`\mu `$. Since the original paper by Yang and Lee  this approach has been widely used in statistical mechanics. For some recent applications we refer to .

With the Glasgow method one obtains information about the full $`\mu `$-dependence of $`Z`$, using an ensemble generated at one fixed value of the chemical potential. Of course, the question of the overlap between whatever ensemble one is using and the fermion determinant for the given $`\mu `$ remains. The issue is even more ominous considering that we are unable to even define the notion of an ensemble at nonzero $`\mu `$, even in this random matrix model, due to the complex action.

In this paper, we analyze the Glasgow method using a random matrix model at nonzero chemical potential. As we have already discussed before, this model mimics the problems of the QCD partition function at nonzero chemical potential. Since the phase structure of this model is known analytically, it is an ideal testing ground for evaluating this algorithm and obtaining a better understanding of its problems.

In section 2 we introduce the random matrix model and derive some of its analytical properties. The bulk of this paper is in section 3. It contains our numerical analysis of the Glasgow method, and an explanation is given of the poor convergence of certain zeros of the partition function. Concluding remarks are made in section 4.

## II Random matrix model

We consider a random matrix model (RMM) defined by the partition function 

$`Z(m,\mu )`$ $`=`$ $`{\displaystyle 𝒟Ce^{N\mathrm{tr}CC^{}}det\left(D(m,\mu )\right)},`$ (2)   
$`D(m,\mu )`$ $`=`$ $`\left[\begin{array}{cc}m& iC+\mu \\ iC^{}+\mu & m\end{array}\right].`$ (5)   
Here, $`C`$ is an $`N\times N`$ complex matrix, and the integration is over the Haar measure,

$`{\displaystyle 𝒟C}={\displaystyle \underset{ij=1}{\overset{N}{}}}{\displaystyle 𝑑C_{ij}𝑑C_{ij}^{}}.`$ (6)   
This model was first formulated  for $`\mu =0`$ in order to describe the correlations of the smallest eigenvalues of the Dirac operator. In this case it has been shown rigorously that the model describes the zero-momentum sector of the low-energy effective partition function of QCD .

In the present context, the matrix $`D(m,\mu )`$ mimics the QCD Dirac operator for quark mass $`m`$ and chemical potential $`\mu `$. The integration over matrix elements replaces the integration over gauge field configurations. The massless part of our random Dirac operator is anti-Hermitian for $`\mu =0`$, but for $`\mu 0`$ it has no definite Hermiticity properties. In order to study the properties of the partition function in the chemical potential plane $`\mu `$, we rewrite the fermion determinant as follows,

$`det\left(D(m,\mu )\right)=det\left(\begin{array}{cc}iC+\mu & m\\ m& iC^{}+\mu \end{array}\right)=det\left(F(m)+\mu \mathrm{𝟏}\right).`$ (9)   
The matrix $`F(m)`$ is analogous to the propagator matrix from lattice QCD  in the sense that its eigenvalues are the values of $`\mu `$ for which the fermion determinant vanishes. In terms of the eigenvalues $`\lambda _k`$ of $`F(m)`$ the RMM partition function is

$`Z(m,\mu )={\displaystyle 𝒟Ce^{N\mathrm{tr}CC^{}}\underset{k}{}(\lambda _k+\mu )}.`$ (10)   
The quantity $`n=_\mu \mathrm{ln}Z(m,\mu )/N`$ is the analog of the baryon number density of QCD. It is equal to the ensemble average of $`\frac{1}{N}_k\frac{1}{\lambda _k+\mu }`$. In QCD at zero temperature, one expects the baryon number density to be identically zero for small $`\mu `$, and then to increase starting from a certain critical value of $`\mu `$ . Similarly, our model shows a phase transition with an increase of the baryon number density. However, $`n`$ is not zero below the critical value of $`\mu `$. See  for an explanation of the relationship between $`n`$ and the baryon number density of QCD.

It is not clear how to define a statistical ensemble of gauge field configurations (or random matrices for that matter), corresponding to the true partition function with finite $`\mu `$. In the quenched approximation one discards the fermion determinant, so the partition function does not depend any more on $`\mu `$ or $`m`$. However, a ‘number density’ can still be computed by taking the average of $`\frac{1}{N}_k\frac{1}{\lambda _k+\mu }`$. Similarly, one can define a quenched ‘chiral condensate’ from the ensemble average of $`\frac{1}{N}\mathrm{Tr}\left(D(m,\mu )^1\right)`$. The quenched approximation can be interpreted either as the limit of a process where one takes the number of flavors to zero, or the quark mass $`m`$ to infinity. The quantity

$`n(\mu _0,\mu )={\displaystyle 𝒟Ce^{N\mathrm{tr}CC^{}}\underset{k}{}(\lambda _k+\mu _0)\frac{1}{N}\underset{k}{}\frac{1}{\lambda _k+\mu }}`$ (11)   
is called the partially quenched baryon number density. The “sea” $`\mu _0`$ defines the ensemble, and the “valence” $`\mu `$ probes the eigenvalue distribution.

### A Quenched eigenvalue distribution

For $`m=0`$ the propagator matrix is block-diagonal, and its eigenvalues are $`i`$ times those of $`C`$. The exact distribution of the eigenvalues of a general complex matrix has been calculated a long time ago by Ginibre . In our normalization, it is given by

$`\rho (\lambda _1,\mathrm{},\lambda _N)=𝒞{\displaystyle d^2\lambda _1\mathrm{}d^2\lambda _Ne^{N\underset{k}{}\left|\lambda _k\right|^2}\underset{k>l}{}\left|\lambda _k\lambda _l\right|^2}`$ (12)   
The corresponding one-point function is

$`\rho (\lambda )=𝒞e^{N|\lambda |^2}{\displaystyle \underset{k=0}{\overset{N1}{}}}{\displaystyle \frac{\left(N|\lambda |^2\right)^k}{k!}}`$ (13)   
In the large $`N`$ limit the eigenvalues are uniformly distributed in the complex unit circle. This follows from the properties of the truncated exponential, to be discussed in more detail later. For $`|\lambda |<1`$, the truncated exponential is a good approximation of the complete one, so $`\rho `$ is a constant. For $`|\lambda |>1`$, the truncated exponential behaves like a power of $`|\lambda |^2`$ which is quickly suppressed by $`\mathrm{exp}(N|\lambda |^2)`$, so the distribution vanishes with a sharp tail of width of order $`1/\sqrt{N}`$.

For nonzero $`m`$ we can calculate the eigenvalue distribution of the propagator matrix in the large $`N`$ limit using the conjugate replica trick . We consider the partition function, where we have replaced the fermion determinant with its absolute value squared,

$`Z(m,m^{},\mu ,\mu ^{})`$ $`=`$ $`{\displaystyle 𝒟Ce^{N\mathrm{tr}CC^{}}\left|det\left(\begin{array}{cc}m& iC+\mu \\ iC^{}+\mu & m\end{array}\right)\right|^{2\stackrel{~}{n}_f}}.`$ (16)   
Here we can use either $`m`$ or $`\mu `$ as a dummy variable probing the eigenvalue distribution of the corresponding operator as a function of the other variable which is made real and therefore has a physical meaning. For the spectrum of the propagator matrix we set $`m=m^{}`$ and probe the eigenvalues with the complex dummy variable $`\mu `$. This is the reverse of what was done in  where the spectrum of the Dirac operator was investigated for given $`\mu `$. The conjugate replica trick also allows us to calculate the eigenvalue distribution in the present case.

Because of the absolute value of the determinant under the integral, the partition function is expected to be a smooth function of $`\stackrel{~}{n}_f`$ and the limit $`\stackrel{~}{n}_f0`$ should be obtainable from the partition function for positive integral values of $`\stackrel{~}{n}_f`$. For any $`\stackrel{~}{n}_f`$, the eigenvalue density is positive definite so it cannot be an analytic function of a complex variable. The resolvent is defined as $`G(z)=\frac{1}{N}_k\frac{1}{\lambda _kz}`$. The average eigenvalue density is then given by $`NRe\left(_z^{}G(z)\right)`$. For the partition function (16) the resolvent in the $`\mu `$ plane is $`G(z=\mu )=_\mu Z(m,\mu ,\mu ^{})`$.

By standard manipulations  the partition function can be rewritten as an integral over an $`2\stackrel{~}{n}_f\times 2\stackrel{~}{n}_f`$ complex matrix. For $`\stackrel{~}{n}_f=1`$, one finds

$`Z(m,m^{},\mu ,\mu ^{})={\displaystyle d^2ad^2bd^2cd^2𝑑e^{N(aa^{}+bb^{}+cc^{}+dd^{})}det\left[\begin{array}{cccc}m+a& \mu & 0& id\\ \mu & m+a^{}& ic& 0\\ 0& id^{}& m^{}+b^{}& \mu ^{}\\ ic^{}& 0& \mu ^{}& m^{}+b\end{array}\right]}.`$ (21)   
The resulting saddle point equations have two kinds of nontrivial solutions, depending on whether the off-diagonal quantities $`c,c^{},d,d^{}`$ are identically zero or not. If they vanish, the partition function factorizes into pieces that depend only on $`m,\mu `$ or on $`m^{},\mu ^{}`$. Then, $`G(z)`$ is always an analytic function of $`z`$. Therefore, the region in the complex plane of $`z`$ ($`m`$ or $`\mu `$) where the eigenvalues are located must be dominated by the other kind of solution, which mixes the parameters and their conjugates. It turns out that this solution is such that $`c=c^{}=d=d^{}`$, and the quantity $`cc^{}`$ is positive. For $`m=m^{}`$ real, it is given by

$`cc^{}={\displaystyle \frac{x^2}{x^2m^2}}{\displaystyle \frac{m^2}{4(x^2m^2)}}+{\displaystyle \frac{y^2m^2}{2x^2(x^2m^2)}}{\displaystyle \frac{(x^2+y^2)y^2m^2}{4x^2}}(x^2+y^2).`$ (22)   
The boundary of the domain of eigenvalues is given by the curve, where the two types of solutions merge, i.e. the condition $`cc^{}=0`$. One can see immediately that for $`m=0`$ the boundary $`cc^{}=0`$ reduces to the unit circle.

In Fig. 1 we show the distribution of the eigenvalues in the complex plane of the propagator matrix of size $`N=96`$ for masses of value $`m=0`$ and $`m=0.0625`$. Each plot consists of eigenvalues from 20 independent configurations. The analytic curve given by $`cc^{}=0`$ is also drawn in each case. We clearly see that this does give the correct result for the boundary of eigenvalues.

### B Unquenched partition function

The unquenched partition function defined in eq. (5) can be computed analytically . It is given by

$`Z(m,\mu )={\displaystyle 𝑑\sigma 𝑑\sigma ^{}\mathrm{e}^{N\mathrm{tr}\sigma \sigma ^{}}det\left(\begin{array}{cc}m+\sigma & \mu \\ \mu & m+\sigma ^{}\end{array}\right)^N}`$ (25)   
where the integration is over the $`N_f\times N_f`$ matrix $`\sigma `$. In the $`N\mathrm{}`$ limit the integrals can be evaluated via saddle point approximation. The saddle points are given by a cubic equation,

$`\sigma ^{}`$ $`=`$ $`\sigma `$ (26)   
$`\sigma (m+\sigma )^2\mu ^2\sigma `$ $`=`$ $`m+\sigma ,`$ (27)   
where the matrix variable is diagonal and $`\sigma `$ is now a number. For fixed real $`m`$, there are four branch points in the complex $`\mu `$ plane, given by four of the six $`\mu `$ values for which the discriminant of the above equation,

$`D={\displaystyle \frac{1}{27}}\left(m^4\mu ^2m^2\left(2\mu ^45\mu ^2{\displaystyle \frac{1}{4}}\right)+\left(1+\mu ^2\right)^3\right),`$ (28)   
vanishes. The branch points are connected by two branch cuts. The derivatives of the partition function are discontinuous across these cuts. The points where the cuts cross the real axis are interpreted as critical values of the parameters. They indicate a first-order phase transition in the thermodynamic limit .

For finite $`N`$, the partition function can be evaluated as a polynomial in either $`m`$ or $`\mu `$. The coefficients are obtained by expanding the determinant and performing the integrals. One exact expansion of the partition function is given by 

$`Z(m,\mu )={\displaystyle \frac{\pi N!}{N^{N+1}}}{\displaystyle \underset{k=0}{\overset{N}{}}}{\displaystyle \underset{j=0}{\overset{Nk}{}}}{\displaystyle \frac{(Nm^2)^k}{(k!)^2}}{\displaystyle \frac{(N\mu ^2)^j}{j!}}{\displaystyle \frac{(Nj)!}{(Njk)!}}.`$ (29)   
The zeros of this polynomial can be readily calculated. They are located along lines in the complex plane and in the limit $`N\mathrm{}`$ they converge to the exact cuts found by the saddle point analysis below (33).

### C Phase structure and zeros of a partition function

When a partition function has a nontrivial phase structure in the thermodynamic limit, the complex plane of one of its thermodynamic parameters is split into regions separated by cuts. Inside each region, the partition function is analytic in the parameters, so that its derivatives are continuous. The (first-order) transitions occur across the cuts, where the (logarithmic) derivatives of the partition function will in general have discontinuities.

Before taking the thermodynamic limit (i.e. at finite $`N`$), the partition function is an analytic function of the parameters. If this partition function at finite $`N`$ is a polynomial in any of its thermodynamic parameters, such as in our case the mass or the chemical potential, it will have zeros in the complex plane of these parameters. Its logarithmic derivatives will have poles at the locations of the zeros. In the thermodynamic limit, the zeros coalesce into the same cuts that define the phase structure.

Our partition function illustrates nicely the connection between the zeros and the cuts. For $`m=0`$ it is a truncated exponential,

$`Z_N(0,\mu )={\displaystyle \frac{\pi N!}{N^{N+1}}}{\displaystyle \underset{j=0}{\overset{N}{}}}{\displaystyle \frac{(N\mu ^2)^j}{j!}}.`$ (30)   
For $`\mu ^21`$ the largest term in this sum occurs well before the truncation so we have a good approximation of the exponential. On the other hand, for $`\mu ^21`$ the series is dominated by the term with $`j=N`$. One can obtain a better estimate using the incomplete gamma function, which is closely related to the truncated exponential :

$`n!e^x{\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \frac{x^k}{k!}}`$ $`=`$ $`\mathrm{\Gamma }(1+n,x)={\displaystyle \underset{x}{\overset{\mathrm{}}{}}}e^tt^n𝑑t.`$ (31)   
For real positive arguments the situation is simple. If $`x`$ is less than the value $`t=n`$ which maximizes the integrand, then the saddle point is integrated over, the integral is a good approximation to the gamma function and the exponential is obtained. For $`x>n`$ the integral is dominated by its endpoint and the integrand is well approximated by the last term of the series.

If $`x`$ is complex, the saddle point dominates if the integration contour can be deformed across the saddle point. For our partition function

$`Z_N(0,\mu )`$ $`=`$ $`\pi e^{\mu ^2N}{\displaystyle \underset{\mu ^2}{\overset{\mathrm{}}{}}}e^{Nu}u^N𝑑u,`$ (32)   
this happens if $`\mathrm{Re}(\mu ^2+\mathrm{ln}\mu ^2)<1`$ and $`\mathrm{Re}\mu ^2>1`$. In this case $`Z_N(0,\mu )\mathrm{exp}(\mu ^2NN)`$. If the integral is dominated by the endpoint then $`Z_N(0,\mu )\mu ^{2N}`$. Nontrivial zeros of the partition function are obtained for values of the chemical potential where the saddle point contribution and the end point contribution are of comparable order of magnitude. The zeros are thus given by the equation 

$`\mathrm{Re}(1+\mu ^2+\mathrm{log}\mu ^2)=0.`$ (33)   
At finite $`N`$ the deviation of the zeros from this critical curve is of order $`1/\sqrt{N}`$.

## III Glasgow averaging in the RMM

The unquenched partition function for given $`\mu `$ and $`m`$ can be thought of as the quenched expectation value of the fermion determinant. If the zeros of the partition function $`\xi _k`$ are known, we can express the partition function as a product over its zeros

$`{\displaystyle \underset{k}{}}\left(\lambda _k+\mu \right)={\displaystyle \underset{k}{}}\left(\xi _k+\mu \right).`$ (34)   
In this identity, the zeros of the partition function appear as a kind of “averaged eigenvalues” of the propagator matrix. The Glasgow method from lattice QCD attempts to perform this averaging. For a given matrix from the quenched ensemble, one writes out the coefficients of the polynomial on the left-hand side. The average coefficients are the coefficients of the right-hand side. The main focus of this paper is to study the convergence properties of the zeros to the exact ones for finite ensembles.

The eigenvalues are in general complex. Because of the structure of $`F(m)`$, the eigenvalues occur in pairs $`\{\lambda _k,\lambda _k^{}\}`$. Also the matrix $`C`$ occurs with the same probability as $`C`$ in (5). Therefore the ensemble average will also contain the pair $`\{\lambda _k,\lambda _k^{}\}`$. Upon ensemble averaging, one can then easily show that the odd coefficients vanish and the remaining coefficients are real.

These simplifications should not mislead us into believing that we have safely avoided the trouble of averaging over the complex phase of the determinant. As we will see shortly, the problem will show up in the form of a very high precision needed for the coefficients in order to calculate the roots reliably. Since the suppression achieved by averaging over the phase of the determinant is on the order of the magnitude of the unquenched partition function, which for $`m=0`$ is $`\mathrm{exp}\left(N(1+\mu ^2)\right)`$, it becomes exponentially difficult to achieve such precision.

The quenched ensemble is not the only way to sample the set of all matrices. One may multiply and divide by any convenient function in the above formula. One factor modifies the ensemble, i.e., the way the individual matrices are generated, and the other one is used to compensate each configuration for the modified weight. One obvious choice is to use the unquenched ensemble at $`\mu =0`$.

In the next section we report the results of numerical experiments performed using Glasgow averaging for the partition function (5), both with the quenched ensemble and the unquenched ensemble at $`\mu =0`$.

### A Numerical simulations

We performed two different kinds of simulations. For $`m0`$, we generated matrix elements with the Gaussian distribution given by (5). We constructed the propagator matrix and obtained its eigenvalues. For each set, the coefficients of the corresponding polynomial are calculated and added to the average. For $`m=0`$, a much more economical procedure is possible, namely, generating sets of eigenvalues directly, using the exact Ginibre distribution (12). In this case we have employed a Metropolis algorithm. A given set is varied using small steps and the modified set is accepted depending on the corresponding value of the weight function. We have found that both cases have similar convergence properties, and in this paper we will only report on the case $`m=0`$.

In Fig. 2 we show the zeros in the complex $`\mu `$ plane for a quenched ensemble of matrices of size $`N=32`$ averaged over $`N_E=10^6,10^7`$ configurations. The averaged roots do converge to the exact values as expected, but the convergence is extremely slow. The unconverged zeros are the ones situated closer to the real axis. They are situated in a cloud of a well defined shape and most of them are located on its edge. As the averaging proceeds the cloud shrinks and finally disappears. All the roots situated outside the cloud are obtained correctly for a given $`N`$, $`N_E`$ combination. This is illustrated in Fig. 3 where we perform ensemble averaging as large as $`N_E=10^8`$. In order to obtain more converged roots, we had to reduce the size of the matrices to $`N=16`$. In this figure we also show results for $`N_E=10^4`$ and $`10^6`$. The roots next to the real axis are always the last to converge. This is unfortunate since the value of the critical chemical potential (in our case $`\mu _c=0.527`$) is determined precisely by discontinuity across the real axis. The same pattern is observed for various matrix sizes $`N`$. The shape of the cloud of unconverged zeros is very similar. However, the number of configurations $`N_E`$ corresponding to a given degree of convergence increases sharply with the matrix size $`N`$. One is able to determine $`\mu _c`$ with reasonable accuracy only for small values of $`N`$. As we can see in Fig. 3, even in that case a very large number of configurations is required. Only after averaging over $`10^8`$ configurations do we get a good estimate of $`\mu _c`$, but the roots are still not completely converged.

Convergence is not improved by including the determinant at $`\mu =0`$ in the statistical weight used to generate the eigenvalues. It was suggested previously  that using an unquenched ensemble at zero chemical potential might improve the efficiency of the averaging. Our results indicate that, at best, using the unquenched ensemble has no effect, but as it appears from Fig. 3 the results are actually worse than for the quenched ensemble. An easy explanation is that a small number of configurations with the determinant close to zero are assigned a very large weight factor (cf. (1)) which has a destructive effect. Using a “negative” number of flavors is another possibility but it is not likely to improve convergence.

We were able to run simulations which show some noticeable degree of convergence with $`N`$ up to $`48`$. There are many ways one could measure the degree of convergence of a given set of roots. The measure we will employ is the distance $`\mathrm{\Delta }(N_E)`$ between where the two curves forming the boundary of the unconverged roots cross the real axis. For the lower boundary we use the smallest imaginary part of the eigenvalue near the real axis. The upper estimate is made by running a circle through the points near the real axis and choosing the one with the largest radius. In both cases the determination of how close to the real axis the points should be is made by eye, choosing a cutoff that appears to give a reasonable fit to the boundary near the real axis. Ideally these are the points where the contour of the cloud of unconverged roots crosses the real axis, corresponding to a lower and an upper estimate for $`\mu _c`$ using the given ensemble.

The dependence $`\mathrm{\Delta }(N_E)`$ on the size of the ensemble, $`N_E`$, is illustrated in Fig. 5 where we show results for $`N=12`$, $`N=24`$ and $`N=48`$. The thickness of the cloud of zeros, $`\mathrm{\Delta }(N_E)`$, shows a logarithmic dependence on $`N_E`$. We have included logarithmic fits for the range $`N_E=10^5`$ through $`10^8`$ ($`10^6`$ for $`N=12`$). Of course, once $`\mathrm{\Delta }(N_E)`$ becomes close to zero, it stays that way. The slopes of the fits for $`N=24,48`$ vary only slightly with $`N`$.

Finally, the issue of most practical interest is how the number of configurations $`N_E`$ required to achieve a fixed value of $`\mathrm{\Delta }(N_E)`$ varies with the size of the matrix. This was estimated by fitting lines to the available data for $`\mathrm{\Delta }(N_E)`$ versus $`\mathrm{log}_{10}(N_E)`$. As we saw in Fig. 5 this does give nice linear fits. We then chose values of $`\mathrm{\Delta }(N_E)`$ and calculated where the linear fits intersected these values for each value of $`N`$. Our results are plotted in Fig. 5. We conclude that the number of configuration required to obtain a given precision increases exponentially with the matrix size $`N`$. This is the most important conclusion of this paper. In the remaining sections we will attempt to reinforce this conjecture by studying the sensitivity of the zeros of the exact partition function to small random perturbations of the corresponding polynomial coefficients.

### B Perturbed exact zeros

It is a well known fact that extreme care has to be taken when calculating zeros of very high order polynomials. In particular, this is the case for a finite representation of a partition function, such as the one under investigation, where one is ultimately interested in large values of $`N`$. In our numerical experiments the polynomial coefficients are obtained as a result of a statistical averaging process. The averaged coefficients are approximations of the exact ones. As the simulation proceeds, the error decreases.

We wish to investigate the sensitivity of the zeros of the partition function to the precision with which the coefficients are obtained. We consider the exact polynomial and add a fixed relative error to each coefficient:

$`\stackrel{~}{c}_k=c_k(1+R_kϵ),`$ (35)   
where $`R_k`$ are random numbers between $`[1,1]`$, and $`ϵ`$ is a small real positive number. We then calculate the zeros of the polynomial $`\underset{k}{}\stackrel{~}{c}_kz^k`$.

The effect of obtaining approximate coefficients by Glasgow averaging is quite similar to that of simply adding “artificial noise” to the exact coefficients. In Fig. 6 we plot the roots obtained for $`N=32`$ via Glasgow from $`N_E=10^7`$ configurations and the roots of the exact polynomial perturbed by noise of magnitude $`ϵ=10^3`$. The patterns of the two sets of roots are hardly distinguishable.

In this subsection we study the effect of such small random perturbations. In particular, we are interested in the dependence of the precision $`\mathrm{\Delta }(N_E)`$ (in obtaining $`\mu _c`$) on $`N`$ and $`ϵ`$. This approach has the advantage that we can consider larger values of $`N`$ than in the case of direct Glasgow averaging. We expect that the relation between the error parameter $`ϵ`$ and the number of Glasgow configurations necessary to achieve the same accuracy is given by the central limit theorem, $`ϵ1/\sqrt{N_E}`$.

One striking feature for larger $`N`$ is the extreme sensitivity of the roots to small perturbations. For $`N=96`$, a noise factor of $`10^{18}`$ already leads to a 20% error in the critical value of the chemical potential. Therefore, computing the zeros for $`N=96`$ is already beyond the capability of standard double precision. In all our calculations we use multiprecision arithmetic, implemented either using the GNU multiprecision package or the one made publicly available by NASA .

In Fig. 8 we plot zeros with different degrees of artificial noise (see label of figure) for $`N=96`$. The zeros for $`ϵ=10^{25}`$ coincide with the exact roots. The noise factor for the remaining four sets of zeros increases with the same factor of $`10^3`$ from one set to the next. The intercepts of the edge of the eigenvalue cloud with the real axis are practically equally spaced for the different sets of zeros. The precision $`\mathrm{\Delta }(ϵ)`$ appears to be a linear function of $`\mathrm{log}(ϵ)`$ until $`\mathrm{\Delta }(ϵ)`$ vanishes.

In Fig. 8 we plot $`\mathrm{\Delta }(ϵ)`$ versus $`\mathrm{log}_{10}(ϵ)`$ for $`N=48,96,192,384`$. For error values not too for away from convergence, we observe a linear relationship between $`\mathrm{\Delta }(ϵ)`$ and $`\mathrm{log}(ϵ)`$. In all cases convergence is approximately achieved when $`\mathrm{log}_{10}(ϵ)N/4`$. This shows that the accuracy required to achieve a given precision increases exponentially with $`N`$. These results are consistent with results obtained with Glasgow averaging where the relative variance of the coefficients is roughly constant and varies as $`1/\sqrt{N}_E`$ with only a weak dependence on $`N`$.

We conclude that adding random noise to the exactly known coefficients has the same effect on the roots of the partition function as Glasgow averaging. This confirms once more that an exponentially large number of configurations is required to obtain a given precision. Extrapolating the $`N`$-dependence of the number of configurations needed given by the linear interpolation in Fig. 5 to lattice simulations leads to extremely large numbers for even a small lattice size. Consider for example a matrix size of $`N=128`$ (so $`D`$ in (5) is $`256\times 256`$), corresponding to a $`4^4`$ lattice with one degree of freedom per site. To achieve a precision of $`\mathrm{\Delta }=0.3`$ we would need approximately $`N_E=10^{28}`$ configurations, and for $`\mathrm{\Delta }=0.2`$ we would need $`N_E=10^{32}`$ configurations. A much more reasonable number of configurations can only be achieved if we consider a small lattice of $`2^4`$ which corresponds to $`N=24`$ for Kogut-Susskind (staggered) fermions. Here we would need $`N_E=10^6`$ configurations for $`\mathrm{\Delta }=0.3`$ or $`N_E=10^8`$ configurations for $`\mathrm{\Delta }=0.2`$. This exponential dependence was noted previously . In , a signal was achieved for a lattice of $`2^4`$ with $`10^5`$ configurations but for a $`4^4`$ lattice only a very weak signal was found after $`2\times 10^6`$ configurations.

### C Further analysis of the convergence of Glasgow averaging

In the previous two subsections we hope to have convinced the reader that the phenomena accompanying Glasgow averaging are nothing but the effect of knowing the polynomial coefficients of the partition function only with limited precision. This follows from the fact that the pattern of the false (unconverged) zeros is reproduced by adding small random numbers to the exact polynomial coefficients.

In our model the zeros close to the real axis are the last to converge. These are also the most interesting in practice, since they determine the critical value $`\mu _c`$. It is conceivable that for certain situations in QCD (such as with finite temperature) the roots close to the real axis are among the faster converging ones, which would provide a glimmer of hope for the Glasgow method.

In this subsection we consider the effect of artificial noise in detail. We wish to understand this phenomenon as well as why the false zeros are always concentrated in a cloud of a well defined shape. We will find that there is a clear correlation between the magnitude of the partition function and the stability of the zeros.

In Fig. 9 we give a topographic map of the absolute value of the RMM partition function on a logarithmic scale as a function of $`\mu `$. That is, we plot the level curves of the free energy per site given (up to a constant) by $`\mathrm{ln}(|Z_N|^2)/N`$. This is the natural quantity to study, since the partition function scales exponentially with $`N`$ . We observe a discontinuity in the derivative of the absolute value of the partition function in the complex $`\mu `$ plane. This is the locus of the zeros of the partition function. The deepest points are those around the (real) critical value of the chemical potential. A quick comparison with the preceding scatter plots should convince the reader that the higher the value of $`\mathrm{ln}(|Z|)`$ the more robust are the zeros in that neighborhood. However, the level curves do not match the contour of the cloud of false zeros. The situation a is little more subtle as we will see below.

#### 1 The noisy partition function

It is useful to separate the approximate (“noisy”) partition function into the exact one and the part due to the error in the coefficients,

$`Z_{\mathrm{tot}}(\mu )=Z_0(\mu )+Z_{\mathrm{err}}(\mu ).`$ (36)   
The partition function $`Z_{\mathrm{err}}(\mu )`$ is a polynomial whose coefficients are the differences between the exact coefficients and the approximate ones,

$`Z_{\mathrm{err}}(\mu )={\displaystyle \underset{k=0}{\overset{N}{}}}\delta c_k\mu ^k;\stackrel{~}{c}_k=c_k+\delta c_k.`$ (37)   
The differences are generated as $`\delta c_k=ϵR_kc_k`$ in the ‘artificial noise’ case. Our initial assumption was that for true Glasgow averaging the relative error in the coefficients is approximately the same for all zeros. This assumption was reinforced by the similarity between the two types of results discussed in previous subsections. In the following, we will discuss mainly the ‘artificial noise’ partition function.

In terms of $`Z_{\mathrm{err}}(\mu )`$, an explanation of the qualitative picture we have observed is the following. In the regions of the complex $`\mu `$ plane where $`|Z_{\mathrm{err}}(\mu )||Z_0(\mu )|`$ one may ignore the $`Z_{\mathrm{err}}(\mu )`$. The roots of the total partition function located in these regions coincide to a good degree with the corresponding exact roots. On the contrary, in regions where $`Z_{\mathrm{err}}(\mu )`$ dominates, the roots are determined by the latter and their general pattern has no similarity to their exact counterparts. This is consistent with the fact that the robust zeros are the ones situated in the region with higher $`\mathrm{ln}(|Z_0|)`$. The shrinking of the cloud of false zeros with decreasing $`ϵ`$ is a consequence of the corresponding decrease in the magnitude of $`Z_{\mathrm{err}}`$.

Throughout most of the complex $`\mu `$ plane, $`Z_0`$ is clearly larger than $`Z_{\mathrm{err}}`$. However, near the roots of $`Z_0`$ the error part has its chance to dominate. This is because the value of the exact partition function is the result of a major cancellation in the region where $`Z_0`$ is well approximated by $`\mathrm{exp}(\mu ^2N)`$. The sum of the (finite) series is much smaller than the general term, which is typically of order $`1`$. By multiplying each term of this series with a random number we spoil the cancellations. The sum of the perturbed series may therefore be significantly larger in absolute value than the exact sum.

#### 2 Sensitivity of the zeros

One arrives at similar conclusions by studying the infinitesimal variation of the roots. Let $`\mu _k`$ be the exact roots (so that $`Z_0(\mu _k)=0`$) and let $`\mu _k+\delta \mu _k`$ be the roots of the total noisy partition function. From our decomposition, we have

$`Z_0(\mu _k+\delta \mu _k)+Z_{\mathrm{err}}(\mu _k+\delta \mu _k)=0{\displaystyle \frac{|Z_{\mathrm{err}}(\mu _k+\delta \mu _k)|}{|Z_0(\mu _k+\delta \mu _k)|}}=1.`$ (38)   
Since the modified roots are not zeros of $`Z_0`$ or $`Z_{err}`$, the above equality is not fulfilled trivially. It indicates that the false zeros should be located in the transition region where $`|Z_0|`$ and $`|Z_{err}|`$ are comparable, i.e., on the border of the region where $`Z_{err}`$ dominates, rather than scattered inside it. In Fig. 10 we show how the contour of the cloud of false zeros is obtained using the formula above.

We can also make a quantitative estimate of how well converged a given root is. If we expand the first half of (38) we get

$`Z_0(\mu _k)+\delta \mu _kZ_0^{}(\mu _k)+Z_{\mathrm{err}}(\mu _k)+\delta \mu _kZ_{\mathrm{err}}^{}(\mu _k)0.`$ (39)   
The first term is simply zero. The last term we will neglect as being of higher order in $`\delta \mu _k`$. The remaining two terms can be rearranged to give

$`\delta \mu _k{\displaystyle \frac{Z_{\mathrm{err}}(\mu _k)}{Z_0^{}(\mu _k)}}.`$ (40)   
In other words, the variation of a given root $`\mu _k`$ is proportional to the ratio of $`Z_{err}`$ to $`Z_0^{}`$. If this quantity is negligible, the root is close to the exact one. If the ratio is close to $`1`$ or larger, the shift in $`\mu _k`$ is large, and the root is not obtained correctly. Therefore we can say that in the region where $`|Z_{\mathrm{err}}||Z_0^{}|`$, the roots of $`Z_{tot}`$ are reliable.

#### 3 Error piece as additional phase

We may check our arguments in the previous sections by comparing the magnitudes of the exact, the total and the noise partition functions (or free energies) in the complex $`\mu `$ plane. Fig. 12 shows the absolute value of the free energy per site, $`\mathrm{ln}(Z_0)/N`$ in one quadrant of the complex $`\mu `$ plane. The location of the exact zeros is on the cusp which runs from $`\mu =i`$ to $`\mu =0.527`$.

In Fig. 12 we have the same type of plot now for the total partition function with a noise factor $`ϵ=10^8`$. The cusp in the noisy result bifurcates. The exact and the noisy surfaces coincide exactly except for the region between the two new branches of the cusp, where the noisy partition function is larger. The new cusps coincide with the locus of the false zeros. A few false zeros are also scattered inside the region. It is clear from Fig. 13 where we plot the surfaces corresponding to $`Z_0`$ and $`Z_{\mathrm{err}}`$ simultaneously that the new cusps are located at the intersection of the free energy surfaces. The exact piece dominates everywhere except for inside the bifurcation, where the error piece dominates. The $`\mu `$ dependence is so steep for both of them that the smaller piece becomes negligible very fast as one moves away from the intersection line.

When we discussed earlier the properties of polynomials as partition functions, we associated the different analytic functions which are approximated by the polynomial in different regions of the complex parameter space with the phases of the partition function. We mentioned that the zeros are typically located along the phase transition lines, since within any given phase the partition function is a smooth analytic function which does not vanish. The partition function $`Z_{\mathrm{err}}`$ has a quasi-analytic behavior similar to that of the partition function itself. The bifurcation of the line of zeros can be interpreted as the presence of an additional spurious ‘phase’ in the partition function, namely, the region where the error piece dominates. Then the fact that the false zeros are located on the $`|Z_{\mathrm{err}}|=|Z_0|`$ line is natural.

#### 4 Analytic approximation of the error partition function

The question is whether the error partition function can be at all approximated by an analytic function. In the present case the answer is very simple. Let us consider the error partition function,

$`Z_N^{\mathrm{error}}(\mu )=ϵ{\displaystyle \underset{k=0}{\overset{N}{}}}R_k{\displaystyle \frac{(N\mu ^2)^k}{k!}}.`$ (41)   
It is the same truncated inverse exponential we have discussed, only now each term is multiplied by a random number $`R_k`$ of order $`1`$. The terms in the original series conspire to achieve a major cancellation, of order $`e^{N\mathrm{Re}\mu ^2}`$. The random numbers spoil this, and each term in the series is left to fend for itself, and the sum is dominated by the largest terms. From the ratio of two consecutive terms in the sum,

$`{\displaystyle \frac{N\mu ^2}{k}}`$ (42)   
we conclude that the largest term (in absolute value) is the one with $`k=k_{\mathrm{max}}=|\mu ^2N|`$. If $`|\mu ^2|>1`$ then the largest term is the one with $`k=N`$. Of course the surrounding terms must have a significant contribution, but the end result should still be proportional to this largest term. For $`|\mu |<1`$ our estimate is therefore

$`|Z_N^{\mathrm{noise}}(\mu )|𝒞{\displaystyle \frac{k_{\mathrm{max}}^{k_{\mathrm{max}}}}{k_{\mathrm{max}}!}}{\displaystyle \frac{e^{N|\mu |^2}}{|\mu |\sqrt{N}}}`$ (43)   
In Fig. 14 we plot the absolute value of the two surfaces corresponding to the continuum limit, $`e^{\mu ^2N}`$ and $`e^N\mu ^{2N}`$, and the one corresponding to our estimate of the error partition function given above. The intersections of the three surfaces follow the pattern of the corresponding zeros. We also checked that (43) approximates the error partition function well.

As an added bonus, we can now explain the scaling with $`N`$ of the required precision. We found that the error partition function, i.e., the exact partition function whose coefficients have been multiplied by random numbers of order $`1`$, is well approximated by an analytic expression. This expression shares an important property of the true partition function, namely, it scales exponentially with $`N`$.

$`{\displaystyle \frac{1}{N_1}}\mathrm{ln}(Z_{N_1}^{exact}(\mu )){\displaystyle \frac{1}{N_2}}\mathrm{ln}(Z_{N_2}^{exact}(\mu ));{\displaystyle \frac{1}{N_1}}\mathrm{ln}(Z_{N_1}^{err}(\mu )){\displaystyle \frac{1}{N_2}}\mathrm{ln}(Z_{N_2}^{err}(\mu ))`$ (44)   
The magnitude of the error partition function is also controlled by the factor $`ϵ`$ which mimics the precision to which the coefficients are approximated in the averaging process. The locus of the false zeros is controlled by the relative magnitude of the true partition function and the error partition function. To have the same pattern of zeros, we must also scale $`ϵ`$ exponentially with $`N`$, $`ϵ(N)=\alpha ^N`$ .

## IV Conclusions

We have investigated Glasgow averaging using a random matrix model at nonzero chemical potential. We have found that in a quenched ensemble the method converges, but that it requires an exponentially large number of configurations.

The roots of the averaged polynomial are initially distributed similarly to the eigenvalues of individual configurations. As the averaging proceeds, the roots approach their exact values. After averaging over a finite number of configurations, the roots clearly separate into two groups. Some roots are close to the corresponding exact ones. Typically, these are the zeros far from the real axis. The remaining roots are situated in a cloud around the intersection of the real axis and the locus of the exact zeros. The zeros outside the cloud are practically exact, while those inside and on the boundary of the cloud are very badly determined. They cannot be traced to individual exact zeros. The shape of the cloud is similar for different matrix sizes. It shrinks as more configurations are taken into account. However it only shrinks as a logarithmic function of the total number of configurations.

By interpolating the number of configurations needed to reach a given precision for several matrix sizes, we were able to estimate the dependence of the number of configurations needed for a given matrix size. Our conclusion is that the number of configurations needed to reach a given level of precision grows exponentially with the size of the matrix.

The results of the corresponding unquenched simulations, using an ensemble generated with $`N_f=1`$ and $`\mu =0`$ are similar to the quenched ones. The unquenched ensemble may have been helpful in improving the overlap between the simulation and the ‘true’ ensemble corresponding to fixed nonzero $`\mu `$. However, this would still necessitate simulating an ensemble at nonzero $`\mu `$, which is precisely what we were trying to avoid. The exponential statistics observed by us are more likely to be the signature of the sign problem itself, i.e., the magnitude of the cancellation brought about by the varying phase of the determinant. Hence the Glasgow method is unable to surmount the sign problem. However, in a problem where the latter is absent, such as $`SU(2)`$ simulations , the Glasgow method – even quenched – should be helpful. It would be interesting to see how the overlap issue manifests itself in this case.

We obtained results very similar to those of Glasgow averaging by perturbing the exact coefficients in the RMM partition function polynomial. We studied empirically the dependence of the reliability of the zeros on the precision with which the coefficients are known and on the size $`N`$ of the model matrix for larger matrices. Just like in the case of averaging, this dependence is exponential. That is, given a desired error bound for the zeros, the necessary precision on the coefficients grows exponentially with $`N`$. The extrapolation of this dependence to Glasgow averaging translates into exponentially large statistics, since the precision on the coefficients should be proportional to the square root of the number of configurations.

There is a correlation between the phenomenon of slower or faster converging zeros and the magnitude of the continuum partition function. The zeros that converge slowly are in the region where the partition function is suppressed. Large cancellations require better precision, hence more statistics. More formally, the effect of perturbing the coefficients can be understood as the addition of an extra (error) polynomial to the true partition function. This extra piece is found to scale exponentially with $`N`$, just like the true partition function. The effect on the phase structure can be seen as the introduction of a spurious phase, which replaces the true ones in the regions of parameter space where the error partition function dominates. The zeros in these regions follow the modified phase boundaries. The scaling property of the error partition function explains the need for exponential statistics in Glasgow averaging.

The negative result regarding unquenched simulations at $`\mu =0`$ is perhaps disappointing. It indicates that the quenched ensemble and the unquenched ensemble at $`\mu =0`$ are equally relevant to the behavior of the model close to the critical $`\mu `$. The upside of the failure of the unquenched ensemble at $`\mu =0`$ is that in future applications of the Glasgow method it might be worth trying to use a quenched ensemble. We remind the reader that using the quenched ensemble is not equivalent to the quenched approximation.

Another glimmer of hope is in the phenomenon of fast converging zeros. In the present random matrix model the zeros that determine the critical $`\mu `$ are the last to converge. From our analysis there is no indication that the sensitive zeros are generally those close to the real (physical) axis. There is a possibility that in QCD or in other interesting non-Hermitian models the critical parameter values are determined by the robust zeros. It is hard to make any statement in this respect from the currently available QCD data . Perhaps a schematic model with features closer to QCD in terms of eigenvalue distribution in the complex $`\mu `$ plane would help clarify this issue. In general, it would be interesting to see how much of the analysis in this work regarding zeros of approximately known polynomial partition functions applies to other models.

## V Acknowledgements

We thank I.M. Barbour, S. Hands, E. Klepfish, M.P. Lombardo and S.E. Morrison for fruitful discussions. R.D. Amado and M. Plümacher are thanked for several critical readings of the manuscript.

This work was supported in part by NSF grants PHY-98-00098 and PHY-97-22101 as well as the DOE grant DE-FG-88ER40388. Part of the calculations presented in this paper have been carried out at the National Scalable Cluster Project Center at the University of Pennsylvania, which is supported by a grant from the NSF. The multiprecision calculations were performed using the GNU multiprecision package or the package made available by NASA .