The effect of stout smearing on the phase diagram from multiparameter reweigthing in lattice QCD

The phase diagram and the location of the critical endpoint (CEP) of lattice QCD with unimproved staggered fermions on a $N_t=4$ lattice was determined fifteen years ago with the multiparameter reweighting method by studying Fisher zeros. We first reproduce the old result with an exact algorithm (not known at the time) and with statistics larger by an order of magnitude. As an extension of the old analysis we introduce stout smearing in the fermion action in order to reduce the finite lattice spacing effects. First we show that increasing the smearing parameter $\rho$ the crossover at $\mu = 0$ gets weaker, i.e., the leading Fisher zero gets farther away from the real axis. Furthermore as the chemical potential is increased the overlap problem gets severe sooner than in the unimproved case, therefore shrinking the range of applicability of the method. Nevertheless certain qualitative features remain, even after introducing the smearing. Namely, at small chemical potentials the Fisher zeros first get farther away from the real axis and later at around $a\mu _q = 0.1 - 0.15$ they start to get closer, i.e., the crossover first gets weaker and later stronger as a function of $\mu$. However, because of the more severe overlap problem the CEP is out of reach with the smeared action.


I. INTRODUCTION
One of the most interesting open problems in the study of QCD is determining its phase diagram in the temperature (T )-baryonic chemical potential (µ B = 3µ q = 3µ) plane. It is established that near µ = 0 there is an analytic crossover [1,2] at T = 150−160 MeV [3][4][5][6]. It is conjectured that at higher chemical potential this crossover gets stronger and above a point it turns into a first order phase transition. The endpoint of the line of first order transitions is called the critical endpoint.
Among these methods reweighting has no other systematic error besides the overlap problem, thus in principle it can lead to the correct results with infinite statistics. With this method, the location of the critical endpoint has been determined for N t = 4 lattice QCD with an unimproved staggered action [38,39], but there are still no results closer to the continuum limit. This is due to reweighting methods getting prohibitively expensive for large lattices. One possible way to reduce the finite lattice spacing effects without increasing the lattice size is to improve the UV behavior of the action, e.g., by introducing stout smearing [41] in the fermion action.
In this paper we first reproduce the old N t = 4 results with the unimproved action with an exact algorithm [42] (not known at the time of Refs. [38,39]) and much larger statistics. Second, we use one step of stout smearing with several values of the ρ smearing parameter and study the behavior of the phase diagram and the severity of the overlap problem at nonzero chemical potential as a function of ρ. Our simulations use N f = 2 + 1 rooted staggered fermions with physical quark masses. The ambiguity of rooting at finite µ was first discussed in detail in Ref. [43] and also recently emphasized in Ref. [44] where the so-called geometric matching method was introduced to deal with the analyticity issues of the rooted staggered determinant. Here we will use the geometric matching method and compare it to standard rooting.

II. MULTIPARAMETER REWEIGHTING
At finite µ B the fermionic determinant detM is generally complex and importance sampling methods cannot be applied. This problem can be circumvented with the reweighting technique. The main idea is to rewrite the grand canonical partition function as: where for a fixed gauge configuration The configurations can now be generated with importance sampling methods since detM(0) is a positive real number. The weight w(µ, β) is instead treated as an observable. Even though Eq. (1) is exact, in practice its application is limited to small enough values of the chemical potential and of the volume. This is due to the fact that the tails of the distribution of the weights w(µ, β) are hard to sample. This is known as the overlap problem and is exponentially severe in the volume. Reweighting in β was introduced to decrease the severity of this problem, with β 0 chosen to be close to the value of the gauge coupling at the crossover at zero chemical potential in order to have a more variable ensemble [37,38]. The analytical form of the staggered fermion determinant on a fixed gauge configuration can be expressed at arbitrary µ with the eigenvalues λ i of the so called reduced matrix [35]: where N s and N t are the linear spatial and temporal size of the lattice in lattice units, and the λ i depend on the gauge fields. This determinant describes four flavors of quarks. In order to describe N f < 4 the rooting trick is invoked, i.e., detM is replaced by (detM) N f /4 . In particular, for N f = 2 the complex square root is taken. This introduces a sign ambiguity of the rooted determinant.
At µ > 0 one usually chooses the complex root that continuously connects to the positive root at µ = 0 for each factor in Eq. (3), i.e., the first factor in the Eq. (2) becomes: In Ref. [44] a new approach, dubbed geometric matching, has been proposed in order to remove the branch point singularities in Eq. (4) configuration by configuration. This makes the partition function an entire function of µ already at finite lattice spacing, which is not the case when one uses Eq. (4). The new definition of the corresponding reweighting factor reads: where the set of the ξ i are obtained from the set of the λ i via geometric matching [44]. Essentially, they are the geometric means of the close-lying pairs of λ i . In the continuum limit the two methods are expected to give the same results. Phase transitions can be studied by means of the zeros of the partition function [45,46]. The zeros as a function of complex β are called Fisher zeros [47]. A phase transition is signaled by an accumulation of the zeros near the real axis in the thermodynamic limit. In the reweighting approach one solves: on the complex β plane for fixed µ, where w(µ, β) is defined in (2). When using rooted staggered fermions the ratio of determinants in Eq. (2) has to be replaced by the appropriate complex root. In particular, for N f = 2 one can use either of Eqs. (4) or (5). As the volume V goes to infinity the imaginary part of the zero closest to the real axis goes to a non-zero constant for a crossover and to zero in case of a true phase transition. At large enough volumes one expects: where for a first order phase transition b = 1 and Imβ ∞ F = 0, for a second order transition b < 1 and Imβ ∞ F = 0, while in the absence of phase transition Imβ ∞ F = 0 and b is in general not known. In our infinite volume extrapolations we will use Eq. (7) with b = 1 in order to determine whether our data is consistent with a first order phase transition or not. Our estimate of the critical endpoint will correspond to the first value of µ where Imβ ∞ F is consistent with zero. When Imβ ∞ F is nonzero it can be considered as a measure of the strength of the crossover.
In this paper we studied the effect of stout smearing [41] on the Fisher zeros. Stout smearing is an analytic map that replaces link variables with a suitable average in order to smooth ultraviolet fluctuations. The procedure depends on a parameter ρ which measures the relative weight of the original and the neighboring links in the average. This procedure has been shown to reduce cutoff effects on several observables. Since stout smeared links are used to compute the fermionic determinant, the Fisher zeros become functions of the smearing parameter ρ.

A. Lattice action and algorithm
We used the plaquette gauge action and 2+1 flavors of staggered fermions with one step of stout smearing with smearing parameter ρ. For ρ ={0, 0.01, 0.02, 0.03} we have performed simulations near the crossover temperature, which for the different values of ρ correspond to β 0 ={5.188, 5.171, 5.154, 5.137} respectively. The bare quark masses were set to m ud = 0.0092, m s = 0.25. These correspond to the physical Goldstone pion and kaon mass and are the same as in Ref. [39]. In this exploratory study we use the same bare quark masses for all values of β 0 and ρ. We checked that this has at most a 2% effect on both the scale setting and the pion mass. For each ρ we simulated lattices of size 6 3 × 4, 8 3 × 4, 10 3 × 4, 12 3 × 4. As in Ref. [38,39] we only introduce a chemical potential for the light quarks, i.e., we use µ s = 0 which corresponds to µ S = µ B /3.
In Ref. [39] the diagonalization required for Eq. (4) was carried out on ∼ 3000 configurations, which were generated with the R algorithm [48] with a fixed step size. This is not an exact algorithm, thus it may have uncontrolled systematic errors. Here we used the Rational Hybrid Monte Carlo algorithm [42], which is exact. Furthermore since the overlap problem occurs, it is worthwhile to repeat the analysis with much larger statistics. For the diagonalization of the reduced matrix, we applied GPU accelerated determinant calculations using the publicly available MAGMA library [49][50][51]. This allowed us to measure the determinant on an order of magnitude more configurations than in Ref. [39]. The number of configurations for each simulation point is shown in Table  I.

B. Results with the unimproved action
First we show in Fig. 1 a comparison of the position of the leading Fisher zero, i.e., the one closest to β 0 , obtained with the standard rooting and geometric matching for the largest volume in our study at ρ = 0. This assesses the systematics related to the rooting ambiguity. The two methods agree nicely in this case. To estimate the position of the leading Fisher zero we performed an infinite volume extrapolation linear in 1/V excluding the smallest volume for both types of rooting. As discussed in the previous section this corresponds to checking whether the data is consistent with a first order phase transition or not. The results are shown in Fig. 2. As can be seen in Fig. 1 and Fig. 2 the leading Fisher zero as well as our estimate of its infinite volume extrapolation first get farther away from the real axis, i.e., the transition gets weaker at small values of µ. Indeed, on the 12 3 × 4 lattice the slope of the imaginary part of the leading Fisher zero at µ = 0 is positive within 3σ. At larger µ the Fisher zero starts to get closer to the real axis, i.e., the crossover strengthens. The value of the chemical potential where the extrapolated imaginary part starts to be consistent with zero is Here, the errors are purely statistical. This is our estimate of the location of the critical endpoint for N t = 4 unimproved staggered fermions. Beyond this point estimate the results with the two types of rooting start to disagree, i.e., the ambiguity in the rooting procedure starts to matter. As we will discuss below, in this region the overlap problem is already strong and therefore it is somewhat unclear whether the difference between the two methods is a genuine effect or only an artifact of an incorrect sampling of the tail of the distribution of the weights w(µ, β).

C. Results with stout smearing
The infinite volume estimates of the leading Fisher zero position with the different values of the smearing parameter ρ can be seen in Fig. 3. Here we only present results obtained with geometric matching. In the left panel we zoom in on the region around µ = 0 to emphasize two effects. First, for all values of ρ considered in this work at small chemical potentials the crossover first gets weaker with increasing µ. Second, already at µ = 0 the leading Fisher zero gets farther away from the real axis as ρ is increased, i.e., the crossover weakens with increasing ρ. This is probably a manifestation of the critical line of the Columbia plot getting farther away from the physical quark masses [52,53]. In the right panel of Fig. 3 we show the full range of chemical potentials. One can see that some qualitative features of the unsmeared case remain: at intermediate values of µ the transition starts to get stronger. However, the musch larger errors make the determination of the CEP unreliable and we preferred not to attempt it.
As can be seen in the right panel of Fig. 3 even though the statistics are comparable the errorbars at finite smearing are considerable larger than at ρ = 0. This is due to two effects. As we will discuss below the overlap problem gets stronger with finite smearing. It is also known that the strength of the sign problem as measured by the fluctuations of the phase of the determinant is weakened by taste symmetry breaking. For a demonstration of this see Ref. [29]. Reducing taste symmetry breaking via smearing is then expected to make the sign problem worse.
To look more closely at the overlap and sign problems we consider the distribution of the real part of the weights w(µ, β) for the case of β = Re β F (µ), i.e., we consider reweighting along the crossover line as defined by the real part of the location of the leading Fisher zero. These weights read: where P i is the average plaquette of a single configuration and the rescaling factor N is chosen such that log W = 0. In Fig. 4 we show the distributions of the positive and negative weights, obtained using geometric matching, separately on a logarithmic scale at ρ = 0 and 0.01. As µ increases the two distributions become comparable, signaling the onset of the sign problem. At the same time the support of the histogram broadens and configurations with large reweighting factors become more and more likely. In particular configurations with very large reweighting factors appear but only rarely ("outliers"), indicating that our µ = 0 simulations are sampling the tail of the distribution poorly. This signals the onset of the overlap problem. The range spanned by the outliers increases as µ is increased, indicating that the overlap problem gets worse at larger chemical potentials. Moreover, it also increases when smearing is introduced, again indicating a worsening of the overlap problem. In Fig. 5 we compare the distributions of the positive weights obtained with geometric matching and standard rooting at the largest value of µ considered in this paper for different values of the smearing parameter. At ρ = 0 and 0.01 the two distributions agree nicely except in the outlier region. However since the contribution of the outliers is substantial the central value for the Fisher zero can (and does)change appreciably, but the statistical errors are also large, since the jacknife samples without the outliers differ considerably from the central value. At ρ = 0.02 the situation is similar.
It seems then that substantial differences in the Fisher zeros obtained with the two rooting methods appear only when the overlap problem is strong. This is most likely also the case for the unsmeared action at aµ q > 0.2.

IV. SUMMARY AND CONCLUSION
The determination of the phase diagram of lattice QCD at finite baryochemical potential, and in partic- ular that of the critical endpoint, are especially difficult due to the infamous sign problem. Bypassing this problem by extrapolating from µ B = 0, where it is absent, to finite µ B is possible in principle, but hindered by the overlap problem, i.e., the difficulty in sampling correctly the probability distribution of the relevant observables. Such a problem is exponentially severe in the volume of the system. Despite these difficulties, multiparameter reweighting techniques were employed successfully in Refs. [38,39] to determine the CEP using unim-proved staggered fermions on lattices of fixed temporal size N t = 4 with physical quark masses. In this paper we have confirmed this result using an exact algorithm [42] not available at the time (so reducing systematic errors), much higher statistics (so reducing the statistical error), and a new definition [44] of the rooted fermion determinant which improves the analyticity properties of the partition function. The approach used here is based on the study of the large-volume limit of the Fisher zeros of the partition function at finite µ B using multiparameter reweighting techniques. We have found for the critical endpoint µ c B /T = 2.28 ± 0.07, a result compatible with that of Ref. [39].
Extrapolation of this result to the continuum is obviously difficult, in particular due to the need of bigger lattices that would make the overlap problem worse. A possible way to reduce UV effects and bring the system effectively closer to the continuum is to employ stout smearing [41] on the gauge links. In this paper we have studied the effect of one step of smearing with small smearing parameter on the Fisher zeros of the partition function at finite µ B . Unfortunately smearing turns out to make the overlap problem worse, making it appear sooner, i.e., at lower values of µ B , and making the determination of the critical endpoint unreliable. This is probably related to the critical line of the Columbia plot moving away from the physical point when smearing is introduced [52,53]. In fact, for the unimproved action the leading Fisher zero(i.e., the one closest to the real β = β 0 where the simulations are performed) is relatively close to the real axis at µ B = 0, thus allowing a sufficiently accurate sampling of the relevant configurations for the finite-µ B physics using simulations at µ B = 0. This reflects the presence of a genuine phase transition at µ B = 0 for quark masses smaller but not much smaller than the physical ones [52]. On the other hand, it is known that smearing pushes the critical values of the masses away from the physical point [53]. This is expected to lead to a weakening of the transition, which should reflect into the leading Fisher zero moving away from the real axis. This is precisely what we observed in our study: even for small values of the smearing parameter the imaginary part of the leading Fisher zero grows appreciably already at µ B = 0. This leads to a poorer sampling at µ B = 0 of the configurations relevant at finite µ B , and an earlier onset of the overlap problem. This is particularly evident if one looks at the distribution of the reweighting factors used to reconstruct the finite-µ B theory, focusing on the presence of outliers with large reweighting factors.
These are a symptom of the poor sampling of the tails of the distribution of the weights, and become more and more relevant as µ B increases. It turns out that smearing makes the problem worse, with outliers appearing already at smaller µ B .
As a side analysis, we have compared the Fisher zeros obtained using the recently introduced geometric matching method [44], with those obtained with the standard rooting. While the two methods are expected to lead to the same continuum limit, finite lattice spacing effects could be very different. Our results with the two methods are compatible whenever the overlap problem is absent or mild, indicating that they have comparable cut-off effects at zero and small µ B , at least for the lattices considered in this work. Incompatibility of the two methods only shows up when the overlap problem is severe, and it is not possible to determine whether it signals a large difference in their finite lattice spacing effects or whether it has to be ascribed to the overlap problem itself.
In conclusion, the results of this paper indicate that the bottleneck in studying QCD at finite µ B with multiparameter reweighting is the overlap and not the sign problem. Extending our results closer to the continuum limit therefore requires a much better control of the overlap problem.