Random matrix analysis of the QCD sign problem for general topology

Random matrix analysis of the QCD sign problem for general topology

Jacques Bloch and Tilo Wettig
Institute for Theoretical Physics, University of Regensburg, 93040 Regensburg, Germany
Email: ,
jacques.bloch@physik.uni-regensburg.de tilo.wettig@physik.uni-regensburg.de

Motivated by the important role played by the phase of the fermion determinant in the investigation of the sign problem in lattice QCD at nonzero baryon density, we derive an analytical formula for the average phase factor of the fermion determinant for general topology in the microscopic limit of chiral random matrix theory at nonzero chemical potential, for both the quenched and the unquenched case. The formula is a nontrivial extension of the expression for zero topology derived earlier by Splittorff and Verbaarschot. Our analytical predictions are verified by detailed numerical random matrix simulations of the quenched theory.

Random matrix theory, Lattice QCD, Quark chemical potential
preprint: December 1, 2008; revised January 15, 2009

1 Introduction

The theory of the strong interactions, also called Quantum Chromodynamics (QCD), describes the interactions between quarks and gluons and is responsible for the existence of hadrons. Lattice-regularized QCD allows for the description of low-energy properties and other nonperturbative phenomena in QCD and has the salient property that it can be systematically improved towards the continuum limit. In lattice QCD, space-time is discretized and the functional integral of the quantum field theory is performed by a Markov-chain Monte-Carlo method.

An important subject of study is the behavior of QCD in an environment exhibiting an abundance of particles over anti-particles. Such conditions arise, e.g., in astrophysical objects like neutron stars, and can be reproduced in heavy-ion collision experiments. Part of the interest comes from the existence of various phases in QCD, which are usually exemplified by means of the QCD phase diagram [Stephanov:2007fk]. To study QCD at nonzero baryon density a quark chemical potential is introduced in the QCD Lagrangian. (In the following, we will omit the qualifier “quark” and only speak of a chemical potential.) In the presence of a chemical potential the QCD Dirac operator is no longer anti-Hermitian, i.e., its eigenvalues spread into the complex plane and its determinant will generically be complex.

In lattice QCD the effect of dynamical fermions can be integrated out, leaving behind the determinant of the Dirac operator. Dynamical lattice simulations for QCD at nonzero chemical potential are problematic because the fermion determinant is complex and hence its real part can be negative, which prohibits its incorporation in the weight of the Monte-Carlo sampling. This is the so-called sign problem, which also occurs in other theories and has been the subject of a large number of investigations in recent years (for an incomplete list see, e.g., Refs. [Ambjorn:2002pz, Ambjorn:2004jk, Troyer:2004, Osborn:2005ss, Imachi:2006mw, Fukushima:2006uv, Ejiri:2007ga, Aarts:2008wh]). While many of these works are concerned with a solution of the sign problem by various clever ideas, here we concentrate on an analytical study of the sign problem, in the hope that the results we derive will contribute to its solution. The severeness of the sign problem depends on the magnitude of the chemical potential, and it is therefore illuminating to investigate the relation between the phase factor of the determinant and the chemical potential.

Chiral random matrix theory (chRMT) is a useful auxiliary in the study of the spectral properties of the Dirac operator in QCD [Shuryak:1992pi, Verbaarschot:2000dy, Verbaarschot:2005rj]. Indeed, to leading order in the -regime of QCD the spectral properties of the Dirac operator are universal and can be described by chRMT [Basile:2007ki]. In the presence of a chemical potential this correspondence is still valid even though the Dirac operator is now non-Hermitian. Appropriate random matrix models have been developed [Stephanov:1996ki, Akemann:2002ym, Akemann:2002js, Osborn:2004rf], and their correspondence with QCD at nonzero chemical potential was verified successfully, see Ref. [Akemann:2007rf] for a review. The agreement of the microscopic spectral properties of the Dirac operator with the predictions of chiral random matrix theory has been confirmed for quenched lattice QCD simulations with chemical potential using the staggered operator [Akemann:2003wg], and more recently using the overlap operator [Bloch:2006cd, Akemann:2007yj]. The latter operator has the interesting property that it satisfies the Ginsparg-Wilson relation and the trace anomaly at finite lattice spacing and can therefore have exact zero modes [Ginsparg:1981bj, Narayanan:1993sk, Narayanan:1994gw, Neuberger:1997fp, Luscher:1998pqa, Bloch:2007xi]. This allowed us to verify the predictions of chRMT at nonzero chemical potential for both zero and nonzero topology. The comparison between lattice QCD and chRMT also allows for a determination of the low-energy constants and of chiral perturbation theory.

Motivated by this agreement, we expect that a study of the behavior of the fermion determinant in chRMT will give us, in certain well-defined limits, important information about the sign problem that is encountered in dynamical QCD simulations at nonzero chemical potential. In Ref. [Splittorff:2007ck] Splittorff and Verbaarschot derived a solution for the average phase factor of the determinant in the microscopic limit of QCD (see Sec. 3.2 for a description of this limit) for the case of trivial topology using chRMT at nonzero chemical potential. However, to compare the overlap data of Ref. [Bloch:2006cd] with chRMT one also needs the RMT predictions for the average phase factor for general topology. The derivation of a formula for general topology is the main goal of this paper. As will be seen, the final expression contains two distinct parts. The first part is the generalization of the integrals representing the solution in Ref. [Splittorff:2007ck] from zero to arbitrary topology. The second part is a low-degree bivariate polynomial in mass and chemical potential which is absent for topological charge . For it gives an important contribution to the average phase factor of the fermion determinant, especially for small mass. As the mass goes to zero only this term remains and completely determines the value of the average phase.

An important ingredient of the derivation is the ability to write the phase factor of the determinant as a ratio of characteristic polynomials. This quantity is recurrent in random matrix studies, both for theories with real [Fyodorov:2002jw, Strahov:2002zu] and with complex eigenvalues [Akemann:2004zu], and its average can be computed in terms of Cauchy transforms of the orthogonal polynomials of the theory. To determine the phase factor of the determinant, the relevant Cauchy transform was computed in Ref. [Splittorff:2007ck] and expressed in terms of one-dimensional integrals for zero topology, i.e., for square random matrices. In the present paper we extend the solution of the Cauchy transform to the case of rectangular matrices. This solution could also be relevant for other applications, like those involving time series, where one matrix dimension is typically much larger than the other.

The structure of this paper is as follows. In Sec. 2 we describe the chiral random matrix model at nonzero chemical potential. In Sec. 3 we show how the microscopic limit of the phase of the fermion determinant, in both the quenched and the unquenched case, can be formally computed for such a matrix model in terms of a complex Cauchy transform integral. This two-dimensional integral is strongly oscillating, and in Sec. 4 we apply and extend the method of Ref. [Splittorff:2007ck] to transform this integral into a much simpler and better behaved expression, involving only one-dimensional integrals and a short double sum (or bivariate polynomial). Explicit results for the quenched and unquenched cases as well as for the chiral and thermodynamic limits are given in Sec. 5. In that section we also verify the analytical predictions for the quenched case by random matrix simulations for various values of the topological charge. We conclude in Sec. 6. A number of technical details are worked out in several appendices.

2 Non-Hermitian chiral random matrix model

To leading order in the -regime of QCD the spectral properties of the Dirac operator can be described by chRMT. In the presence of a chemical potential the Dirac operator is no longer anti-Hermitian, and in the non-Hermitian chiral random matrix model introduced by Osborn [Osborn:2004rf] it takes the form


where the matrices and are complex random matrices of dimension , distributed according to the Gaussian weight function


For a detailed analysis of this model, see also Ref. [Akemann:2004dr]. For the conversion of random matrix units to physical units, see the beginning of Sec. 3.2.

The parameter will be taken to infinity when computing the microscopic limit (see Sec. 3.2). The matrix in Eq. (2.1) has exact zero modes, which allows us to identify with the topological charge. In the following, we keep fixed as and assume without loss of generality that . (For we can simply replace by in the analytical results that will be computed below in the large- limit.) The nonzero eigenvalues of come in pairs . For , the are purely imaginary.

For fixed , the partition function of the random matrix model is given by


where the integration measure is defined by


is the number of dynamical quarks, and the are the quark masses. The quenched case corresponds to , i.e., the fermion determinants are absent.

To perform the integration over and , it is convenient to go to an eigenvalue representation of the random matrix . As shown in Ref. [Osborn:2004rf], the partition function can be rewritten, up to a normalization constant that depends on and , as an integral over the ,


where we introduced , the integrals over the are over the entire complex plane,


is a Vandermonde determinant, the weight function is given by


and is a modified Bessel function. The quenched partition function will be denoted by .

The ensemble average of an observable is given by


When there is no danger of confusion we will omit one or both of the subscripts on .

Our derivation will follow the general line of arguments given in Ref. [Splittorff:2007ck] for , with the necessary generalizations to arbitrary topology . To analyze the spectral properties of the random matrix model it is useful to introduce the orthogonal polynomials corresponding to the weight function (2.7) [Osborn:2004rf],


where is the generalized or associated Laguerre polynomial of order and degree . These orthogonal polynomials satisfy the orthogonality relation


with norm


For later use we also introduce the Cauchy transform of the orthogonal polynomials defined by


3 Phase factor of the fermion determinant

3.1 The phase factor as a complex Cauchy transform

The Dirac operator describing a massive fermion is defined as , where we assume that is real. If we write its determinant as , the phase factor can be extracted by forming the ratio


In this expression, is the mass of a valence quark. From the physics point of view, an interesting quantity is the ensemble average of with two light dynamical quarks that have the same mass as the valence quark. This quantity tells us how the two-flavor determinant in the weight function oscillates. For simplicity, we shall refer to as the phase factor of the determinant, even though it is really the phase of the square of the determinant.

Because of the symmetries of (2.2), each matrix appears in the ensemble average with the same probability as its Hermitian conjugate. As the corresponding determinants are complex conjugate, the ensemble average of the phase factor is real. For strongly oscillating determinants the average phase factor will be close to zero, and the sign problem will be severe. On the other hand, for values of the chemical potential for which the average phase factor is close to unity one should still be able to perform dynamical simulations.

For each topic that is treated here and in the following sections, we will first address the quenched case and then generalize to the unquenched case. The virtue of this approach is that the quenched case already contains the essential ingredients, but the arguments and the notation can be kept simple. The generalization to the unquenched case is straightforward but leads to somewhat more complicated expressions.

3.1.1 Quenched case

The quenched ensemble average for the phase factor is given by


The quantity is the partition function of a random matrix model with one fermionic quark and one conjugate bosonic quark, see Ref. [Splittorff:2007ck] for a detailed discussion.

Using the formalism developed in Refs. [Bergere:2004cp, Akemann:2004zu], the quenched average of ratios of characteristic polynomials can be written in terms of the orthogonal polynomials (2.9) and their Cauchy transforms (2.12). Applying this formalism to the quenched average phase factor (3.2) gives the compact expression


which is a complex integral over the orthogonal polynomials due to the Cauchy transforms. This expression (and its analog for the unquenched case, see Sec. 3.1.2) will prove to be very useful to compute the phase factor. Inserting the Cauchy transform (2.12) in Eq. (3.3) yields an integral over orthogonal polynomials,


where we also used . The well-known recurrence relation for the generalized Laguerre polynomials,


translates into a recurrence relation for the defined in Eq. (2.9),


Since the determinant remains unchanged when forming linear combinations of its columns, we can rewrite the phase factor (3.4) using the recurrence (3.6) as


where all the orthogonal polynomials are now of equal degree.

3.1.2 Unquenched case

In the presence of dynamical fermion flavors with masses , the phase factor for a valence quark of mass is given by


where is given by Eq. (2.5). This can be written as a ratio of two partition functions,


where, in analogy to Eq. (3.2), is the chRMT partition function with fermionic quarks and one conjugate bosonic quark. Both partition functions can be computed using the results of Ref. [Akemann:2004zu], but to apply these results we need to change the normalization and divide both and by the quenched partition function. Then the partition functions can be interpreted as averages of ratios of characteristic polynomials in the quenched ensemble. Applying the formalism of Ref. [Akemann:2004zu] to the numerator of Eq. (3.9) (with modified normalization) we find


where the matrix in the determinant is of size . Substituting the definition of the Cauchy transform (2.12) gives


With the recurrence relation (3.6) for the orthogonal polynomials this becomes


The denominator in Eq. (3.9) (with modified normalization) can be written as [Akemann:2004zu]


and using the recurrence relation (3.6) for the orthogonal polynomials this can be rewritten as


3.2 Microscopic limit

Universal results, i.e., results that also apply to QCD, can be obtained from chRMT in the so-called microscopic regime. This regime is obtained by taking while keeping the rescaled parameters , , and fixed, and rescaling the spectrum using . The rescaled random matrix parameters can be converted to the physical parameters , , and using the relations , , and , where is the four-volume.111To be more precise, one should distinguish random matrix parameters and physical parameters in the relations , , and . This distinction makes it explicit that the limits and can be decoupled. Furthermore, the pion mass can be introduced through the combination , where we have used the Gell-Mann–Oakes–Renner relation (assuming equal quark masses).

We now introduce the microscopic limits (denoted by a subscript ) of the orthogonal polynomials, the norm, and the weight function, respectively. They are worked out in App. A, and we obtain


3.2.1 Quenched case

In terms of the above definitions, the microscopic limit of the quenched average phase factor (3.7) is given by


As expected, the dependence on has dropped out, leaving a finite microscopic limit for the average phase factor. Substituting the asymptotic results from Eqs. (3.15)–(3.17) yields


where we renamed the integration variable back to . This equation can be rewritten compactly as


where we defined the integral


which is closely related to the microscopic limit of the Cauchy transform (2.12). For the quenched case, this integral is only needed for , but as we shall see in the next subsection, in the unquenched case it will be needed for .

3.2.2 Unquenched case

We now take the microscopic limit of Eqs. (3.12) and (3.14). In Eq. (3.12), the Vandermonde determinant is a product of factors for which the microscopic limit yields a dependence on , while the explicit mass and factors in the determinant yield a factor . After also introducing the microscopic limits (3.15), (3.16), and (3.17) for the orthogonal polynomials, their normalization factor, and the weight function, respectively, one finds


where we have again renamed the integration variable from back to and introduced the notation


In the microscopic limit of Eq. (3.14), the Vandermonde determinant yields a factor which exactly cancels the factor coming from the explicit mass factors in the determinant. After introducing the microscopic limit (3.15) of the orthogonal polynomials we find




The microscopic limit of the average phase factor (3.9) is given by the ratio of Eqs. (3.22) and (3.24). The dependence on drops out to give


which can be rewritten using the integral definition (3.21) as


We have thus reduced the problem of calculating the phase factor to the calculation of the two-dimensional integral in Eq. (3.21) for . This integral will be computed in Sec. 4.

3.2.3 Equal mass fermions

We now consider Eq. (3.27) for the special case in which all dynamical fermions have the same mass as the valence quark. To simplify the general expression we perform a Taylor expansion of the entries of the determinant around ,