Strong-Coupling Lattice QCD on Anisotropic Lattices

Anisotropic lattice spacings are mandatory to reach the high temperatures where chiral symmetry is restored in the strong coupling limit of lattice QCD. Here, we propose a simple criterion for the nonperturbative renormalisation of the anisotropy coupling in strongly-coupled SU($N$) or U($N$) lattice QCD with massless staggered fermions. We then compute the renormalised anisotropy, and the strong-coupling analogue of Karsch's coefficients (the running anisotropy), for $N=3$. We achieve high precision by combining diagrammatic Monte Carlo and multi-histogram reweighting techniques. We observe that the mean field prediction in the continuous time limit captures the nonperturbative scaling, but receives a large, previously neglected correction on the unit prefactor. Using our nonperturbative prescription in place of the mean field result, we observe large corrections of the same magnitude to the continuous time limit of the static baryon mass, and of the location of the phase boundary associated with chiral symmetry restoration. In particular, the phase boundary, evaluated on different finite lattices, has a dramatically smaller dependence on the lattice time extent. We also estimate, as a byproduct, the pion decay constant and the chiral condensate of massless SU(3) QCD in the strong coupling limit at zero temperature.


I. INTRODUCTION
For all practical purposes, the sign problem in lattice QCD with staggered fermions at finite density has been solved at strong coupling. By integrating out the gauge degrees of freedom exactly -which allows replacing Grassmann integration by a sum over fermionic colour singlets -the sign problem becomes mild enough to allow for controlled numerical results at moderate volumes, by combining importance sampling and reweighting methods. As a result, the phase diagram of lattice QCD in the strong coupling limit [1] and at first order in the strong coupling expansion [2] can be completely mapped.
In practice, however, it is not sufficient to simulate the strongly-coupled theory directly on rectangular lattices, because the critical temperature of chiral symmetry restoration is higher than what can be reached using the smallest lattice time extent. 1 In order to study the thermodynamical properties of staggered lattice QCD, in particular across the chiral phase transition, it is therefore necessary to simulate the theory on anisotropic lattices.
On anisotropic lattices, one assigns independent lattice spacings to the spatial and temporal directions, respectively a and a t . The corresponding physical extents of the lattice can then be varied continuously, and independently. A more useful parameterisation of the lattice geometry uses the spatial lattice spacing, a, and the anisotropy parameter ξ, which becomes unity when the lattice is isotropic, and diverges in the continuous time limit a t → 0. In this parameterisation, the lattice temperature is given by: where N t is the lattice time extent. Hence, the lattice temperature can be varied continuously, through ξ.
In lattice gauge theory, the physical parameters a and ξ can only be varied implicitly, through independent bare parameters: the bare gauge coupling β and the bare anisotropy coupling γ. These bare parameters couple differently to the spatial and temporal plaquettes in the Wilson action of SU(N c ) or U(N c ) pure lattice gauge theory in d + 1 dimensions [3]: where U x,µν is the ordered product of link variables around a plaquette parallel to theμ andν directions.
For a single flavour of staggered fermions in the strong coupling limit (β = 0), the anisotropic lattice action is given by: where a t m q and a t µ q are the bare quark mass and quark chemical potential, respectively, and η xµ = ±1 are the staggered phases. In the case of U(N c ), gauge invariance dictates that colour singlets are independent of a t µ q , hence we may set a t µ q to zero without loss of generality.
How a and ξ depend on the bare parameters of the theory is unknown a priori. This knowledge is, however, essential for precision measurements on anisotropic lattices, e.g. bulk thermodynamic quantities, and any uncontrolled approximation can easily be the main source of systematic errors.
In the weak gauge coupling regime (β → ∞) of the SU(N c ) pure gauge theory Eq. (3), perturbation theory and the non-renormalisation of the speed of light can be used to calibrate the anisotropy coupling [4]. In that regime, it is found that ξ pert (γ) = γ (as expected classically).
Using mean field techniques, the behaviour of the renormalised anisotropy at strong coupling (β 1) and at large values of γ is predicted to be quadratic, with unit prefactor [5]: In the nonperturbative regime, however, the relation between bare and renormalised anisotropy couplings can only be determined numerically. This has been done, for example, in pure gauge theory [3,6], in lattice QCD with staggered fermions [7] or Wilson fermions [8]. The nonperturbative renormalisation of the bare parameters requires fine-tuning, guided by some physical criterion which controls the recovery of Euclidean symmetry.
In this Letter we present a simple, precise, and nonperturbative method to calibrate the anisotropy coupling in lattice QCD with massless staggered fermions, in the limit of strong gauge coupling.

II. DIAGRAMMATIC REPRESENTATION OF LATTICE QCD
The partition function of SU(N c ) or U(N c ) QCD on a bipartite N t ×N d s lattice, with a single flavor of staggered fermions, in the strong coupling limit (β → 0) factorises into a product of solvable fermionic one-link integrals: In the SU(N c ) case, the group integration of the link variables, followed by the Grassmann integration of the fermionic degrees of freedom, yields the partition function of a monomer-dimer-loop system [9]: This partition function is a constrained sum over in-teger occupation numbers of monomers and dimers, n x , k xµ ∈ {0, 1, . . . , N c }, and of oriented baryon links, xµ ∈ {0, ±1}, which combine to form oriented baryon loops. The global quantities: enumerate the monomers, temporal dimers, and temporal baryon links on the lattice, respectively. σ( ) = ±1 is a geometric sign associated with the configuration of baryon loops , | | is their length, and w( ) is their winding number around the Euclidean time direction. The monomers represent fermion condensates, M nx x , dimers represent meson hoppings, (M x M x+μ ) kxµ , and baryon links represent baryon hoppings, B x B x+μ or B x+μ B x , where M x is a meson and B x is a baryon: In order for a configuration of occupation numbers to contribute non-trivially to the partition function Eq. (7), the Grassmann integrals over the corresponding fermionic degrees of freedom must be non-trivial on each lattice site.
Due to their Grassmann nature, such configurations must necessarily represent arrangements of exactly N c fermions and N c anti-fermions on each lattice site. 2 This imposes the following local constraints on the integer occupation numbers: Eq. (10b) is a local discrete conservation law for baryon links, which formalises our statement above that baryon links in admissible configurations form closed oriented loops. In the U(N c ) case, since xµ = 0, the partition function Eq. (7) reduces to a sum over monomer-dimer configurations: with the same Grassmann constraint for monomers and dimers on each site: Likewise, the U(N c ) observables are defined in the same way as the observables (in the mesonic sector) of the SU(N c ) theory.

III. CONSERVED CURRENTS AND CONSERVED CHARGES
Let σ x = ±1 be the parity of the site x on a bipartite lattice. From Eq. (10b), it is easy to construct baryonic currents: (13) which are conserved at every site: The corresponding conserved charges are integrals of the baryonic currents Eq. (13) over a codimension-1 lattice slice S µ , perpendicular toμ: Similarly, by rewriting Eq. (10a) as: it is easy to construct the corresponding (pion) currents: from which a local discrete Gauss' law for dimers results: Thus, monomers are sources of the pion currents. Using Grassmann variables, the source term on the r.h.s. of Eq. (18) corresponds to −a t m q ψ x γ 5 ψ x . Only in the chiral limit, i.e. in the absence of monomers, are the pion currents conserved. In the chiral limit, the corresponding conserved charges are integrals of the pion currents over a lattice slice S µ : In the U(N c ) theory, since xµ = 0, the pion currents simplify to:

IV. NONPERTURBATIVE ANISOTROPY CALIBRATION
In this Section, we show how the conserved pion charges can be used to calibrate the anisotropy coupling in lattice QCD with staggered fermions, at zero temperature, in the strong coupling limit.
In the strong coupling limit, the partition functions of SU(N c ) and U(N c ) lattice QCD with staggered fermions have monomer-dimer-loop representations, Eqs. (7) and (11), with no dependence on the spatial lattice spacing, a. In order for the pion charges Q µ to be conserved, we take the lattice fermions to be massless, a t m q = 0. In the SU(N c ) case, we only consider the case of zero chemical potential, a t µ q = 0. 3 The corresponding partition functions thus depend only on a single parameter: the bare anisotropy coupling γ.
Let us consider the theories to be defined on anisotropic N t × N d s lattices. In order to calibrate the anisotropy, we compare the fluctuations of the conserved pion charges in different directions.
Due to spatial isotropy, the expectation values of fluctuations of the spatial pion charges Q i , i = 1, . . . , d must coincide. Therefore, it is convenient to quantify spatial fluctuations using the expectation value of: while the temporal fluctuations are quantified using the expectation value of Q 2 t = Q 2 0 . Now, when the lattice is hypercubic, i.e. N t = ξN s , the fluctuations of the spatial and temporal conserved charges must be equal. This provides a simple, nonperturbative criterion for the renormalisation of the anisotropy coupling: the value of the bare parameter, γ np , corresponding to the renormalised value, ξ(γ np ) = N t /N s , is that for which the fluctuations of the spatial and temporal conserved charges are equal: In Fig. 1, we give a practical example. In a numerical simulation of U(3) lattice QCD on a 32 × 16 3 lattice, we evaluate Q 2 s and Q 2 t for a few values of the bare parameter γ, about the correct nonperturbative value γ np associated with the renormalised anisotropy parameter, ξ = 2. Using Ferrenberg-Swendsen multihistogram reweighting, we interpolate the measurements of the fluctuations, and estimate with high precision the value of the bare parameter for which the two curves  intersect, i.e. when the lattice is hypercubic. In this particular case, γ np = 1.55725(29). This value is to be compared with the commonly accepted mean field prediction, γ mf = √ ξ = 2 ≈ 1.41421.

V. RUNNING ANISOTROPY
It is also possible to estimate the running of the anisotropy parameter, 1 ξ dξ dγ , using extra information from the intersection point in Fig. 1. This quantity -the strong-coupling analogue of Karsch's coefficients [4] -is important for computing e.g. bulk thermodynamic quantities, like the energy density and pressure [11].
The fluctuations of the conserved charges scale with the volume of the lattice slices on which the corresponding conserved currents are integrated over: The ratio of temporal and spatial fluctuations then becomes directly related to the renormalised anisotropy: We have already explained the fact that this ratio is 1 when the lattice is hypercubic. Now, taking the derivative of Eq. (24) with respect to the bare parameter γ, at the intersection of the curves in Fig. 1, yields the value of the running anisotropy at that point: Therefore, in order to estimate the value of the running anisotropy at γ np , we also need the value of the fluctuation of the conserved pion charges on a hypercubic lattice: and the values of the slopes of the tangents to the curves at the intersection point: Q 2 t γnp and Q 2 s γnp .

VI. NUMERICAL RENORMALISATION
The Monte Carlo sampling of the U(N c ) partition function Eq. (11) is highly efficient when using directed path algorithms [10]. In the SU(N c ) case, observables must be reweighted because of the occurrence of negative-weight baryonic configurations, even at zero chemical potential. However, this sign problem is mild and controllable for moderate lattice volumes [1,2,11].
We simulate massless U(3) and SU(3) lattice QCD in the strong coupling limit, using the directed path algorithm [10], for several values of the renormalised anisotropy ξ. For each ξ, we estimate the corresponding value of the bare parameter γ np on a (ξN s ) × N 3 s lattice, for several values of N s , using the method described in Section IV. We also measure the running of the anisotropy coupling Eq. (25). The results for U(3) and SU(3) are summarised in Tables I and II, respectively.
In these tables, rather than storing the estimators of Eq. (25), we instead store the estimators of its reciprocal, the reason being that the latter enters linearly in the definition of important bulk thermodynamic quantities, e.g. the energy density: The nonperturbative relation between the renormalised and bare anisotropy parameters, in the thermodynamic limit, is presented in Fig. 2a. At large anisotropies, the renormalised parameter depends quadratically on the bare parameter. Such a behaviour is expected from mean field arguments. However, the corresponding prefactor differs significantly (≈ 25%) from that of the mean field relation Eq. (5). This introduces a significant systematic error in any numerical study of strongly-coupled lattice QCD.  (2) 0.43408 (2) Table I. Values of the bare anisotropy coupling γnp associated with the renormalised anisotropy ξ, from numerical simulations of massless U(3) lattice QCD on (ξNs) × N 3 s lattices. Also, the corresponding values of the running anisotropy (derivative), the helicity modulus a 2 Υ, and the chiral susceptibility density a 6 χ/N 4 s . The quantity γnp exhibits small finite-volume corrections, and is consistent (within errors) with its thermodynamic limit, even on the smallest lattices. This rapid convergence justifies using small-lattice measurements as thermodynamic estimators for γnp. This is particularly useful in simulations at large ξ, for which significant statistics can only be obtained on small volumes.
We find that the whole range of measurements is well described by a simple, one-parameter rational Ansatz (see Fig. 2b): where κ is a constant, and λ ! = κ/(1−κ), from the requirement that ξ(1) ! = 1. The approach to the continuous time limit is better captured by Taylor expanding Eq. (28) to quadratic order in 1/ξ 2 (see Fig. 2c): The fitted values of κ using the Ansatz Eq. (29) -consistent with those obtained using the Ansatz Eq. (28) -  Figure 2. Nonperturbative relation between the bare and renormalised anisotropy parameters, for U(3) (green) and SU(3) (purple) massless lattice QCD, in the thermodynamic limit, presented in 3 different ways. Fig. 2a shows that, as predicted by mean-field, the renormalised anisotropy at large γ is ξ(γ) ∝ γ 2 , but with a smaller prefactor than predicted. Fig. 2b shows the ratio ξ/γ 2 for a wide range of γ, larger and smaller than 1. A simple one-parameter Ansatz Eq. (28) describes the data well. Fig. 2c shows the approach to the continuous time limit, i.e. ξ → ∞. In that regime, ξ/γ 2 approaches a constant κ, with quadratic corrections in 1/ξ 2 . The behaviours of U(3) and SU (3) where errors are statistical only. This prefactor is significantly different from the mean field value 1. Values for U(3) and SU(3) are statistically consistent with each other. This is to be expected: in the continuous time limit, baryons become increasingly static, and their effect on pion currents vanishes at T = 0.
The Ansatz Eq. (28) is also consistent, after differentiation, with the Monte Carlo data for the running anisotropy. In particular, for the isotropic case, instead of the mean field value: we find nonperturbative corrections consistent with:

VII. APPLICATIONS
In this Section, we use the nonperturbative relation between ξ and γ, determined above, in order to control the convergence of several physical quantities to their continuous time limits.
First, we examine the N t -dependence of the phase boundary of the (µ q , T ) phase diagram of massless SU(3) lattice QCD, and its sensitivity to the anisotropy prescription. Then, we estimate the continuous time values of the static baryon mass am B , the pion decay constant aF π , and the infinite-volume chiral condensate a 3 Σ, in massless U(3) or SU(3) lattice QCD. We use a quadratic Ansatz in 1/ξ 2 , consistent with O(a 2 ) discretization errors of staggered fermions, to model the anisotropy corrections to the continuous time limit: where O is one of the physical quantities listed above, and O CT is the corresponding continuous time value. Table II. Values of the bare anisotropy coupling γnp associated with the renormalised anisotropy ξ, from numerical simulations of massless SU(3) lattice QCD on (ξNs) × N 3 s lattices. Also, the corresponding values of the running anisotropy (derivative), the helicity modulus a 2 Υ, the chiral susceptibility density a 6 χ/N 4 s , and the average baryonic sign. The quantity γnp exhibits small finite-volume corrections, and is consistent (within errors) with its thermodynamic limit, even on the smallest lattices. This rapid convergence justifies using small-lattice measurements as thermodynamic estimators for γnp. This is particularly useful in simulations at large ξ, for which significant statistics can only be obtained on small volumes. In the continuous time limit, the baryons become increasingly static, which explains the lack of fluctuations that contribute to the sign problem at large ξ. For example, the pion decay constant at T = 0 can be shown to be related to the helicity modulus Υ [13]: which corresponds, in the diagrammatic representation, to the variance of the conserved pion charges Q µ on a hypercubic lattice [12]: In turn, the chiral condensate Σ at T = 0 can be estimated from the finite-size scaling of the chiral susceptibility χ, evaluated on hypercubic lattices. This has been done in d = 3 + 1 at finite temperature [12]. In our case where T = 0, chiral perturbation theory of the O(2) model predicts the leading finite-size corrections to be of the form [13]: where β 1 = 0.140461 and α is given by: with β 2 = −0.020305, and Λ Σ , Λ M are renormalisation group invariant scales. The average value of the chiral susceptibility is estimated using intermediate configurations -generated with the directed path algorithm -which sample the mesonic two-point function, as described in [10].

A. Phase diagram
An example of a study that is sensitive to the choice of an anisotropy prescription is the mapping of the phase diagram of massless SU(3) lattice QCD, in the strong coupling limit [1].
The phase boundary separating the chirally broken phase at low (µ q , T ) and the chirally symmetric phase at high (µ q , T ) is determined by monitoring the chiral condensate a 3 Σ during Monte Carlo simulations, using directed path algorithms and sign reweighting for importance sampling on moderate volumes (see Fig. 3).
For fixed N t , the temperature is varied implicitly through the bare coupling γ [2]. Assuming the mean field relation Eq. (5), the observed phase boundary has a strong dependence on N t (see Fig. 3, top), which makes its interpretation questionable. This systematic error is dramatically reduced by using the nonperturbative prescription Eq. (28) for the renormalised anisotropy (see Fig. 3, bottom). Note that, under the nonperturbative prescription, the tricritical couplings on the temperature and chemical potential axes both decrease by ≈ 25%.
Moreover, analytic studies of the phase diagram generally consider Euclidean time as continuous [15], and should be compared with the N t = ∞ data only.  Figure 3. Phase diagram of SU(3) lattice QCD with massless staggered fermions, in the strong coupling limit, in which the anisotropy is set using mean field (top) [2], or using the present nonperturbative prescription (bottom). Under the nonperturbative prescription Eq. (28), the Nt-dependence of the phase boundary and of the tricritical point decreases substantially. Also, the tricritical couplings on the horizontal and vertical axes both decrease by ≈ 25%. The Nt = ∞ data is produced from simulations directly in the continuous time limit [14].

B. Static baryon mass
The static baryon mass am B is another observable for which the inexact calibration of anisotropy can have a strong effect. This observable can be determined using the "snake algorithm" [16], which samples partition functions Z k describing the system with an open baryonic segment of length k: We simulate massless SU(3) lattice QCD for different anisotropies using the snake algorithm, and estimate am B as a function of ξ (see Fig. 4). Under the two anisotropy prescriptions, Eqs. (5) and (28), baryon masses differ by ≈ 25% at large ξ. In this regime, the On an isotropic lattice, static baryons have mass am B ≈ 2.88 [1], and become heavier with anisotropy. In the continuous time limit, the baryon mass is only ≈ 20% heavier than the isotropic case, when using the nonperturbative prescription for the anisotropy, as compared with the ≈ 50% difference when using mean field.

C. Pion decay constant
Using our nonperturbative prescription for ξ, we can obtain reliable estimates of several physical quantities in the continuous time limit, e.g. the pion decay constant, aF π . In order to estimate this quantity, we measure the helicity modulus Eq. (35) for several finite hypercubic lattices and values of ξ. The results are summarised in Tables I and II, and displayed in Fig. 5 (top). The pion decay constant (squared) corresponds to the thermodynamic limit of the helicity modulus, in accordance with Eq. (34).
Again, the numerical data can be suitably fitted using the Ansatz Eq.

D. Chiral condensate
We also estimate accurate values for the infinitevolume chiral condensate a 3 Σ, by analysing the finitesize scaling of the chiral susceptibility a 6 χ, using chiral 4 New, direct measurements of aFπ in the continuous time limit [17] are consistent with our extrapolation. perturbation theory, and by using our nonperturbative prescription for the lattice anisotropy.
To this end, we estimate the chiral susceptibility density a 6 χ/N 4 s (as in [10]) for several finite hypercubic lattices and values of ξ (see Tables I and II). We estimate a 6 Σ 2 at finite ξ by extrapolating a 6 χ/N 4 s to the thermodynamic limit, modelling the finite-size corrections in accordance with chiral perturbation theory, see Eq. (36).
The dependence of a 3 Σ on ξ is again well described by the Ansatz Eq. (33) (see Fig. 5, bottom). The vertical intercepts give the values of the chiral condensate in the continuous time limit at T = 0: As before, U(3) and SU(3) are equivalent in the thermodynamic and continuous time limits, within errors. We also observe that, when keeping β 1 as a free parameter in Eq. (36), the finite-size fits are consistent with its theoretical value.

CONCLUSION
It is very important to have a precise scale for the lattice anisotropy. Even though mean field captures the correct power scaling of the renormalised anisotropy for asymptotically large values of the bare anisotropy, namely ξ ∼ γ 2 , it fails to predict the nonperturbative prefactor. The discrepancy between the mean field and nonperturbative prefactors introduces systematic errors of the same magnitude in many physical quantities of interest, particularly in the continuous time limit. This should be kept in mind when comparing strong-coupling Monte Carlo results and analytic mean field results, since the latter are usually formulated in continuous time.
In the dimer representation of the strong coupling limit of lattice QCD with massless staggered fermions, we have proposed a simple method to determine the nonperturbative dependence ξ(γ) between the bare and renormalised anisotropy couplings. The method is amenable to Monte Carlo simulations using very efficient directed path algorithms which, together with the multihistogram reweighting method, allows us to determine ξ(γ) with high precision. In the end, the nonperturbative prefactor is observed to be off by ≈ 25% with respect to the mean field prefactor.
As an application, we revisit the phase diagram of SU(3) lattice QCD [1], and update it using our nonperturbative relation ξ(γ). A strong dependence of the phase boundary on N t , introduced by the mean field anisotropy, essentially vanishes. The new locations of the phase boundary and of the tricritical point reveal corrections of ≈ 25%, in the chemical potential and temperature, compared with the old mean field values. We also compute the mass of the static baryon in the continuous time limit, which again receives corrections of ≈ 25% compared with the mean field value. These corrections are the direct consequence of the ≈ 25% correction to the mean field prefactor to ξ(γ) mentioned above.
We also estimate the values of the pion decay constant, aF π , and of the infinite-volume chiral condensate, a 3 Σ, in massless lattice QCD in the strong coupling limit at T = 0. The anisotropy corrections to these quantities are small, and provide a reliable extrapolation to their continuous time limits.
Even though the strong coupling limit of lattice QCD is unphysical, it may still be of interest to compare its predictions with those of continuum QCD, in the regime where chiral symmetry is spontaneously broken. For example, the strong-coupling SU(3) lattice value of the pion decay constant Eq. (40), in units of the critical temperature aT c ≈ 1.089 [11], is F π /T c ≈ 0.72, which is about 15% above the continuum QCD value.
Our approach can be generalized to the case of massive quarks. As before, in a hypercubic box the variances of the spatial and the temporal pion charges Eq. (19) can be required to be equal. Since they still scale as in Eqs. (23a) and (23b), the renormalisation criterion Eq. (22) is justified. What changes is that the pion charges are no longer conserved as per Eq. (18), i.e. have different values on parallel codim-1 lattice slices. A sensible observable is the average over such parallel slices of the variance of the pion charge on each slice. Thus, the setting of the anisotropy should be performed in a fixed volume L 4 , characterised by the value of m π L. This implies a fine-tuning of the quark mass, in order to keep fixed m π L while the bare anisotropy γ is varied. Alternatively, the anisotropy may also be set by keeping ξ and m q L = N t a t m q fixed while varying γ [17].
It may also be possible to extend the present study to finite β, in the framework of the O(β) partition function defined in [2]. The new occupation numbers (associated with plaquettes) introduce new Grassmann constraints on the extended configuration space. Such constraints may be used to construct analogues of the pion current, which would include plaquette corrections. In the chiral limit, we expect such currents to be conserved. The associated conserved charges could then be used to define nonperturbative renormalisation criteria for the (independent) spatial and temporal gauge couplings. An extension of this program to finite quark mass would be similar to the above proposal for β = 0.