Dual formulations of Polyakov loop lattice models

Dual representations are constructed for non-abelian lattice spin models with U(N) and SU(N) symmetry groups, for all N and in any dimension. These models are usually related to the effective models describing the interaction between Polyakov loops in the strong coupled QCD. The original spin degrees of freedom are explicitly integrated out and a dual theory appears to be a local theory for the dual integer-valued variables. The construction is performed for the partition function and for the most general correlation function. The latter include the two-point function corresponding to quark-anti-quark free energy and the N-point function related to the free energy of a baryon. We consider both pure gauge models and models with static fermion determinant for both the staggered and Wilson fermions with an arbitrary number of flavours. While the Boltzmann weights of such models are complex in the presence of non-zero chemical potential the dual Boltzmann weights appear to be strictly positive on admissible configurations. An essential part of this work with respect to previous studies is an extension of the dual representation to the case of 1) an arbitrary value of the temporal coupling constant in the Wilson action and 2) an arbitrary number of flavours of static quark determinants. The applications and extensions of the results are discussed in detail. In particular, we outline a possible approach to Monte-Carlo simulations of the dual theory, to the large N expansion and to the development of a tensor renormalization group.


Introduction
Dual representations of lattice gauge theories (LGTs) and classical spin models are a useful non-perturbative tool that allows to study many aspects of lattice quantum field theories. In the early days of LGT the dual transformations proved very efficient in the studies of the confinement and related problems, especially in the abelian gauge theories [1,2]. Also, dual representations appear to be very efficient for numerical simulations both at zero [3] and at finite temperatures for U (1) LGT [4]. The following years have seen many attempts to extend the duality transformations to non-abelian models using different approaches and strategies. In the pure gauge case the dual representation can be constructed starting from the plaquette formulation [5,6]. Dual variables are introduced as variables conjugate to local Bianchi identities [7,8]. The dual model appears to be non-local due to the presence of connectors in the Bianchi identities for gauge models. An analogue of the plaquette formulation for the principal chiral model is so-called link representation [9,10]. In this case one can construct a local dual theory for all U(N) and SU(N) principal chiral models [11]. Another approach is based on 1) the character expansion of the Boltzmann weight and 2) the integration over link variables using the Clebsch-Gordan expansion [12,13]. The resulting theory is the local dual theory written in terms of invariant 6j symbols. Several attempts to simulate this dual version have been undertaken (see Ref. [14] and references therein). In the opposite case, the strong coupling limit, the SU(N) LGT can be mapped onto monomer-dimer and closed baryon loop model [15].
During the last decade the dual representations have been applied to solving, fully or partially, the sign problem appearing in the lattice QCD in the presence of non-zero chemical potential and/or non-trivial topological term, like the θ-term. While it is still too early to say unambiguously if this approach can solve the sign problem in QCD, some advances in simpler models are encouraging. For example, the dual form of the massless two-dimensional U (1) LGT with one or two flavours of staggered fermions is free of the sign problem [16]. The same was proven in the strong coupling limit of the scalar QCD with one, two or three scalar flavours [17].
In general, there are two strategies attempting to construct positive Boltzmann weight for QCD or QCD-like theory at finite density. The first one relies on full integration over original degrees of freedom, i.e. gauge and fermion fields. The form of the final result strongly depends on the method of integration [18], [19], [20], [21], [22]. We do not discuss details of this approach here because in this paper we use a second strategy. It consists, first, in construction of an effective model for gauge loops winding around the lattice in the temporal direction, i.e. for Polyakov loops.
Only in the second step, the integration over Polyakov loops is accomplished. This strategy was successfully applied for the SU(3) Polyakov loop model in the strong coupling region for both temporal and spatial couplings of the Wilson action and in the heavy quark regime [23], [24], [25]. Similar results for U(N) models in the same approximations have been presented in [21]. More discussion on the effective Polyakov loop models can be found in Refs. [26], [27].
In this paper we calculate the dual representations for two Polyakov loop models. The Boltzmann weight of the first model is the same as the weight studied in [23]. We extend the results of [23] in several directions. First, our calculations are done for all values of N and in any dimension. Second, we consider the full static quark determinant with an arbitrary number of staggered or Wilson fermion flavours of different masses and chemical potentials. Finally, the result is given for the most general correlation function. These include, as particular cases, the partition function, two-point function related to the free energy of quark-anti-quark pair and for N-point function which gives the free energy of a baryon state. The Boltzmann weight of the second model is defined for all values of the temporal coupling constant, so that the strong coupling limit is imposed only with respect to the spatial coupling. Again, we treat all SU(N) models with an arbitrary number of static quark flavours and compute both the partition and correlation functions. We shall also explain how the results obtained can be easily transformed into results for U(N) and Z(N) Polyakov loop models. Boltzmann weights of all dual representations are non-negative, therefore our formulation could be used for Monte-Carlo simulations of the models at finite baryon or other chemical potentials.
This paper is organized as follows. In Sect. 2 we define the Polyakov loop models and introduce our notations. In Sect. 3 we derive dual representations for spin models in the strong coupling region of the temporal coupling constant. In Sect. 4 the result is extended to the arbitrary values of the temporal coupling. The possible applications and perspectives are discussed in Sect. 5. In the Appendix we explain all definitions and our notations related to the group representation theory. Also, we evaluate all group integrals appearing in the main text.

Polyakov loop models
We work on a d-dimensional hypercubic lattice Λ = L d with linear extension L and a unit lattice spacing. x ≡ x = (x 1 , ..., x d ), x i ∈ [0, L−1] denotes the site of the lattice, l = ( x, ν) is the lattice link in the ν-direction and p = ( x, µ < ν) is the plaquette in the (µ, ν)-plane. e ν is a unit vector in the direction ν. Periodic boundary conditions are imposed in all directions. Let G = U(N), SU(N); U(x) ∈ G, and dU denotes the (reduced) Haar measure on G. TrU will denote the fundamental character of G. The character of the irreducible representation λ will be denoted by s λ (U). The dimension of the representation is s λ (I).
In this paper we shall study some spin models on G with a local interaction in the external field and whose degrees of freedom are labelled by s λ (U). These models describe an effective interaction between Polyakov loops in (d + 1)-dimensional LGT with N f flavours of static quarks at finite temperature and non-zero quark chemical potential µ. The general form of the partition function of the models is given by On an anisotropic lattice and in the limit of vanishing spatial gauge coupling β s one can explicitly integrate out all spatial-like fields in any number of dimensions to get the following Boltzmann weight describing the Polyakov loop interaction (see, for instance Ref. [28] and references therein) The coefficients of this weight depend on the temporal gauge coupling β t ≡ β and can be expressed as Here I n (x) is the modified Bessel function and N t is the lattice size in the temporal direction. In the strong-coupling region β ≪ 1, the leading contribution comes from the fundamental character with coefficient D F (β), therefore the whole Boltzmann weight is approximated as The Boltzmann weight corresponding to the contribution of static staggered fermions can be presented, for N t even, as where the determinant is taken over group indices and The similar formula for static Wilson fermions reads In this case one has The unessential constants will be omitted in the following. When m f ≫ |µ f | or κ f ≪ e ±µ f one usually replaces these exact expressions with their approximation where h ± = s f h f ± , s = 1 for the staggered and s = 2 for the Wilson fermions. The Boltzmann weight of all these models is complex if µ f = 0 or, in general, if In what follows we assume h f ± in (6), (8) and (10) are arbitrary complexvalued variables. If h f ± are positive, the obtained dual weight is positive, too.

Dual of spin models I
In this section we consider the partition function (1) with the weight B g (β) given by (4). The static fermion contribution B q (m f , µ f ) will be taken either in its approximate form (10) or in exact forms (5), (7). The former case has been analyzed in [23] for SU(3) by making use of an exact parameterization of the SU(3) characters and measure. Consider the following Taylor expansion of the Boltzmann weight B g (β) exp TrU(x)TrU † (x + e ν ) For the fermion weight (10) we use the similar expansion To deal with exact static determinants (5) and (7) we use an expansion of the determinant in the Schur functions (Eq.(103) in Appendix), which is valid, in such generality, both for the staggered and for the Wilson fermions. Notations and some explanations regarding this formula are given in Appendix. We shall calculate the dual expression for the most general correlation function . (13) The partition function equals In what follows we analyze separately two cases: 1) heavy quark approximation (12) and 2) exact static determinant (103). All formulas below will be given for SU(N) models. In the end, we shall explain how one can easily obtain the corresponding dual representations for U(N) and Z(N) models using SU(N) results.

Heavy quarks
The original partition function in the presence of sources η(x),η(x) is given by Using (11) and (12) it can be written after some rearrangement as Here Q N (n, p) is a group integral defined and calculated in Appendix, Eqs.(83), (109) where j = min(n, p), d(λ) is the dimension of the representation λ of the symmetric group S n and the notation λ + |q| N is defined in the Appendix after Eq.(104). The integers n(x) and p(x) are given by where l i , i = 1, ..., 2d are 2d links attached to a site x and s(l

Pure gauge theory
Strictly speaking, the conventional duality transformations can be carried out only in the pure gauge theory, i.e. when h + = h − = 0 and, hence m(x) = k(x) = 0. Then, if j(x) = min(n(x), p(x)), the expression (15) takes the form while the constraint (20) reads This constraint can be solved in terms of dual variables in any dimension. It is important to emphasize that only Z(N) invariant correlation functions are nonvanishing due to above constraint. Indeed, taking into account that on the periodic lattice x d ν=1 (r ν (x) − r ν (x − e ν )) = 0 one can be assured that Eq. (23) implies that only invariant, i.e. mesonic and baryonic correlators of the Polyajov loops are non-vanishing in the absence of the external field (dynamical quarks).
In the following we consider, for the sake of simplicity, the two-point correlation function, corresponding to the free energy of the quark-anti-quark pair and the Npoint correlation function, corresponding to the N-quark (or baryon) potential. In the first case the sources are given by In the second case we introduce sources as We give below explicit formulas for d = 1, 2, 3 which follow from Eqs. (21) and (22).

One-dimensional model:
One-dimensional model is especially simple because we get from (20) where r ∈ [0, N − 1] becomes a global variable, k(l) ∈ [−∞, ∞] and η(l) = η for a set of links between sites x = 0, x = R and η(l) = 0 for links lying outside of the interval [0, R]. The delta-function in the 1st line of (21) is now δ k(l)−k(l−1),q(x) . Making a shift in q(x), the partition function with sources can be presented as where links l 1 , l 2 have a site x in common. Signs " + " and " − " corresponds to n(x) and p(x), correspondingly.

Two-dimensional model:
The solution of the constraint (20) in the two-dimensional model and in the presence of sources for the quark-anti-quark potential is given by the dual variables as (sites are placed in the center of original plaquettes, links are dual to links and sites become dual plaquettes) Here, η(l) = η if l ∈ S R , where S R is some path connecting points 0 and R, and η(l) = 0, otherwise. The partition function on the dual lattice takes the form where we have introduced notations Four links l i form a dual plaquette p with vertices and signs " ± " correspond to duals of n(x) and p(x) defined in Eqs. (17), (18). The solution of the constraint (20) in the presence of the baryon sources (25) can be constructed as follows. Let us take an arbitrary point x 0 and connect all N points x i with x 0 by some path S i consisting of dual links. Introduce dual variables as in (29), where η(l) = η if l ∈ S i and η(l) = 0, otherwise. The N-ality constraint (20) becomes where the plaquette p 0 is dual to the site x 0 and k(p) is the same as in (31). Strictly speaking, the solution of the form (29) is only valid in two dimensions if N ≤ 4. Then one can take all paths S i consisting of non-intersecting links and solution (29) holds. Though it is not a problem to extend the solution (29) to arbitrary N we restrict ourselves here to the case N ≤ 4. We thus conclude that the partition function in the presence of such baryon sources is of the form (30), where one has to substitute k(p) → k(p) + ηδ p,p 0 .

Three-dimensional model:
In the physically most relevant three dimensional case one obtains the solution of (20) in the following form Here, four links l i form a plaquette p dual to the original link l.
where S R is some path consisting of dual plaquettes and connecting points 0 and R, and η(p) = 0, otherwise. The partition function on the dual lattice reads c is a product over all cubes of the dual lattice and the notations are used Six plaquettes p i form a dual cube c and signs " ± " correspond to duals of n(x) and p(x) defined in Eqs. (17), (18). Extension of this result to the N-point correlation function is done precisely like in two-dimensional theory. In particular, if N ≤ 6 the solution of (20) can be taken as in (34). Then, defining paths S i , i = 1, · · · , N that connect points x i with some reference point x 0 (on the dual lattice path S i is formed out of plaquettes and connects cubes c i and c 0 which are dual to the corresponding sites) and introducing sources η(p) = η on plaquettes belonging to S i one finds that the N-point correlation function is described by Eq. (35) where one has to take the corresponding sources η(p) and make the substitution k(c) → k(c) + ηδ c,c 0 .
We can conclude that all three-dimensional SU(N) spin models are dual to the gauge models whose partition function is given by Eq. (35) with η(0) =η(R) = 0.

Full theory
Here we proceed with the full theory given by Eq. (15). Using N-ality constraint (20) one can sum up over k(x). With the help of notation we obtain after some manipulations the following expression We used here the property x r(x) = 0 and introduced parametrization The expression (39) is our final dual representation for SU(N) Polyakov loop models valid for all N and in any dimension. The functionQ N, ) and j ± (x) is given by Some comments are in order: • As follows from Eq.(39) and exact expression for the function Q N,q (j) given in (16), the dual Boltzmann weight is non- Hence, in this region the dual formulation can be used for the numerical simulations of the model with non-vanishing chemical potentials.
• Most thermodynamical functions and local physical observables, like the energy density, the baryon density, the quark condensate, etc. can be easily translated into the dual form by taking the corresponding derivatives with respect to β ef f , h f ± or µ f . This amounts to a local shift in a corresponding summation variable and can be presented as an expectation value calculated over the dual partition function.
• The long-distance observables, like two-and N-point correlation functions can also be written as an expectation values in the dual form. This follows directly from (39).

U(N) and Z(N) models
Here we explain briefly how the general result for SU(N) models can be used to compute the corresponding dual representations for U(N) and Z(N) models. The latter is equivalent to vector Potts models and can be obtained from SU(N) models by replacing U(x) matrices with their center elements. For simplicity we restrict ourselves here to the partition functions, i.e. η(x) =η(x) = 0.
U (N ) model: As explained in the Appendix, the only term contributing to U(N) group integrals is the term with q(x) = 0. Therefore, from Eq.(39) one gets for the partition function Z(N ) model: where I r (x) is the modified Bessel function, the partition function appears to be Let us add some more comments here: • Clearly, all comments given in the end of Sec.(3.1.2) remain valid for U(N) and Z(N) models.
• It follows from (42) that the partition function and invariant observables do not depend on the chemical potential for U(N) models with one fermion flavour.
In case of two flavours, the Boltzmann weight depends only on the difference of chemical potentials µ 1 − µ 2 .
• In the pure gauge case the corresponding representations for U(N) and Z(N) models can be straightforwardly obtained from Eqs. (27), (30) and (35).
• The dual form of the XY model can be calculated either taking N = 1 in U(N) case or as a limit N → ∞ in Z(N) case. E.g., in the pure gauge threedimensional case one recovers the following dual gauge-like form of the XY model

Exact static determinant
In this subsection, we compute the dual representation for the theory with the exact static determinant with an arbitrary number of flavours of the staggered, Eq.(5), or the Wilson, Eq,(7), fermions. As in the previous subsection, we shall calculate the dual expression for the most general correlation function defined in (13). The original partition function in the presence of sources η(x),η(x) is given by The gauge part of the Boltzmann weight B g (β) is treated as in the previous subsection using the expansion (11). Substituting this expansion into (47) one gets after some rearrangement Here, the function R N,N f (r, s; m f , µ f ) is a group integral defined in Eq.(85) of Appendix. The integers n(x) and p(x) are given by where l i , i = 1, ..., 2d are 2d links attached to a site x.
To deal simultaneously with staggered and Wilson fermions we use the representation (103) for the N f -flavour static determinant proven in Appendix. With this representation and making use Eq.(102) the group integral can be calculated exactly. This is done in Appendix, formulas (115)-(118). Presenting the N-ality constraint as we write down the final result (118) in the explicit form where c = 1 for the staggered and c = 2 for the Wilson fermions. The explicit form of the N-ality constraint is The variables H ± , depending on m f , µ f , and other notations are defined and described in Appendix. Combining last expressions with (48) the final result for the partition function with arbitrary sources gets the form We conclude this subsection with few comments: • All factors entering the Boltzmann weight of (54) are positive. This is true also for the Schur functions appearing in (52). Hence, this representation is suitable for the numerical simulations of the theory.
• Both long-distance and local observables can be written as expectation values over the dual partition function.
• An explicit form of the group integral for one flavour of the staggered fermions is presented in Appendix, Section (5). Even more detailed formula is given there for N = 3. This presentation can be directly used for the Monte-Carlo simulations.

Dual of spin models II
Here we investigate the partition function (1) with the weight B g (β) given by (2). The static fermion contribution B q (m f , µ f ) will be taken in its exact forms (5), (7). Like in the previous section we calculate the dual form both for the partition function and for the most general correlation function. In this case such correlations are given by an expectation value of the product of SU(N) characters taken in arbitrary representations η(x),η(x). Precisely, one has The partition function is recovered by taking trivial representations in all lattice sites. Exchanging order of the summations and integrations and rearranging product over links in the second line of (55) we write down the result in the form The coefficients D λ (β) are defined in (3) In the next subsections we specify this general formula for several important cases.

Pure gauge theory
In the pure gauge theory the group integral H d N,N f simplifies to G d N (λ i , γ i ) given by Eq.(84). The result of the integration is given by Eqs.(110)-(113). Denoting r(l) ≡ r ν (x) = |λ(l)|, the sum over all representations λ can be written as and the N-ality constraints (112) and (114) in the presence of sources read As before, we analyze each dimension separately. The result will be given for the partition function and for the correlation function of the general form. We shall also explain the particular solution of the N-ality constraints for the two-and N-point correlation functions. The corresponding sources can be taken like in Eqs. (24) and (25), respectively.

One-dimensional model:
One-dimensional model is exactly solvable. Using orthogonality relation (110) and expression for the Littlewood-Richardson coefficients C η λ 1 λ 2 (97) one finds for the partition function and for the two-point function Summations in the last expressions goes over SU(N) representations

Two-dimensional model:
In two-and three-dimensional cases the final result can be significantly simplified if we multiply the Schur functions in the integrand in a special way. Namely, in two dimensions we divide the lattice into a set of even and odd plaquettes and couple the Schur functions as shown in the left panel of Fig.(1). In this way the summations over original representations λ(l) are factorized inside every even plaquette. Letλ be a representation conjugate to λ, see Eqs.(89), (90). Extending the integration result (111) to the correlation functions and using decomposition (58), we obtain where the Boltzmann weight B p (ρ i (x)) and the function F (η(x),η(x)) are . . .
For the partition function F (0, 0) = δ ρ 2 (x),ρ 1 (x)+q(x) N and (62) reduces to The Boltzmann weight B p (ρ(x)) coincides with (63) up to replacement ρ 2 (x) → ρ 1 (x) + q(x) N . The N-ality constraint in (62), (64) has the same form as in (22). Therefore, both for quark-anti-quark sources and for baryon sources, one can use the solution (29). Because the link variables r(l) in (62) take on the non-negative values, the dual variable k(l) on the right-hand side of Eq. (29) is also non-negative and the difference r(x) − r(x + e n ) should be defined modulo N. The choice of η(l) for each case also remains as has been described after Eq. (29).

Three-dimensional model:
In three dimensions we divide the lattice in a set of even and odd cubes and couple the Schur functions as shown in the right panel of Fig.(1). In this way the summations over original representations λ(l) are factorized inside every even cube. Moreover, we first couple representations lying in the horizontal plane, then the resulting representations are coupled with representation sitting on the vertical links. The final step is to couple the representations obtained with representations η(x),η(x) where the Boltzmann weight B c (ρ i (x)) is The function F (η(x),η(x)) is of the form (63). The specification of this result for two-and N-point correlation functions is essentially the same as in the twodimensional case. The solution of N-ality constraint is taken as in (34) with restrictions described after (64). The resulting dual model is a three-dimensional model possessing local Z(N) invariance in analogy with (35).

Strong coupling limit
In the strong coupling limit, β = 0, the gauge part is absent and only static fermion contribution appears in the partition and correlation functions. Essentially, the model is nothing but the one-dimensional lattice QCD. The partition function is given by the integral (119) The result of the integration for the correlation function can be easily extracted from (120) These two formulas give an exact solution for one-dimensional QCD with arbitrary number of flavours of different masses and chemical potentials, both for the staggered and for the Wilson fermions. The detailed investigation of these solutions will be presented elsewhere.

Full theory: one-dimension
One-dimensional model corresponds to two-dimensional QCD. It is important to emphasize that our approach takes into account the full Wilson action in this case and the only though essential approximation is that we neglect the fermion propagation in one spatial direction. The corresponding one-site integral is given by Eq.(120). Substituting this into (56) one gets for the correlation function The Boltzmann weight and the function F (η(x),η(x)) are found to be The partition function is easily recovered putting η(x) =η(x) = 0.

Full theory: two-dimensions
To get a dual form for two-dimensional theory we use the result of the integration presented in Eqs.(120), (121) and follow the strategy used in the pure gauge theory. Namely, we first multiply characters from the gauge Boltzmann weight according to Fig.(1), resulting representations are coupled with representations from the fermion weight and, finally with representations η(x),η(x) from the correlation function.
After some algebraic manipulations we present the expression for the correlation function in the following form The function F (η(x),η(x)) is the same as in Eq.(72). The Boltzmann weight reads (74)

Full theory: three-dimensions
In three-dimensional theory we use the same strategy as above. The original link representations are coupled as shown in Fig.(1). With the help of Eqs.(120), (121) and using the same notations we obtain after long algebra the following representation for the correlation function The function F (η(x),η(x)) is the same as in Eq.(72). The Boltzmann weight reads Notation for sites of a cube are: As an important example, let us write down an explicit formula for the partition function with one flavour of the staggered fermions. In this case one has α = 1 k ,

In this case the Boltzmann weight (76) becomes
Let us add some remarks on resulting dual formulations: • All factors entering the dual Boltzmann weights obtained in this Section are positive including the Schur functions. Hence, this representation is suitable for the Monte-Carlo simulations at non-vanishing chemical potentials, at least in principle. Possible approaches to such simulations are discussed in the end of this section.
• One technical remark concerns the function F (η(x),η(x)) which appears in the presence of external sources for the correlation function. From its explicit expresssion, Eqs. (63) and (72), it follows that the product over lattice sites can be rearranged in a way that allows to include all Littlewood-Richardson coefficients from F (η(x),η(x)) into the Boltzmann weights. In this way the variables ρ i (x) become internal summation variables in each even plaquette (cube). Instead, variables σ(x) can be made dynamical variables of the dual theory.
• We have not presented explicit expressions for the thermodynamical quantities. They can be easily obtained by taking the corresponding derivatives with respect to β, h f ± or µ f . As in the region of the strong temporal coupling, this simply amounts to a local shift in a corresponding summation variable and can be presented as an expectation value calculated over the dual partition function.
• An explicit form of the group integral H d N,N f for one flavour of the staggered fermions is presented in Appendix, Section (5).
• The dual forms of U(N) and Z(N) models can be obtained following the lines described in the previous section. It follows from the N-ality constraint in Eq.(77) that the partition function and invariant observables do not depend on the chemical potential for U(N) models with one fermion flavour. In SU(N) and Z(N) models the dependence appears in the form e N qµ , as is expected on the general grounds.
Finally, we would like to discuss shortly some approaches to Monte-Carlo simulations of the dual formulations and related problems. Since the Boltzmann weights for the dual model obtained in this paper are nonnegative, they can be used for direct Monte-Carlo simulations. Dual weights based on the expansion of the group integrals into the Littlewood-Richardson coefficients might look quite complicated at first sight. Let us remark, however that if a positive dual weight at finite density exists, for the full theory it will certainly be much more complicated. Therefore, it is desirable to have a working algorithm for this complicated but still simplified dual theory. Here we give our thoughts on the way Monte-Carlo simulations can be performed. We will address explicitly the case of the partition functions (70), (73), (75) but the approaches described can also be applied to other cases.
The partition function as it is written includes summation over many sets of variables. While a direct Monte-Carlo simulation is possible, the convergence of the averages in this case can be slow. This can be overcome by dividing the variables into two groups -dynamical variables which are sampled in the Monte-Carlo simulation, and the variables over which the summation is done explicitly. For N = 3 such summation can be done by using the explicit values for the Littlewood-Richardson coefficients [29], [30]. If we take this approach, we can, for example, leave only |α(x)|, |β(x)| and ρ(x) variables as dynamical.
Another problem is that, while each configuration has nonnegative weight, many configurations formally allowed in the summation will have zero weight due to Littlewood-Richardson coefficients inside them becoming zero. Note that the Nality condition explicitly written in the partition functions is a necessary but not sufficient condition for the corresponding Littlewood-Richardson coefficients to be nonzero. This problem reduces the acceptance and convergence rate for any simple Metropolis-like update scheme and, on the more fundamental level, raises the question of ergodicity of the update process -one has to be sure that the whole acceptable configuration space can be probed.
For fixed small N values the full set of conditions for the Littlewood-Richardson coefficient to be nonzero can be written in the form of inequalities and one can either try to explicitly resolve them, or, at least, to build the update process respecting these inequalities. Another approach is to use the worm update algorithm [31], developed just for resolving such problems. In our case the worm has to propagate on an auxiliary lattice, that has the Littlewood-Richardson coefficients of the partition function as vertices, which are connected by a link if they share a common partition. Like for the Metropolis-like update, an explicit summation over part of partition variables can be done to reduce the phase space of the system -here it would amount to dividing the auxiliary lattice into blocks that are connected only by the links corresponding to the dynamical variables and treating each block as a site of a new lattice, while preserving the probabilities of the worm leaving the block through a given link.
If one wants to calculate the correlation functions, the set of acceptable configurations for the partition functions with sources and the one without sources will become different, requiring either rewriting the correlation function in terms of the valid configurations for the partition function without sources, or, if one uses the worm update, sampling the correlation functions using the worm algorithm in a way similar to the one described in [32].

Discussion
In this paper we have presented calculations of the dual representations for several Polyakov loop models. All these models have been derived in the strong coupling limit for the spatial coupling of the Wilson action. Contribution of the fermions is taken into account via the static determinant for an arbitrary number of the staggered and Wilson fermions of different masses and chemical potentials. Our results are valid in any spatial dimension and for all relevant groups, i.e. for SU(N), U(N) and Z(N). The main motivation is the construction of a positive Boltzmann weight in the presence of the baryon chemical potential suitable for numerical simulations. Some versions of representations from Sect. 3 have been derived before [23]. These formulations have been used for numerical computations of various local observables in [24], [25]. We have already applied our formulation (30) for studying two-and three-point correlation functions in two dimensional SU(3) spin models [33]. In [34] we have used the representation (54) to simulate the three-dimensional model with one flavour of the staggered fermions. In addition to local observables we have computed many two-point correlation functions in the presence of baryon chemical potential.
Let us briefly discuss other applications of the dual formulations.
1. One-dimensional model (70) with one flavour of the staggered fermions can be studied by the transfer-matrix method. Such study reveals the existence of an oscillating (or the so-called liquid) phase in some regions of the (h + , h − )plane [35]. This means the correlation function of the Polyakov loops while decaying exponentially is modulated by a periodic function. In other words, the mass spectrum of the theory becomes complex. The transfer matrix approach reveals the similar behaviour in one-dimensional Z(N) spin model in the external complex field [36], [37]. Monte-Carlo simulations of the same three-dimensional model also show the presence of such phase [37]. Detecting the liquid phase with the existing simulation methods at non-zero baryon chemical potential seems an extremely difficult problem. The formulations given in Sect. 3 might help to clarify if the oscillating phase exists in the threedimensional SU (3) LGT, at least in the region of validity of the Polyakov loop models used here.
2. It turned out that the dual formulations of Sect. 3 are well suited for the studies of the models in the large N limit. We have accomplished such studies and arrived at quite unexpected results: the large N 't Hooft limits of U(N) and SU(N) models are different in the presence of the chemical potentials. These results will be published elsewhere.
3. The partition function in Eq.(56) can be written at zero sources as where 2d links l 1 , . . . , l 2d are attached to a site x and the tensor T reads The rank of the tensor is 2dN. The trace can be done by properly contracting the indices labelling group representations. Such a formulation enables one to use the tensor renormalization group methods to study the theory at finite density and is certainly worth a separate investigation. (68) and (69)  Finishing this paper we would like to address the question of how one could systematically improve the dual formulations of Sect. 4? In the abelian case the dual construction can be extended to the full Wilson action. The underlying reason for this is that the exact and positive dual Boltzmann weight is known for all Z(N) and U(1) pure gauge LGT. Adding fermions in the form of the static determinant with any number of flavours does not destroy this property. Details of this formulation will be reported elsewhere. Much more difficult is the case of non-abelian models and the inclusion of the corrections to static determinant. In these cases one expects that the effective Polyakov loop model becomes non-local. For example, one such model describing non-local interaction between Polyakov loops has been derived in [38] (and Refs. therein) via the relative weight method

It would be interesting to investigate analytically solutions
The action of the model involves only fundamental characters, therefore the dual representations for this and similar models can be calculated by using the integration methods of Sect. 3. Clearly, if the kernel K(x − y) is positive for all distances considered the dual weight will be also positive but highly non-local. In general, the full effective action will contain all irreducible representations. More general effective action can be written in the form Even in this case the dual theory could be calculated with the help of integration methods of Sect. 4. The real challenge is to determine coefficients K λ,γ . One strategy is to expand the Wilson action at large spatial coupling and expand the fermion determinant around static contribution in powers of a lattice anisotropy. This will be the subject of future investigations.

A1: Definitions, notations and expansion formulas
First, we introduce some notations and definitions. Let λ = (λ 1 , λ 2 , · · · , λ N ) be a partition λ ⊢ r, λ 1 ≥ λ 2 ≥ · · · ≥ λ N ≥ 0 and l(λ) i=1 λ i ≡ |λ| = r, where l(λ) is the length of the partition λ. As a shorthand we will sometimes use a notation λ = a b to denote a partition consisting of b parts equal to a (i.e., λ i = a, 1 i b), and use λ + µ to signify elementwise addition of two partitions and λµ to signify a union of parts of two partitions. χ λ (σ) denotes a character of σ ∈ S r in representation λ. d(λ) = χ λ (1) is the dimension of the representation λ. The Schur function s λ (U) = s λ (u 1 , · · · , u N ) is a character of the unitary group G and u i -the eigenvalues of the matrix U ∈ G. s λ (I) is the dimension of the irreducible representation λ of G. One has The representation dual to λ will be denoted by λ ′ . The dual representation is defined by exchanging raws and columns in the corresponding Young diagram, i.e. λ ′ i = j 1 λ j ≥i . One has the following identity between the Schur function and its conjugate for U(N) group The similar identity for SU(N) group reads Given the complete symmetric functions h k and the elementary symmetric functions e k in m variables u 1 , · · · , u m h k = 1≤n 1 ≤···≤n k ≤m u n 1 · · · u n k , e k = 1≤n 1 <···<n k ≤m the Schur functions can be computed with the help of identities where λ ′ is a partition dual to λ and the following rule is understood s (λ 1 ,··· ,λn) (u 1 , · · · , u n−1 , 0) = 0, if λ n = 0 , s (λ 1 ,··· ,λ n−1 ) (u 1 , · · · , u n−1 ), if λ n = 0 .
Two other useful expressions for the Schur function read with τ being a permutation such that the number of cylcles of length j in τ is τ j .
The Weingarten function is used to evaluate the polynomial integrals over unitary groups. In the case of G = SU(N) this function is defined as where λ + q N = (λ 1 + q, · · · , λ N + q) and the sum in (95) is taken over all λ such that l(λ) ≤ N. For U(N) group one has to put q = 0. The Littlewood-Richardson coefficients C ν λ γ can be defined as coefficients appearing in the expansion of the product of two Schur functions From the orthogonality of the Schur functions one gets C ν λ γ are positive integers for unitary groups U(N) and SU(N) satisfying certain conditions. One such important condition on C ν λ γ to be non-zero arises from the integration over U(1) subgroup if G = U(N), while in SU(N) case it follows from the summation over Z(N) subgroup. Let U = Z be a center element of G. Then s λ (Z) = z r d(λ) , r = |λ| .
As follows from (97) the necessary condition for C ν λ γ to be non-vanishing is We refer to these conditions as to U(1) and N-ality constraints, respectively. More information on the Littlewood-Richardson coefficients as well as their closed forms for N = 2, 3, 4 can be found in [30]. Finally, we need the following formulas to treat the models with static fermion determinant. The first ones are the Cauchy identity and its dual where X = (x 1 , · · · , x L ), Y = (y 1 , · · · , y N ). The summation over λ runs over all partitions of NL such that l(λ) ≤ N and l(λ ′ ) ≤ L. The second one is an expansion of powers of the fundamental character into series over the Schur functions (TrU) r = (u 1 + u 2 + · · · + u N ) r = λ⊢r d(λ)s λ (U) .
With the help of Eq.(101) the fermion contribution given in (5) and in (7) for the staggered and the Wilson fermions, respectively, is presented in the form For the staggered fermions one has H ± = (h 1 ± , · · · , h N f ± ). The summation over α and β is taken over all partitions such that l(λ) ≤ N, l(β) ≤ N and l(λ . The summation over α and β is taken over all partitions such that l(λ) ≤ N, l(β) ≤ N and l(λ (6) and (8) for the staggered and the Wilson fermions, correspondingly.

A2: Group integrals
The Schur functions realize representations of U(N) group. Therefore, all integrals (83)-(86) are evaluated over the U(N) Haar measure. If G = SU(N) one should introduce an additional constraint into the measure This constraint can be implemented into the group integrals by multiplying the integrand with the delta-function ∞ q=−∞ (det U) q . Taking into account that one can easily prove with the help of (94) that If q < 0 one should replace the eigenvalues u i by u * i . Here and further we use the short-hand notation λ+q N = (λ 1 +q, · · · , λ N +q). The SU(N) constraint is enforced in the formulas below by summation over q. Furthermore, we shall present results only for the SU(N) group. The U(N) case is easily recovered by omitting all sums over q and taking q = 0 in all formulas below. More relevant information on the group integration and similar integrals can be found in Refs. [39], [40], [41].
Q N (r, s): To evaluate Q N (r, s) given by Eq.(83) we expand the traces in the integrand as sums over diagonal elements For the integral (83) this leads to The last integral can be calculated with the help of the Weingarten function. The details of the derivation can be found in [21,41]. Performing summation over group indices one gets Another way to compute (83) is to use the expansion (102). Then, the result (109) follows from the orthogonality of the Schur functions.
Three-dimensional case is treated similarly. One finds, e.g.

A3: One flavour of staggered fermions
For fixed values of N f many formulas given above can be specified and simplified by using explicit values for the Schur functions which, in turn, can be calculated from Eq.(92). As the simplest but important example, let us consider the integral R N,N f (r, s; m f , µ f ) with one flavour of the staggered fermions. Taking into account that in this case α = 1 k , β = 1 l , 0 ≤ k, l ≤ N and s α ′ (H + ) = h k + , s β ′ (H − ) = h l − , one gets from (118) the following simple answer Here the Littlewood-Richardson coefficients can be calculated using the following formula Finally, we specify the result (122) for the physically relevant case, namely N = 3.