Dual Formulation and Phase Diagram of Lattice QCD in the Strong Coupling Regime

We present the computation of invariants that arise in the strong coupling expansion of lattice QCD. These invariants are needed for Monte Carlo simulations of Lattice QCD with staggered fermions in a dual, color singlet representation. This formulation is in particular useful to tame the finite density sign problem. The gauge integrals in this limiting case $\beta\rightarrow 0$ are well known, but the gauge integrals needed to study the gauge corrections are more involved. We discuss a method to evaluate such integrals. The phase boundary of lattice QCD for staggered fermions in the $\mu_B-T$ plane has been established in the strong coupling limit. We present numerical simulations away from the strong coupling limit, taking into account the higher order gauge corrections via plaquette occupation numbers. This allows to study the nuclear and chiral transition as a function of $\beta$.


Introduction
The finite baryon density sign problem in lattice QCD hinders a direct evaluation of the phase structure of QCD in the µ B − T plane. In particular, the existence of a critical endpoint (CEP) that is sought for in heavy ion collision experiments at RHIC and LHC could not be established yet via lattice simulations. Although the well established methods for small µ B /T , such as Taylor expansion, reweighting and analytic continuation from imaginary chemical potential can in principle make statements about the existence of the CEP, it is likely that the CEP, if it exists, has a quite large µ crit B , such that it is not within reach with the aforementioned methods.
In recent years, many alternative methods have been proposed and tested to circumvent the finite density sign problem. Most notably, the complex Lagenvin method together with gauge cooling could address full QCD in the deconfined phase [1,2]. Another method based on complexified QCD, the Lefschetz thimbles, are currently applied to QCD-inspired models with few degrees of freedom, but the method is far from being applicable to full lattice QCD [3][4][5].
A promising alternative strategy is to change the degrees of freedom of the original partition function. Since the sign problem is representation dependent, it may be possible to find a different set of variables that are closer to the true eigenstates of the Hamiltonian. Finding such a basis would reduce the sign problem significantly, or even solve it. Changing the degrees of freedom can be for example obtained by a Hubbard-Stratonovich transformation [6], or by introducing auxiliary fields. Another way is to integrate out some of the degrees of freedom to obtain a "dual" representation in terms of world lines. This strategy has been successfully applied to address sign problems in models with an abelian gauge group (such as the massless Schwinger model [7], and the gauge-Higgs models [8]) It is however quite non-trivial to find a dual representation for non-abelian gauge groups. A recent attempt is to decompose the non-abelian components into abelian "color cycles" [9].
Our attempt to perform Monte Carlo simulations on the QCD phase diagram is based on the strong coupling expansion. The starting point is the well-established partition function of staggered fermions in the strong coupling limit. Here, the phase diagram is well established. We then propose a dual representation in terms of world lines and world sheets that incorporates some contributions of the gauge action. For small β, we are able to determine the phase boundary between the chirally broken and chirally restored phase. The leading order correction has been addressed via reweighting from the strong coupling ensemble to β > 0 in [10]. We go beyond this scope by directly sampling the partition function including next to leading order gauge corrections.

Lattice Action and Partition Function
We consider the standard lattice action for staggered fermions (no rooting, no improvement) together with the Wilson gauge action: with a t µ = 1 Nc a t µ B the quark chemical potential. The only modification is that we introduced a bare anisotropy γ, favoring temporal fermion hoppings over spatial fermion hoppings, giving rise to an anisotropy of the lattice spacings a at = ξ(γ). This will allow us later to vary the temperature continuously in the strong coupling regime.
The standard approach for lattice simulations is to integrate out the Grassmann-valued staggered fermions χ andχ to obtain the fermion determinant. However, the fermion determinant becomes complex for finite quark chemical potential, resulting in the finite density sign problem. Our strategy is to expand the action S = S F + S G both in the fermion hoppings and in β = 2Nc g 2 . Then we exchange the order of integration, i.e. integrate out the link variables analytically first, and afterwards the Grassmann variables. The remaining degrees of freedom will be color singlets on the links, and the plaquette occupation numbers n P (from the moments of the fundamental plaquettes tr[U P ] n P ) andn P (from the moments of the anti-fundamental plaquettes tr[U † P ]n P ).
The fermions can be gathered into matrices All elementary plaquettes P from the expansion of S G that share a given link U µ (x) need to be taken into account when integrating out the link U ≡ U µ (x): with P + the subsets of plaquettes in forward and P − in backward direction, as illustrated in Fig. 1. Hence the one-link integral over gauge group G = SU(N c ), U (N c ) that we will consider has the fermion matrices M, M † and the set of staples S P with U P = U µ (x)S P as external sources: where we expand in the forward hoppings κ, backward hoppingsκ, and plaquette and antiplaquette occupation numbers n P ,n P . In the last line, we have decomposed the traces to separate the staples from the gauge link and summation over the set of indices i, j, k ,l is implied. It is the tensor C(β, {S P , S † P }) j i ,l k which leads to non-local color contractions and can be related to the set of plaquette occupation numbers {n P ,n P } when contracting the m open color indices from U andm open color indices from U † with the one-link integrals from the neighbor links: The remaining integral can be related to integrals over the link matrices only [12]: Here, a = κ + m and b =κ +m is the number of U -matrix and U † -matrix elements. In this one-link integral, only the color indices from the quark matrices will be contracted. The contraction of the remaining indices can in general not be carried out easily, however in certain cases link integration on the complete lattice will be possible to give rise to a color singlet partition function: where the admissible graphs G are such that they fulfill the constraint .

Link Integration in the Strong Coupling Limit
For β = 0, link integration factorizes: The corresponding one-link integrals K 0 will not depend on any external gauge links: Hence, link integration can be carried out analytically. Only a finite number of integrals have to be evaluated due to the Grassmann nature of the fermions: since they come in N c colors, 0 ≤ κ ,κ ≤ N c . Moreover, integral Eq. (7) will only be non-zero if κ −κ = qN c with q = 0, ±1 (see next section). The corresponding result for Eq. (5) was first addressed in [13] when deriving the strong coupling partition function for N f = 1: Here, M x =χ x χ x are the mesonic and i Nc are the baryonic degrees of freedom. After the final Grassmann integration, where also the expansion of e 2amqχχ enters, the partition function is exactly rewritten in terms of integer variables: where k b ∈ {0, . . . , N c } are the so-called dimers, i.e. multiplicities of bonds b that represent meson hoppings, n x ∈ {0, . . . , N c } are the so-called monomers and representχ x χ x not being part of dimers, and the baryon world lines form oriented self-avoiding loops, with loop weight Here, N 0, is the number of temporal baryon segments on . The sign σ( ) of a baryon loop is due to geometry: number of backward directions N −, , winding number r and staggered phases along the loop. The sign of a configuration is the product of the signs of all baryonic loops. The sign problem of sampling this partition function is however very mild for any value of the chemical potential, because the baryons are heavy and hence tend to have simple geometries which contribute with positive signs.

Weingarten Functions
In order to obtain the partition function away from the strong coupling limit, we will make use of Weingarten functions [14,15]. This is particularly useful since when some of the link matrices emerge from the Wilson gauge action, we also need contributions to Eq. (7) for a > 0 and b > 0. For n ≡ a = b, the result is expressed via permutations σ, τ ∈ S n on the color indices that go into 2n Kronecker deltas, and are multiplied by the Weingarten functions, which sums over all irreducible representations (irreps) λ of SU(N c ) that are tensors of n fundamental irreps: with D λ (N ) the dimension of the irrep λ of U(N ) and f λ the dimension of the irrep λ of S n . The irreps of both the unitary and the symmetric groups are labeled by integer partitions and due to the finite number of available color indices, the corresponding partitions have a finite number of parts, l(λ) ≤ N . The Weingarten functions contain the character χ ρ λ of the symmetric group S n , which only depends on the conjugacy class ρ = [π] of a permutation π ∈ S n , given by the cycle structure of π. The conjugacy class ρ n is also labeled by an integer partition. Some examples of Weingarten functions are: .
The Weingarten functions for a − b = qN c with q = 1 has been addressed in [16]. For q 0, also epsilon tensors enter Eq. (17), which leads to lengthy expressions. The generalization for q > 1 will be addressed in a forthcoming publication. Here we simply want to illustrate that we recover the strong coupling limit and the leading order gauge correction within this formalism.

Link Integration via Weingarten Functions
The Weingarten functions are a powerful tool to address gauge corrections and integrals for many flavors, given that the matrices M, M † are generalized to N f > 1. Depending on how many fermion hoppings contribute to the link integral, we restrict the sum over irreps within the Weingarten function to those consistent with the fermion content: This restriction is possible due to the orthogonality of characters: for any λ [n] (i.e. with the exception of the completely symmetric irrep which has χ ρ [n] = 1 for all ρ) it holds that with h ρ the number of elements in the conjugacy class ρ. However, due to the additional minus signs from the ordering of the Grassmann variables, there are other irreps λ ∈ Λ which are non-zero. At strong coupling, where all sources are fermionic, only the completely anti-symmetric irrep is non-zero, This agrees for N f = 1 with the result in Eq. (13) For N f = 1 this reproduces the known result [10,11]: With this result, one address gauge corrections as shown in Fig. 3. Similarly, other gauge corrections can be addressed, which we plan to do in a forthcoming publication.

Grassmann Integration
Given that all link integrals K i j ,k l are computed, the remaining task is to organize the fermions such that they can be integrated out. If the integrals K i j ,k l have more than two open indices, the Grassmann integration gives rise to a tensor network that is difficult to evaluate. For U(N c ) gauge theory, the contractions are however possible as there are exact cancellations, as shown in Fig. 3. Only integrals with two open indices K 1,0 give nonzero contributions. Since Grassmann integration results in one incoming and one outgoing loop per site if the site is on the boundary of a plaquette surface, these have to be contracted along loops. The resulting simplifying constraint, exact for U(N c ) and valid for the q = 0 sector of SU(N c ), is that plaquette surfaces are bound by quarks which form self-avoiding loops. However, for q 0, quark loops can intersect such that the constraint is no longer valid. We will nevertheless apply this constraint, resulting in systematic errors for fermionic observables at O(β Nc ).

The Partition Function
With the above simplification, the resulting partition sum is a sum over monomers, dimers, world lines and world sheets defined as surfaces of constant plaquette occupation numbers. To do so, we have to introduce two auxiliary variables which are completely determined by the plaquette configuration: where f b counts the number of fermion fluxes through a bond b, f x counts the number of fermion fluxes through a site x, and f are the self-avoiding loops that are defined on the boundary of the plaquette surfaces of constant plaquette occupation numbers. With this, the partition function reads n P +n P n P !n P ! (30) Figure 3. Simplification due to Grassmann integration within the U(Nc) sector due to exact cancellations. Top: plaquette configurations that result in fx > 1. Bottom: plaquette configurations that result in f b > 1. Hence, in U(Nc), the quark fluxes around the plaquette surfaces form self-avoiding loops. We will also restrict to that in SU(Nc), where this simplification is no longer applicable and introduces systematic errors at O(β Nc ).
Due to restriction discussed Sec. 3.1, we however only sample plaquette surfaces where either n P = 0 or n P = 0, resulting in a net plaquette occupation number n P − n P = 0 ∈ Z.
The color constraint, a modification of the Grassmann constraint, is The N c -flux loops Nc have the same role as baryon loops at strong coupling, but they are now not necessarily made up of N c quarks. Likewise, also dimers are not necessarily mesons, but can be composed of a quark-gluon combination. The bond weights are modified in case a bond is both part of a loop Nc and a loop f : with B 1 a N c -flux bond without and B 2 with an additional dimer. Also the site weights are modified in case fermion flux is reoriented, i.e. when f x = 1, with v 1 = (N c − 1)! the weight when it merges into a dimer, and v 2 = N c ! when it merges with a N c -flux. We sample the partition function Eq. (31) by extending the mesonic and baryonic worm algorithm used at strong coupling. In particular, we update the plaquette occupation numbers on closed loop configurations, and the 0-flux and N c -flux worms take modified weights on edges with f b 0. A detailed discussion of the algorithm will be left for a forthcoming publication. Figure 4. Typical 2-dimensional configuration at finte β, atµ and amq. Left: degrees of freedom that are sampled: monomers (blue), dimers (black), 3-fluxes (red) and plaquette occupation numbers (green). Right: the same configuration but with the substructure of color singlets and triplets along excited plaquettes: quarks (red) and gauge fluxes (green). Baryons becomes extended objects.

Sign Problem
Although the finite density sign problem has been made very mild in the strong coupling limit, this is not necessarily the case away from the strong coupling limit, as fermion hoppings on the boundary of plaquette surfaces take place. Single fermion hoppings are however not suppressed by a large mass. In fact, the sign problem in the dual representation due to finite β even arises for the U(N c ) gauge theory, which is sign problem-free in the conventional fermion determinant representation, as the depenence on the chemical potential drops out. The sign of a configuration factorizes in the N c -flux sign and the fermion flux sign: For N c = 3, the combination of fermion loops and 3-flux loops lead to the following identification, as shown in Fig. 4: dimers on bonds with fermion f b 0 are fermionic, whereas 3-fluxes on bonds with fermion f b 0 are bosonic.
The example of a negative configuration, Fig. 5 (right), illustrates that in two dimensions, negative contributions are related to frustration of monomers: a loop trapping an odd number of monomers has negative sign. This is known from the dual representation of the Schwinger model at finite quark mass. But for dimensions d > 2, even without monomers, a sign problem is induced as dimers and N c -fluxes can be perpendicular on a plaquette surface, giving rise to topologically inequivalent configurations with opposite signs.

Crosschecks
We have made extensive crosschecks on small 2-dimensional volumes where exact enumeration is possible. In Fig. 6 some gauge observables, the average plaquette and the Polyakov loop, are shown as obtained from the dual representation, as a function of am q for µ = 0, and for various gauge groups. They agree well both with the exact result and with hybrid Monte Carlo (HMC). Another important crosscheck where HMC and Meanfield results [19] are available is the phase boundary in the β-T plane for SU(3) at µ = 0. Fig. 7 shows that the results from direct sampling agree well with extrapolations of HMC.

Strong Coupling Regime at Finite Temperature
We have derived the dual representation in the strong coupling limit by taking into account the bare anisotropy γ in order to continuously vary the temperature independent of β. In a recent publication [17], one of us has determined with collaborators the non-perturbative anisotropy a/a t as a function of the bare anisotropy in order to unambiguously define the temperature: We adopt this non-perturbative definition of the temperature, which differs significantly from the previously used mean field result aT = γ 2 Nt . Likewise we convert the chemical potential:

Phase Diagram in the Strong Coupling Regime
Lattice QCD with staggered fermions has a residual chiral symmetry even in the strong coupling regime, since there is an exact Goldstone mode in the spin⊗taste basis γ 5 ⊗ γ 5 . The lattice action at zero quark mass, and likewise partition function Eq. (31) has the symmetry i.e. even and odd sites transform independently. The chiral symmetry is spontaneously broken at low temperatures, but restored at some phase boundary aT c (aµ B ). The transition in the chiral limit is second order for small and intermediate aµ B and turns into a first order transition at low temperatures, separated by a tri-critical point. This point at (aµ tric B , aT tric ) = (1.56(4), 0.73(4)) turns into a critical end point as soon as the quark mass becomes finite. The ratio µ CEP B /T CEP > 2 becomes even larger as a function of the quark mass. The phase boundary for the chiral transition in the strong coupling regime can be measured by finite size scaling of the chiral susceptibility, as shown in Fig. 9 (top). The nuclear transition can be obtained from the position of the gap in the baryon density. In the strong coupling limit, the first order chiral and nuclear transition coincide. The reason is that the nuclear liquid phase is actually a Pauli saturated phase of a baryon crystal, such that no quarks are left for the formation of a chiral condensate. This finding seems to be independent of the quark mass [18]. We restrict in the following to the chiral limit, where simulation via the Worm algorithm are even faster than with finite quark mass, in contrast to HMC.
Via reweighting from the β = 0 ensemble, Fig. 8 (left), it was found that the chiral transition aT c (aµ B ) for small chemical potential indeed decreases, as expected since the lattice spacing a(β) becomes smaller. However, the chiral and nuclear first order transition still coincide with reweighting and mean field theory. This results makes use of the mean field value of a/at = γ 2 for better comparison. The direct simulations favor the scenario of extrapolating the phase boundary via an exponential ansatz (right) rather than a linear ansatz (left), as has been discussed in [10].
with the strong coupling result for small β. This may be very likely a reweighting artifact, as it is impossible to reweight from one phase to another phase across a first order transition. We only found that the nuclear critical end point separates from the chiral tri-critical point, but does not split from the first order line. The expectation is however that the chiral and nuclear transition split, a possible scenario is shown in Fig. 8 (right). It is however a priori not clear how much µ nuclear c and µ chiral c are separated in nature, and how large β needs to be to observe that splitting.
In order to understand the relation between nuclear and chiral transition, we need to sample the partition function Eq. (31) directly at finite β. With the direct simulations at finite β, based on local plaquette updates together with the worm to update the dimers and 3-flux world lines, we find that the chiral first order transition indeed depends on β, as shown in Fig. 9 (bottom). Our lattices were N s times4 with N s = 4, 6, 8, and for various temperatures and baryon chemical potentials, which suffices to determine the chiral phase boundary quite accurately. These preliminary results still needs to be reconciled with the first order nuclear transition, which requires larger volumes.

Conclusion
We have presented a partition function that includes higher order gauge corrections with the constraint that the plaquette world sheets are bound by fermion loops. Plaquette occupation numbers are in principle unbounded, such that we sample contributions of the gauge action at arbitrarily large order in β. However, due to the complicated non-local structure of the tensors C(β, {S P , S † P }) j i ,l k , it is not yet possible to write down a partition function that is correct for all orders in β. Hence we restrict to the limit where plaquettes form surfaces  (3) in the chiral limit as a function of small β, obtained from reweighting [10] but with the non-perturbative anisotropy a/at to convert to aT and aµB.
Contrary to the expectation, the nuclear and chiral transition did not split, which is likely an artifact from reweighting. Right: one of several possible scenarios on the β-dependence of the chiral and nuclear transition for unrooted staggered fermions in the chiral limit.
bounded by quark flux. This restriction is no longer valid for SU(N c ), and our approximation will result in systematic errors in fermionic observables at O(β Nc ). However, in the strong coupling regime with β 2N c , these systematic errors are expected to be small. Due to the sign problem induced by the boundaries of the plaquette surfaces, simulations are restricted to β 1. We presented first direct measurements at non-zero β and µ, which are consistent with the previous results from reweighting. It will be essential to improve on the sign problem further to apply these methods for β > 1. A systematic error on the phase boundary as shown in Fig. 9 is due to the anisotropy ξ = a at . We only considered the bare anisotropy γ F ≡ γ in the Dirac coupling, but one should also introduce an anisotropy in the Wilson action, γ G = β t /β s . Then the lattice anisotropy is a non-perturbative function of both bare anisotropies, ξ(γ F , γ G ), that can in principle be determined in a similar way as in [17].
In this work we have only studied the gauge corrections of the phase diagram in the chiral limit. We plan to study the gauge corrections also at finite quark mass.

Acknowledgement
We would like to thank Philippe de Forcrand and Hélvio Vairinhos for stimulating discussions. This work is supported by the Emmy Noether Program under the grant UN 370/1-1.
Computations have been carried out on the OCuLUS cluster at PC2 (Universität Paderborn).