New approach to lattice QCD at finite density; results for the critical end point on coarse lattices

All approaches currently used to study finite baryon density lattice QCD suffer from uncontrolled systematic uncertainties in addition to the well-known sign problem. We formulate and test an algorithm, sign reweighting, that works directly at finite μ = μB /3 and is yet free from any such uncontrolled systematics. With this algorithm the only problem is the sign problem itself. This approach involves the generation of configurations with the positive fermionic weight |Re det D(μ)| where D(μ) is the Dirac matrix and the signs sign(Re det D(μ)) = ±1 are handled by a discrete reweighting. Hence there are only two sectors, +1 and −1 and as long as the average 〈±1〉 ≠ 0 (with respect to the positive weight) this discrete reweighting by the signs carries no overlap problem and the results are reliable. The approach is tested on Nt = 4 lattices with 2 + 1 flavors and physical quark masses using the unimproved staggered discretization. By measuring the Fisher (sometimes also called Lee-Yang) zeros in the bare coupling on spatial lattices L/a = 8, 10, 12 we conclude that the cross-over present at μ = 0 becomes stronger at μ > 0 and is consistent with a true phase transition at around μB /T ∼ 2.4.

More precisely, the Taylor expansion method for non-infinitesimal µ requires the computation of high order µ-derivatives at µ = 0. It has the advantage that it provides well-defined physical quantities, namely the cumulants of the baryon number distribution at µ = 0 directly. However, the measurement itself leads to ever growing cancellations among fermion contractions as the order of the derivative increases. Furthermore even if a potentially large number of Taylor coefficients are computed with acceptable statistical uncertainty, the best case scenario is a reliable estimate of the radius of convergence. The Taylor expansion method will only provide information within this radius and extrapolation beyond it necessarily will involve uncontrolled systematics.
The second extrapolation method mentioned above involves simulating at imaginary µ where there is no sign problem, but the extrapolation from negative µ 2 to positive finite µ 2 requires assumptions about the functional form of the µ 2 dependence. This leads again to uncontrolled systematic uncertainty in the extrapolation, similar to the case of the Taylor method.

JHEP05(2020)088
The third popular method, the complex Langevin approach, is appealing because it is set up at finite µ directly but the precise set of necessary and sufficient conditions for it to give the correct result in QCD is so far unknown. A set of sufficient conditions for the correctness of the algorithm in general (some a priori, such as the holomorphicity of the action, and some a posteriori, such as the quick decay of the field distributions at infinity) has been proven [51,52], however these conditions are not satisfied in lattice QCD. Although one may formulate various tests of incorrectness and the lack of observed such signals may boost confidence in the correctness of the results, the systematic uncertainties associated with the potential breakdown of the algorithm cannot be estimated quantitatively. Numerical investigations indicate that present incarnations of the method break down at low temperatures. Whether an extension of the method capable of simulating also at low temperatures exists is a matter of ongoing research.
Finally the fourth method, reweighting from µ = 0, leads to the well-known overlap problem at some finite µ. This means that if a suitable weight is found, w(µ), which may depend on any number of further parameters [40,41] beyond µ, and expectation values are computed via O µ = Ow(µ) 0 / w(µ) 0 , then the histogram of w(µ) becomes wider and wider for increasing µ. Sampling the tail of the histogram becomes eventually prohibitively expensive and a reliable error estimate at finite statistics impossible. Furthermore, there is no sharply defined condition which would signal the presence of the overlap problem or absence thereof. In practice one may attempt to confirm the lack of the overlap problem from various statistical observations and may very well obtain reliable results, but the inherent systematic uncertainty will nevertheless linger.
Our motivation for the present paper is to devise an algorithm which is free of uncontrolled systematic uncertainties and has a well-defined set of conditions for its applicability. In other words we would like to have a trustworthy algorithm in the sense that results obtained with it are reliable with well-defined statistical uncertainties and have quantifiable, controlled systematic uncertainties. We will not solve the sign problem and do not aim to. Our approach involves the generation of configurations with the positive fermionic weight |Re det D(µ)| where D(µ) is the Dirac matrix and the signs sign(Re det D(µ)) = ±1 are handled by a discrete reweighting.
As an application of the method we perform a study of the conjectured critical end point in the µ − T phase diagram. At µ = 0 QCD has a cross-over thermal transition and it is expected that as µ is increased the transition gets stronger and eventually at some µ = µ c it becomes a second order phase transition, beyond which at µ > µ c the transition is first order. We would like to unambiguously observe this strengthening of the transition in a manner which is free of uncontrolled systematic uncertainties. The present paper will be limited to the unimproved staggered discretization at fixed N t = 4 hence we do not claim to arrive at continuum results. Nonetheless even at fixed N t the lattice system, as a well defined statistical physics system, may or may not possess a critical end point. This latter question is the one we attempt to address in our paper. Note that the mere idea of using |Re det D(µ)| as a positive weight to generate configurations is not new [53][54][55]. Actual numerical simulations with this method were nevertheless only carried out in the canonical approach in the past [56][57][58].

JHEP05(2020)088
The organization of the paper is as follows. In section 2 we formulate the relevant path integrals in the presence of a chemical potential and reorganize them in a form which allows for a numerical simulation. We present our numerical results in section 3 including our Monte-Carlo algorithm directly at non-zero µ as well as the details of our analysis of the leading Fisher (sometimes also called Lee-Yang) zeros of the partition function. The volume scaling of the leading Fisher zeros is used to infer the order of the phase transition at any given non-zero µ. Finally in section 4 we end with some conclusions and outlook for future work.
2 Path integral at finite µ At finite chemical potential the partition function and expectation values are computed as, where D(U, µ) is the fermionic Dirac matrix involving all flavors and mass terms and S g (U ) is the gauge action. As is well-known det D(U, µ) is complex for real µ = 0, but Z(µ) is nonetheless real. Hence we may equivalently write, It is worth emphasizing that taking the real part above is exact and does not introduce any approximation, as Z(µ) in (2.1) and (2.2) are exactly identical if charge conjugation invariance holds. For a large class of observables we may further write, or if the observable is related to derivatives of Z(µ) with respect to a real µ or mass, etc. In this work we will only be concerned with observables of this type and (2.3) will hold. Although the weights are real now the sign problem of course persists as they can be negative. However one may split the sign ε(U, µ) = sign Re det D(U, µ) of the weights from their absolute values and arrive at Clearly, |Re det D(U, µ)|e −Sg(U ) is positive and can be used as a weight in importance sampling. Configurations will be generated using this weight and the corresponding expectation values will be denoted by . . . abs,µ . The signs ε(U, µ) = ±1 will be dealt with by a discrete reweighting, leading to which is meaningful if the denominator is non-zero. Furthermore, if indeed the denominator is non-zero then the result is trustworthy as there cannot be any overlap problem, since the only reweighting we need to deal with is a reweighting with respect to a discrete set; there are only two sectors, those with ε(U, µ) = +1 and −1. The sign problem is of course still present and it will be signified by the denominator being zero within errors. To summarize the above, what we have achieved by the formulation (2.5) is that the only problem is the sign problem, there is no other uncontrolled systematic which may spoil the result even when the sign problem is not prohibitively severe. Consequently, we have both a sufficient and necessary condition for the correctness of the results: if at a given set of parameters and lattice volumes ε(U, µ) abs,µ is consistent with zero within statistical uncertainties then we have no result, if on the other hand it is non-zero then whatever the result is, it is reliable with well-defined statistical errors.
Let us denote the sets of configurations with ε(U, µ) = ±1 by U ± (µ), which of course depend on µ. Then we have, where the inequalities are meant as exact results at infinite statistics while at finite statistics the left hand sides may of course be consistent with zero within errors.

Numerical results
In our simulations we employ the Wilson plaquette gauge action and 2 + 1 flavors of rooted staggered (unimproved) fermions on N t = 4 lattices. The spatial lattice sizes are L/a = 8, 10, 12 and the fermion masses are set to their physical values am ud = 0.0092 and am s = 0.25. The chemical potential is introduced for the light quarks only, µ u = µ d = µ and µ s = 0 is set for the strange. The setup is identical to [42]. At each µ and spatial volume the bare coupling was set to β c as follows. For each µ, initial β c0 values were taken from [59]. The leading Fisher zeros (see section 3.2), β 1 + iβ 2 were measured in shorter runs and β c0 was modified by ∆β = β 1 − β c0 if necessary. Then all further production runs were performed at these β c = β c0 + ∆β. The resulting values are shown in table 1. From [59] we also glean that the spatial volume dependence of β c is rather mild and in this first exploratory work we set the same bare coupling for all of our 3 spatial volumes.

JHEP05(2020)088
3.1 Monte-Carlo with µ > 0 We would like to generate configurations with the weight |Re det D(U, µ)|e −Sg(U ) . This is a non-trivial problem and to our knowledge no pseudo-fermion type construction can be found. What one may still do is rewrite the weight as since det D(U, 0) is real and positive, and utilize a standard (R)HMC algorithm at µ = 0 and include the µ-dependent ratio in the Metropolis accept/reject step at the end of the trajectory. This will clearly be an expensive algorithm because the full determinant needs to be computed, but with the help of the reduced matrix construction the cost is still manageable for the lattice volumes we will consider in this paper.
Since we are working with rooted staggered fermions we need to compute the full determinant, its square root, its real part and then its sign and absolute value. At finite temperature these steps can most easily be done with the help of the reduced matrix [38]. This has 6(L/a) 3 eigenvalues, λ i (U ), and the main utility of them is that the full staggered determinant can be given at finite chemical potential as, For the precise definition of the reduced matrix see [38]. We will define the square root branch factor-by-factor in the above product by requiring continuity in µ or in other words by requiring that in the µ → 0 limit each factor below goes to unity, The branch cut of the square root is placed on the negative real axis. This procedure fully fixes the complex determinant ratio. The real part, sign and absolute value can then be taken straightforwardly. Notice that with this procedure the partition function remains real since det D(U * , µ) = det D(U, µ) * for real µ, and so our approach maintains its validity. Clearly, if µ is small the ratio |Re det D(U, µ)|/|Re det D(U, 0)| is close to unity and hence will not affect the Metropolis step much, i.e., a tuned (R)HMC algorithm at µ = 0 will perform just as well. On a given spatial volume as µ increases the ratio will influence the Metropolis step more and more and will decrease the acceptance rate. This can be compensated by employing shorter (R)HMC trajectories as this will change the links less and consequently the change in the ratio with respect to the beginning and end of the trajectory will decrease. In this way we are able to keep the acceptance rate above 50% for all runs. The shorter trajectories will of course lead to larger autocorrelation times. Concretely, our estimate of integrated autocorrelation times of our key observable (3.5) are between 50 and 500 depending on µ and L/a. The total number of configurations are between 5 · 10 4 and 2 · 10 5 , leading to a few hundred independent configurations for each simulation point. We observe that "tunnelling" between the +1 and −1 sectors are JHEP05(2020)088 frequent, i.e., the change in the µ-dependent ratio is small enough so that even if the trajectory changes sector we observe good acceptance. The crucial measure of whether the results are reliable or not is given by ε abs,µ , i.e., the average sign, which at the same time measures the strength of the sign problem itself. Since ε abs,µ → 1 as µ → 0 we parametrize it as, with the 4-volume V = L 3 /T and show in figure 1 both ε abs,µ as well as f (µ, V ). Clearly, f (µ, V ) depends mildly on V but does depend non-trivially on µ. As can be seen the volumes L/a = 8, 10, 12 and chemical potentials aµ ≤ 0.2 are safely in the region where ε abs,µ is several standard deviations away from zero, hence the sign reweighting (2.5) can be performed without issues. In particular, as emphasized, there is no overlap problem to contend with. It is worth exploring what the effect of the sign reweighting is on some observables, more precisely how different some observables are in the +1 and −1 sectors. As an example we show the gauge action per unit space time volume in figure 2 as a function of the chemical potential separately for the two sectors.

Fisher zeros
Once it has been determined which volumes and chemical potentials allow for the application of the sign reweighting (2.5) we are able to compute observables. Since our primary interest is the order of the phase transition as a function of µ we will compute the Fisher zeros in the bare coupling β, i.e., we will look for complex bare couplings such that Z(µ, β) = 0 at given µ and volume; see (2.6). This amounts to measuring the observables   several zeros as a function of complex β, we will be looking for the one closest to the real axis, which in every run happens to coincide with the one closest to (β c , 0) in the complex plane as well. This zero will be called the leading zero.
The volume scaling of Im β determines the order of the transition: if Im β → const as L/a → ∞ the transition is a cross-over, if Im β ∼ a 3 /L 3 the transition is first order and finally if Im β ∼ (a/L) α with a non-trivial exponent α > 0 the transition is second order. Although these leading order expressions are unambiguous in all three cases, the subleading terms are a priori not known. Since we know that at µ = 0 the transition is a cross-over for physical quark masses, it is generally expected that for small µ it will stay a cross-over. At fixed µ > 0 the imaginary part of the leading Fisher zero is then extrapolated to infinite volume via, where the exponent 3 in the subleading term is merely an ansatz. In this first study we only simulated at 3 volumes, L/a = 8, 10, 12 and hence we are unable to fit the exponent of the subleading term simultaneously with A and B. Empirically, we do find that the above fit function provides acceptable statistical fits for our choice of chemical potentials.
JHEP05(2020)088  The existence of a critical end point would suggest that Im β ∞ (µ) = A(µ) is a decreasing function of µ and as µ → µ c we have A(µ) → 0.
The real part of the leading Fisher zero on the other hand may be used to define the critical coupling. The simulations we performed at particular values of β c = β c (µ) and we have checked that for aµ > 0.1 the differences ∆β = Re β − β c are deviating from zero less than 3σ and rarely beyond 1.5σ. Note that a smooth cross-over means that different observables may lead to different definitions of the pseudo-critical coupling.
The measured imaginary parts of the Fisher zeros are shown in the left panel of figure 3. The extrapolations to infinite volume using L/a = 8, 10, 12 are shown in figure 4 together with the resulting χ 2 /dof values of the fits. Out of the 10 extrapolations the largest χ 2 /dof values are at aµ = 0.1, 0.1875, 0.2 and are 4.3, 2.37, 2.95. Note that dof = 1 and even the largest χ 2 /dof = 4.3 leads to a q-value of 4%. The resulting Im (β ∞ (µ)) as a function of µ is finally shown in the right panel of figure 3.
The most important result from our investigation can be gleaned from figure 3. Both at finite volumes and correspondingly in infinite volume the imaginary part of the relevant Fisher zero is decreasing as the chemical potential becomes sufficiently large. More precisely, the infinite volume extrapolated result shows that the imaginary part of the leading Fisher zero is more or less flat up to aµ ∼ 0.15 and a sharp decrease is observed for 0.15 ≤ aµ ≤ 0.2. The observed flatness agrees within errors with the slight increase seen in [42,59], and cannot be significantly distinguished from it with the currently available statistics. This means that in the range of chemical potentials where our results are reliable with trustworthy statistical errors, i.e., ε abs,µ = 0, we are able to conclude with high statistical significance that the leading singularity of log Z is eventually moving closer and closer to the real axis. In fact, the location of the singularity is consistent with a real value at aµ ∼ 0.2. This in turn means that the strength of the transition is eventually increasing and very suggestive that a true phase transition occurs at around aµ ∼ 0.2. This corresponds to µ/T c ∼ 0.8, in agreement with [42], however the latter result should be interpreted with caution since, as we explained, we do not expect the fixed N t = 4 results JHEP05(2020)088

JHEP05(2020)088
to be particularly close to the continuum with our chosen discretization. Nonetheless our results are trustworthy in the well defined statistical model given by the N t = 4 lattice system.

Conclusion and outlook
In this paper we have introduced a new technique for evaluating the path integral at finite baryon chemical potential. The approach involves generating configurations by the absolute value of the real part of the fermionic determinant and taking the signs into account by a discrete reweighting. The first step necessitates the evaluation of the full determinant during the Monte-Carlo simulation which makes the algorithm rather costly but still manageable for 8 3 ×4, 10 3 ×4 and 12 3 ×4 which are the volumes we used. The second step, the discrete reweighting by the sign of the real part of the fermionic determinant, is a fully controlled step provided the average sign is several standard deviations away from zero, i.e., the sign problem is not too severe. This requirement can be easily monitored and once it is fulfilled, the results are completely trustworthy with well-defined statistical errors. This feature is the main advantage of our method. It improves on traditional reweighting in µ and/or some other continuous parameter because in that case the notorious overlap problem may invalidate the results even though a naive application of the reweighting formula O µ = Ow 0 / w 0 seemingly presents no problems. Since our main interest was the order of the thermal phase transition as a function of the chemical potential, we have determined the first few Fisher zeros and the volume scaling of the leading one (the one closest to the real axis) for 10 choices of µ in the range 0.025 ≤ aµ ≤ 0.2. We have observed that the strength of the phase transition stays flat within errors for 0 < aµ < 0.15 and increases sharply for 0.15 < aµ < 0.2, signified by the decrease in the imaginary part of the leading Fisher zero. The infinite volume extrapolation of the leading Fisher zero at aµ ∼ 0.2, corresponding to µ B /T ∼ 2.4, is in fact consistent with a true phase transition, i.e., the imaginary part is consistent with zero.
There are however several avenues to improve on our work in the future. First, we have performed simulations at fixed N t = 4, i.e., we have not addressed the continuum limit at all; simulations with larger temporal extents would be necessary in order to do so. Once the continuum behavior is investigated it might be worthwhile to use an improved action, both for the gauge and fermionic actions. In the present work we have used the Wilson plaquette gauge action and unimproved staggered fermions. The motivation was to replicate the setup of [42] where the critical end point was investigated using traditional reweighting in µ. It is worth noting that even though the unimproved staggered discretization on N t = 4 lattices is far from the continuum, it is a well-defined lattice statistical physics model with a sign problem. Hence it makes perfect sense to study it in order to gain valuable insight into the sign problem in general.
Second, the volume scaling of the Fisher zeros is of central importance and our ansatz (3.6) was simply motivated by empirically observing good statistical fits as well as the fact that we only had data on 3 volumes. Hence we were unable to fit all 3 param-JHEP05(2020)088 eters, A, B and C in the general form, Im β = A + B (a/L) C , (4.1) which would otherwise be the justified procedure. Once an additional volume 14 3 × 4 is added, the exponent C could be determined or at least constrained. Third, we have set the quark masses to their physical values at β = β c at µ = 0 and have not changed them for µ > 0 along the line of constant physics. Even though the effect is expected to be negligible relative to other sources of errors, in future work we do plan to follow the line of constant physics for µ > 0.
Fourth, we have included the chemical potential at the quark level as µ u = µ d = µ B /3 = µ and µ s = 0 which corresponds to µ S = µ B /3. Nonetheless our method can be trivially modified to include other chemical potential assignments, e.g., strangeness neutrality S = 0 or µ S = 0.
Finally we mention that the recently introduced geometric matching procedure [43] provides a new rooting procedure at finite N t which is nonetheless expected to agree with the one followed in this paper towards the continuum limit. We repeated the determination of the leading Fisher zeros using geometric matching and found that they agree with the ones presented in this paper within statistical uncertainties. This type of cross-check will be especially useful for future studies targeting the continuum limit.