A new method to study lattice QCD at finite temperature and chemical potential

Due to the sign problem, it is exponentially difficult to study QCD on the lattice at finite chemical potential. We propose a method --an overlap improving multi-parameter reweighting technique-- to alleviate this problem. We apply this method and give the phase diagram of four-flavor QCD obtained on lattices 4^4 and 4\cdot6^3. Our results are based on {\cal{O}}(10^3-10^4) configurations.

: The average of the quark condensates at β=5.085 as a function of Im(µ), for direct simulations (squares; their sizes give the errors), our technique (crosses) and Glasgow-type reweighting (dots).
configurations for these partition functions an additional Ferrenberg-Swendsen [17] reweighting can be done in the other parameter, thus in β, too. For Re(µ) =0 and/or Im(β) =0 significant cancellations of the complex phases appear, but exactly this feature is used in the determination of the zeros of the partition functions (Z), when looking for the transition point. Simultaneously reweighting in the two parameters µ and β we can keep the system on the transition line, which can be controlled e.g. by the inspection of the Lee-Yang zeros [18] of Z at complex β. (Note, that the idea of performing a reweighting near the QCD critical line was already suggested in [8]). This technique gives a good overlap with the original transition-like states. In principle any other parameter can be used for this type of reweighting.
We test this method and illustrate its success compared to the Glasgow method. We present the exploratory results for the µ-T phase diagram of the dynamical n f =4 staggered QCD. Simulations were done on L t = 4 lattices of 4 4 and 4 · 6 3 . We estimate the phase diagram in physical units using the ρ mass (m ρ ) as the definition of the scale.
Due to the small lattices and large spacings our estimate has systematic uncertainties, which can be reduced by approaching the continuum limit. The study of this limit is clearly not the goal of the present Letter.
The Letter is organised as follows. Section II. presents the overlap ensuring multi-parameter reweighting technique. Section III. illustrates the applicability of the technique using n f =4 dynamical QCD and gives the phase diagram in physical units. We conclude in Section IV.

Overlap improving multi-parameter reweighting
Let us study a generic system of fermions ψ and bosons φ, where the fermion Lagrange density isψM (φ)ψ. Integrating over the Grassmann fields one gets: where α denotes a set of parameters of the Lagrangian. Eg. in the case of QCD with staggered quarks α consists of β, the quark mass (m q ) and µ (which is included as exp(µ) and exp(−µ) multiplicative factors of the forward and backward timelike links, respectively). For some choice of the parameters α=α 0 importance sampling can be carried out (e.g. for Re(µ)=0). Rewriting eq. (1) one obtains Now we treat the terms in the curly bracket as an observable -which is measured on each independent configurationand the rest as the measure. It is well known that changing only one parameter of the set α 0 the ensemble generated at α 0 provides an accurate value for some observables only for very high statistics. This is ensured by important but rare fluctuations as the mismatched measure occasionally sampled the regions where the integrand is large. This is the so-called overlap problem. Note however, that we have several parameters and the set α can be adjusted to α 0 to ensure much better overlap than obtained by varying only one, single parameter. Since the calculation of the determinants is expensive the most straightforward way to obtain a good overlap is to fix the parameters of the matrix M and adjust the parameters which appear only in the bosonic action (in other words perform a Ferrenberg-Swendsen reweighting based on the bosonic part of the curly bracket). By simulating at a given set of bosonic couplings and using Ferrenberg-Swendsen reweighting we can get information on the system at other values of the couplings, even for complex ones. At finite volumes Z(α) has zeros -thus the free energy singularities-for complex values of these couplings. Standard finite size scaling techniques can be used to analyse the volume dependence of the Lee-Yang zeros. When looking for these zeros we use reweighting for some bosonic couplings. Simultaneously changing two -or more-parameters we can ensure that the system is reweighted along a transition line. This can be monitored by inspecting the Lee-Yang zeros of Z.
3 Illustration: n f =4 dynamical QCD results For the case of lattice QCD at nonvanishing µ the above idea can be applied as follows. We will use two parameter reweighting, namely reweighting in β and µ. One performs the simulations at some β, m and µ with Re(µ)=0 (note that purely imaginary µ can be directly simulated). For each independent gauge configuration one calculates the average value of the plaquettes and the ratio of the determinants det M (µ ′ )/ det M (µ; Re(µ) = 0). For each µ ′ some δβ can be used to reweight with the measured plaquette values. By this way a much better overlap can be ensured than by reweighting only in µ (Glasgow method).
We have tested these ideas in four-flavor QCD with m q =0.05 dynamical staggered quarks. The molecular dynamics Monte-Carlo code of the MILC collaboration was used [19]. We checked that the determinants were calculated on independent configurations. Our statistical errors were obtained by a jackknife analysis.
In order to directly check the applicability of our method we first collected 1200 independent V=4·6 3 configurations at Re(µ)=Im(µ)=0 and used the Glasgow-reweighting and also our technique to study Re(µ)=0, Im(µ) =0. For our method we used the transition β (5.040) to generate the configurations, while for the Glasgow method β=5.085 was used. At Re(µ)=0, Im(µ) =0 and at β=5.085 direct simulations are possible. After performing these direct simulations as well, a clear comparison can be done. Figure 1 shows the predictions of the three methods for the quark condensate. The prediction of our method is in complete agreement with the direct results, whereas the predicition of the Glasgow-method is by several standard deviations off. Figure 1 indicates that the two-parameter reweighting is far more trustworthy than the single parameter one. Note, that imaginary chemical potential is a useful check on the proposed method, though it is different from a real chemical potential in biasing the ensemble towards non-zero baryon density.
Based on these experiences we expect that our method is superior also at Re(µ) =0. Next we study the physical Re(µ) =0 case. In order to have a better overlap and to check further the applicability of our reweighting technique we carried out simulations on 4 4 lattices at four different imaginary µ values Im(µ)=0,0.1,0.2,0.3. The results obtained by the different runs are in complete agreement after reweighting. In the following we use our largest sample generated at Im(µ)=0. On this smaller V we used 13000 independent configurations. On the 4 · 6 3 lattice Im(µ)=0 was used with 1200 independent configurations. The runs were carried out at the transition β values at Re(µ)=0.
Let us illustrate that we are really at a transition point for nonvanishing µ=0.3. Figure 2 shows the reweighted Polyakov-line and chiral condensate as a function of β. Furthermore, we check the correctness of our Lee-Yang reweighting approach by showing the structure of a histogram. As it can be seen the turning point (indicating the coexistence of the two phases, thus the transition) is at ≈ β=4.94. At µ=0.3 the histogram method predicts a equal-weight double-peak stucture at β c =4.939 (5), which is in complete agreement wit the prediction of the Lee-Yang zero technique (see later).
Re(µ) Re(β 0 ) for V = 4 4 Re(β 0 ) for V = 4 · 6 3 0. 4 Table 1. gives the real parts of the Lee-Yang zeros. Only zeros with the smallest imaginary parts are listed. The real parts of these zeros are usually used as a definition of the transition β at finite V. (Note, that the V→ ∞ limit of the imaginary parts tells the order of the transition, cf. [18,20].) Based on V=4 4 and 4 · 6 3 we estimated the V→ ∞ limit by 1/V scaling. The critical β in this limit is used to transform the results to physical units.
It is of particular interest to determine the phase diagram of QCD on the T -µ plane. Though the lattices we used are absurdly small and the spacing is quite large, it is illustrative to transform β, m q and µ into physical units. Several parameters can be chosen to fix the scale, they give quite different values at our β couplings. In the present analysis we fixed the scale by m ρ =770 MeV. For small β values, studied by the present paper, m ρ can be obtained by interpolating between the strong coupling regime [21] and the early measurements [22]. Figure 3 shows the phase diagram in physical units. The errorbars indicate the statistical uncertainties reached on our sample of only O(10 3 − 10 4 ) configurations. Note, that L t =4 lattices with the above definition of the scale restrict T to be larger than approximately 100 MeV (for clarity we used this value at the origin).

Conclusions, outlook
We proposed a method -an overlap improving multi-parameter reweighting technique-to numerically study non-zero µ and determine the phase diagram in the T -µ plane. We applied this technique for n f =4 QCD with dynamical staggered quarks. We showed that for Im(µ) =0 the predictions of our method are in complete agreement with the direct simulations, whereas the Glasgow method suffers from the well-known overlap problem. Based on rather small statistics we were able to determine the critical gauge couplings as a function of the real chemical potential, which result was transformed into physical units.
In this exploratory study we concentrated on the transition line separating the two phases. Clearly, the same technique can be applied to T -µ values somewhat below or above the line, for which the configurations should be collected at an appropriately chosen β -below or above the transition one-and at Re(µ)=0.
Note, that the factor in curly braces in (2) is of order exp[−V δ], with δ, of course, dependent on the parameters α and α 0 and the configuration φ. Making δ small by a judicious choice of α 0 seems to work reasonably well at finite temperature. Nevertheless, there are two apparent limitations of the method. For large enough volumes the reweighting factor will always be exponentially suppressed, and at zero temperature δ is larger than at finite T, thus the method most probably can not be applied to locate the transition at T=0. Our method can be easily applied to any number of Wilson fermions, whereas for n f =2 or n f =2+1 staggered fermions the situation is more complicated due to the ambiguity of the roots of the determinants [20].
Note, that the present reweighting does not provide a general solution to the sign problem in QCD, but could be useful to locate the critical point of QCD, for which we are interested in a relatively small µ and high T region of the phase diagram and may get away with relatively small volumes. Another application is to determine the curvature of the phase diagram [23].