Equation of state for cold and dense heavy QCD

A previously derived three-dimensional effective lattice theory describing the thermodynamics of QCD with heavy quarks in the cold and dense region is extended through order $\sim u^5\kappa^8$ in the combined character and hopping expansion of the original four-dimensional Wilson action. The systematics of the effective theory is investigated to determine its range of validity in parameter space. We demonstrate the severe cut-off effects due to lattice saturation, which afflict any lattice results at finite baryon density independent of the sign problem or the quality of effective theories, and which have to be removed by continuum extrapolation. We then show how the effective theory can be solved analytically by means of a linked cluster expansion, which is completely unaffected by the sign problem, in quantitative agreement with numerical simulations. As an application, we compute the cold nuclear equation of state of heavy QCD. Our continuum extrapolated result is consistent with a polytropic equation of state for non-relativistic fermions.


Introduction
A fully non-perturbative description of QCD at finite baryon density, and in particular of nuclear matter at zero temperature, remains an outstanding challenge because of the infamous "sign problem" of lattice QCD, which prohibits direct Monte Carlo simulations of that regime. Approximate methods like reweighting, Taylor expansion about quark chemical potential µ = 0 or analytic continuation of results at imaginary chemical potential, for which there is no sign problem, all require µ/T < ∼ 1 in order to work [1]. So far, no signal of a critical point or a first order phase transition has been found in this controlled region, where the different approaches are in quantitative agreement. Simulations by complex Langevin algorithms do not suffer from the sign problem and a lot of progress in extending their range of operation has been achieved [2,3]. Nevertheless, coarser lattices still appear inaccessible and so far mostly the heavy dense limit of QCD has been studied [4]. Up to now, no non-analytic phase transitions have been reported from this approach.
It is therefore desirable to further develop alternative approaches where the sign problem is fully controlled, even if those are restricted to certain parameter regions of QCD. Chiral fermions reformulated to a flux representation can be simulated in the strong coupling regime by means of a worm algorithm, and gauge corrections can be successively included [5]. Also analytic expansion methods are attempted in that regime [6]. Here we pursue the complementary approach and consider QCD thermodynamics with heavy quarks but much closer to the continuum. This situation is described by a 3d effective theory derived by a combined character and hopping expansion and features a mild sign problem only, allowing for the simulation of real chemical potentials and the mapping of both the hot and cold regions of the phase diagram [7][8][9]. This approach has also been successfully tested against two-colour-QCD [10], where there is no sign problem.
The phase diagram of QCD with heavy quarks is sketched in Fig. 1. At zero density there is a first order deconfinement transition related to the spontaneous breaking of the centre symmetry in pure gauge theory. Real chemical potential for quarks weakens that transition, which then features a critical end point. The thermal phase transition and the dependence of the critical end point on the quark or pion mass has been calculated in [7] on the lattice. The same picture emerges in studies of continuum Polyakov loop models [11] or Dyson-Schwinger equations with heavy quarks [12]. At low temperatures and higher chemical potentials, also the nuclear liquid gas transition is accessible [8]. The temperature of its critical endpoint is of the order of the nuclear binding energy and also depends on the quark mass. In the infinite mass limit it moves to zero temperature [9].
In the present work, we extend the results of [9] in two ways. First, we push the derivation of the effective action for the cold and dense regime through order u 5 κ 8 to leading order in N −1 τ . Second and most importantly, we apply linked cluster expansion methods [13] to our effective theory and demonstrate that its thermodynamic functions and equation of state can be computed entirely analytically in the domain of its validity. We then devise a resummation scheme to sum up a particular class of diagrams and finally compute the equation of state for heavy nuclear matter.

Derivation
The derivation of the effective theory has been discussed in previous publications [7,9,14], so we only outline the procedure and give our results. Starting point is lattice QCD with the Wilson plaquette and fermion actions on an N 3 s × N τ lattice, e aµ f δ ν0 (1 + γ ν ) αβ U ab ν (x)δ x,y−ν + e −aµ f δ ν0 (1 − γ ν ) αβ U ab −ν (x)δ x,y+ν , and U ab −ν (x) = U †ab ν (x −ν). We denote colour indices with Latin characters and Dirac indices with Greek characters throughout the article. The effective action is defined by integrating over the spatial link variables, and we have split the effective action into contributions coming from the pure gauge theory and the fermion determinant. We now specify N f = 2 degenerate quarks with κ u = κ d = κ. In our approach we approximate the exponential in (2.4) by truncated expansions in the fundamental character of the gauge group, u(β) = β/18 + O(β 2 ) [24], and the hopping parameter κ = (2am + 8) −1 . After the expansion, the gauge integration can be done analytically. The resulting effective theory is three-dimensional and only depends on temporal Wilson lines or, equivalently, Polyakov loops, with the one-point coupling constants to leading order The static determinant has a particle-hole symmetry about half-filling akin to the Hubbard model [15]. The expansion is then in spatial hops of the remaining kinetic determinant, For our physics region of interest, the cold and dense regime µ T , considerable simplifications arise. At fixed lattice spacing the zero temperature limit corresponds to N τ → ∞. The centre-symmetric couplings, λ i , have been calculated in previous publications and tested against the full Yang-Mills theory [14,16,17]. They are suppressed as λ i ∼ u nNτ with n ≥ 1 and u(β) < 1 always. In this work we employ β ≤ 6.2 and N τ ≥ 116, such that λ 1 < ∼ 10 −18 and λ i ≤ λ 1 . Thus the pure gauge sector plays no role in the cold and dense regime and can be safely neglected. Similarly,h 1 → 0 in the zero temperature limit. The summation of all temporal windings produces the basic building blocks of the effective action, We have calculated the effective action through order κ 8 u 5 in the low temperature limit, i.e. to the leading power of N τ . Because of its length we will give the result here in a compact, graphical representation and relegate the full expression to appendix A. We symbolise factors of W n,m ( x) by vertices, where n is the number of bonds entering a vertex, and m is the number indicated on the node. Furthermore, vertices which are connected by one or more bonds are nearest neighbours on the lattice.
The sums over the "degrees of freedom" constitute the traces in coordinate space. The effective couplings to the order computed here are Note that higher order corrections in κ, u to h 1 and h 2 are subleading in the high N τ limit.

Observables
We are interested in the thermodynamical functions, which are directly related to the partition function, in particular the baryon number density, pressure and energy density, 14) Here a is the lattice spacing and z = exp(aN τ µ) the fugacity.

Hadron masses and the physical scale
In order to interpret the results in the following sections, we use the masses of mesons and baryons for N f = 2 to all orders in the hopping expansion [18] and to resummed next-to-leading order in the strong coupling expansion, with M = 1 2κ . Explicit evaluation shows that these formulae are remarkably convergent in the heavy mass regime for β < ∼ 6.2, which we employ in our analysis. Since heavy quarks have little influence on the running of the coupling we use the beta-function of pure gauge theory for the lattice spacing in units of the Sommer parameter, a(β)/r 0 with r 0 = 0.5 fm [19]. Temperature is then set via T = (aN τ ) −1 .

Convergence of the effective action
We simulate our effective theory as described in [9] by cross checking complex Langevin simulations with simulations using standard Metropolis updates and reweighting. Our first task is to assess the range of validity of our new action. One expects the additional orders in κ to extend the convergence region, within which the description of thermodynamic functions by the effective action is reliable. We test this by computing the baryon number density at fixed values of the coupling h 1 and N τ . Varying κ then allows us to assess the convergence of the expansion of the kinetic quark determinant. Fig. 2 (left) shows the results obtained with effective actions of increasing order in κ. One observes clearly how two adjacent orders stay together for larger values of h 2 (κ) as the order is increased, thus extending the range where our effective action is reliable. Fig. 2 (right) shows the same exercise for the largest κ considered in this work, this time increasing the orders of the character expansion. We observe good convergence up to β ∼ 6, which is a sufficiently weak coupling to allow for continuum extrapolations. It is interesting to note that the convergence properties are not determined by the size of the expansion parameters alone. Even though the u(β)-values far exceed the κ-values employed in the figures, convergence in u(β) appears to be faster. The gain in convergence region by the additional orders in the effective action can be exploited to study the systematics of our effective theory.

Continuum approach
An important question for any lattice investigation concerns the continuum limit. Fig. 3 (left) shows the baryon number as a function of chemical potential and highlights a severe issue of lattice QCD at finite baryon density, irrespective of the sign problem or the accuracy of effective actions: cut-off effects at finite density cause not only quantitative systematic errors, but alter the qualitative behaviour of the system. Because of the finite number of lattice sites available, the Pauli principle leads to a saturation density of n sat B = 2N f baryons per site, which does not exist in the continuum. Once lattice saturation is reached, a further increase of chemical potential makes no sense. Thus lattices have to be made finer before higher densities can be addressed. On finer lattices the saturation density in physical units grows and in the continuum limit moves to infinity. This lattice artefact starts to make itself felt already quite early, as is also apparent in the numerical behaviour of the Polyakov loop [8] and related to the half-filling symmetry of the static action [15].
The difficulty is also reflected in Fig. 3 (right), where the slopes of the continuum approach rapidly increase with growing chemical potential, such that a continuum extrapolation is increasingly difficult to control. The figure shows results from our previous simulations obtained with the κ 4 action at two values of µ > µ c , i.e. beyond the nuclear onset transition, and compares it with the new κ 8 action. The baryon density just about reaches the domain with leading cut-off effects linear in a, which are expected for standard Wilson fermions. In this context it should prove particularly valuable to work with an im-proved action with O(a) lattice corrections removed. For still finer lattices the data points break away from this behaviour, signalling the limit of validity of our finite series. We conclude that the hopping expansion is systematic and controlled, with additional orders in the action allowing to access finer lattices, but the progress is very slow. For sufficiently heavy masses and not too high densities, we have thus attempted a continuum extrapolation [9], which we reproduce in Fig. 7 below.
Another comment regarding our extrapolation is in order. The region with linear cut-off dependence starts at a ∼ 0.1 fm, resulting in am π ∼ 10. Thus our lattices are too coarse to resolve the structure of hadrons, which effectively appear as point particles, and one might wonder how this could possibly be consistent with continuum physics. Indeed, the hadronic mass values corresponding to the formulae (2.17, 2.18) are afflicted by large cut-off effects and do not represent the true mass values in the continuum. In principle this could be repaired by the methods of heavy quark effective theory [20], which however is beyond our present interest. On the other hand, the nuclear physics in this parameter region is effectively governed by the interactions between baryons and not within baryons. Moreover, in the case of very heavy mesons the Yukawa potential between nuclei is extremely short ranged, i.e. in that limit the nucleons really do interact as point-like particles as in our setup. Thus our extrapolation should reflect continuum physics, though we only roughly know the hadron masses this limit corresponds to.

Mass dependence
A second way to benefit from the additional orders in the hopping expansion is to keep the lattice spacing fixed and study smaller masses. This is shown in Fig. 4 for two different lattice spacings. The error bars in these plots are systematic and give the difference between results obtained by the action to the highest two orders in the hopping expansion. Growing error bars thus indicate the loss of good convergence and control. Again, this behaviour is in complete accord with qualitative expectations, with increasing orders in the hopping expansion making smaller quark masses accessible and coarser lattices allowing for lighter quarks. However, these results also illustrate the fundamental difficulties and limitations of an effective theory based on the hopping expansion. While the systematics appears to be controllable and reliably tell us about its breakdown, the gain in mass range per additional order in the hopping expansion appears to be too small to envisage an extension to the physical quark masses of QCD at present.

Linked cluster expansion for the effective theory
So far we have derived an effective, three-dimensional theory for QCD thermodynamics in the cold heavy mass regime and used it for numerical simulations. However, we have observed in previous publications [9,16] that the couplings of the effective theory are sufficiently small to suggest a perturbative calculation of thermodynamic functions. In this section we develop a systematic expansion scheme by applying the linked cluster expansion, well-known from spin models [13]. In order to apply it, we change variables to Polyakov loops L( x), which are complex numbers and resemble continuous spins. For the transformation of the measure, see [14]. The rational expressions of Wilson lines in the effective action can be converted using the generating functional given in appendix B.

General framework
We begin by summarising the basic features of the linked cluster expansion, for a more thorough review, see [13]. Consider an N -component scalar field with a 2-point coupling, which may also extend over larger distances than nearest neighbour, All information on the interaction is encoded in v ij (x, y), which we assume to be small. We will see later that in our case v ∼ κ 2 . Our goal is to study thermodynamic quantities, so we are interested in the free energy rather than the partition function, The linked cluster expansion is thus defined by the series expansion of W in powers of the coupling, A systematic way of taking the derivatives with respect to the coupling is by introducing source terms to define the generating functionals .
The derivatives of the free energy with respect to the sources are the cumulants, or connected n-point functions, e.g.
Finally, setting the interaction to zero means that the cumulants only give a contribution for fields on the same site, which we use to define the n-leg expressions M i 1 ...in , . . .
Thus we get the series expansion of W,

Graphs and embeddings
The last expression suggests a graphical notation where the M 's are n-legged nodes and the v's are bonds connecting them. It is apparent that the order of a node is determined by the number of bonds entering it, e.g., (4.10) With this notation we can express the expansion (4.9) by graphs, The prefactors give the symmetry factor of a graph, which is the inverse number of ways one can label the bonds and the nodes while keeping the same mathematical expression (i.e. connecting the same pairs). The expansion of W to some power n of v now requires computing all graphs with n bonds. So far the graphical notation does not contain information about the spatial dependence of a graph. For a translationally invariant theory all spatial dependence can be summed up in a single embedding number, which counts the number of ways to put a graph on the lattice. It depends on the type and lattice distance of the interaction, as well as the dimension and geometry of the lattice we are working on. For example the v 3 term cannot be put on a square lattice, and thus its embedding number is 0, while on a triangular lattice it would be non-zero. A quick summary of the lowest order graphs with symmetry factors and embeddings is given in table 1.

Graph Symmetries Embeddings Graph Symmetries Embeddings
We are now ready to map our effective theory to O(κ 2 ) onto this computational scheme. The partition function to this order contains a nearest neighbour interaction between two W 1,1 -terms, We can thus apply the results from the linked cluster expansion for a one component field The building blocks for (4.9) are then

Higher order couplings
At O(κ 4 ) we are confronted with 3-point couplings. Fortunately, introducing higher npoint interactions to the linked cluster expansion is straightforward. In our case we need a generalised partition function which has a cluster expansion where the derivative with respect toũ is once more given by the cumulants, . (4.17) The geometry of the interaction term is contained in u ijk (x, y, z). For example if we take φ as a two-component field, φ = W 1,1 , W 2,1 , the three-point O(κ 4 ) term has an interaction tensor corresponding to a wedge. In this case the linked cluster expansion of W is the sum of all diagrams which can be made out of these two components, where the three new diagrams come from the 2-point and 3-point terms in the κ 4 -action.
Note that directions are necessary to distinguish a node W 2,1 from W 2 1,1 . This also changes the symmetry factor. It is thus possible to compute all graphs from combining elements up to a certain order, carefully calculating symmetry factors as one proceeds to higher and higher orders.
Alternatively, and as an independent check, one can use the idea of embedding graphs from the effective action onto the basic graph topologies of the cluster expansion. As an example, consider the square graph symmetry: 8 . (4.20) The following O(v 4 ) terms can be embedded on it, with an embedding number, that counts the number of ways this is possible.
Graph Embeddings Symmetry This will modify the symmetry factor to be # of unique embeddings graph symmetry factor . Thus we can systematically get the linked cluster expansion for our full effective theory, by writing down all topologically distinct diagrams and then embed our effective theory terms onto the graphs. The result of this endeavour is too lengthy to include in this publication, but the interested reader may obtain the full result by contacting the authors.

Results
We are now in a position to evaluate thermodynamical functions completely analytically, presently we have computed through order κ 8 . Fig. 5 shows the evaluation of the κ 2action (left) and the κ 8 -action (right) in various orders of the linked cluster expansion in comparison with the numerical evaluation, again in the strong coupling limit β = 0. It is interesting to note the faster convergence for the lower order action, as one might expect. For the higher order action, the linked cluster expansion first has to "catch up" to the order of the action before it can start resumming its contributions. Comparing with Fig. 2 we observe that, to the orders computed, the linked cluster expansion converges roughly as far as the derivation of the effective theory. If we allow for 10% deviation between different orders and from the full result, the combined calculation is valid up to h 2 ∼ 0.08.
Having coefficients order by order in κ 2 , we next attempt to improve the convergence of the series by constructing Padé approximants. These have been used previously to effectively resum the strong coupling and hopping expansion of the QCD deconfinement transition [21] and considerably improve convergence compared to the straight expansion [22]. Padé approximants are rational functions constructed from power series of order N and are defined as In order to uniquely determine the coefficients a i , b i , it is necessary to have L + M ≤ N , if N represents the highest available order of the original expansion. In this way an [L, M ] approximant is correct up to, but not including O((κ 2 ) L+M +1 ), and larger approximants represent more expansion coefficients than smaller ones. Because of the choice of L, M for longer series, there is no unique Padé approximant for a given series. We discard approximants that produce unphysical singularities through a vanishing denominator. Of the remaining ones, the diagonal or close to diagonal ones are expected to be most reliable [23]. Fig. 6 (left) shows a considerably improved convergence when using Padé approximants for consecutive orders of the linked cluster expansion. An important quantity characterising nuclear matter is the binding energy per nucleon. It can be defined thermodynamically by the energy density minus mass density in the zero temperature limit, In previous work we have shown numerically as well as analytically to leading order in the hopping expansion, that this quantity displays the silver blaze property, i.e. it is zero until the onset transition, where it becomes negative [9]. We can now extend this study to slightly larger densities. Fig. 6 shows the binding energy extracted from the Padé approximants to the partition function at various orders. While quantitative convergence breaks down    shortly after the onset transition near 3µ ∼ m B , we obtain a new qualitative result: in higher orders we see the binding energy becoming positive again with growing chemical potential, as is expected from nuclear physics. A minimum characterising bulk nuclear density appears, which however is not yet settled quantitatively at the available orders.
Finally, we evaluate our effective theory by linked cluster expansion in a regime where we are able to fully control it, which we again monitor by systematic errors taken as the difference between consecutive orders. This is the heavy and cold regime which we have studied numerically in [9]. We evaluate the equation of state for 6 different lattice spacings in the range of 0.079 fm < a < 0.136 fm and perform a continuum extrapolation. The resulting baryon density as a function of chemical potential is shown in Fig. 7 and quantitatively agrees with the numerical results.

A chain resummation
So far we have managed to reproduce most simulations of the effective theory with analytic calculations. In this section we will present a resummation scheme for the analytic approach which extends beyond the reach of the numerical evaluation, in that it generates and then sums additional terms in the effective action.

General idea
We start with a motivating example. Consider the following four terms in the effective action Eq. (2.11), It is easy to see that these four terms follow a common pattern to generate a chain. Each term extends the length of the chain by one node while maintaining a common prefactor.
Looking at the equations we see that every link in the chain adds a factor h 2 W 2,1 to the term, along with the necessary spatial geometry. One can check with the terms up to order κ 8 in Eq. (2.11) that this holds not only for the simple one-string-chain shown above, but for all terms in the action with a singly connected node, meaning nodes that correspond to a factor of W 1,1 . The "chain" resummation can schematically be represented as   where n is the total number of links attached in the chains, and in the end we sum over all n. In the formulae this amounts to the substitution where G({x n }) contains the geometry of the chain. Although the pattern is simple, to show that the prefactors come out in a way that is summable is quite involved. A more thorough analysis of the resummation is given in appendix C. We now have a resummation of a class of diagrams to all orders in κ for every order in the effective theory, but evaluating the final gauge integral is impossible. This is because we need to sum over all geometries for all the terms in the resummation, which cannot be evaluated analytically. To proceed, an additional constraint must be introduced, namely that only embeddings with the same basic geometry as the starting structure are included. This implies that all nodes of the chain will be at separate lattice points, and an n + 1 long chain will result in the following integral (cf. appendix C), ≡ (2dh 2 ) n+1 I n 2,1 I 1,1 . The I n,m denote the integrals over W on a lattice point occupied by a single W n,m factor. To carry out the integral one needs again to convert from W to L using the generating functional presented in appendix B. The full chain will then give

Validity of the combinatorics
Carrying out the chain resummation and embedding it on the simplest cluster expansion graph as outlined above does introduce small systematic errors. This is because the embedding factor of (2d) also counts graphs belonging to self-overlapping chains and thus is too large. However, the difference of the overcounted contributions and the actual selfoverlapping embeddings, which represents the error, results in cumulants of the participating factors. In Fig. 8 examples of these cumulants from the overlapping embeddings of a specific term have been plotted. We see that the non-overlapping configuration is orders of magnitude larger than the overlapping ones, especially as the fugacity h 1 approaches 1.
The same holds for the other types of graphs as well. This behaviour is due to the fact that the integrals dL W n 1,1 and dL W 1,1 n are of the same magnitude, resulting in cancellations in the cumulants. Therefore, the error introduced by our embedding is small.

Results
With a new and extended effective action we can redo the calculations from section 4.4. In Fig. 9 (left) the increase in convergence due to the resummation scheme is clearly visible, more than doubling the convergence region in h 2 . Matching the plot with that of Fig. 6,   one sees a comparable increase in convergence to that from the Padé approximation. This is both expected and reassuring as both approaches produce rational expressions, and the superior convergence of the Padé is expected due to the fact that it is not restricted to a particular class of diagrams and might therefore predict higher order behaviour. In Fig. 9 (right) we have repeated the pion mass convergence plot and one can see that the resummation extends the convergence region in a natural way.
We now give our final result, the equation of state for nuclear matter with heavy quarks calculated fully analytically, Fig. 10. The error bars represent the uncertainty resulting from continuum extrapolations including a varying number of points. The line represents a fit to a polytropic equation of state for non-relativistic fermions, This might of course be expected on physical grounds: the bosonic baryons condense and do not contribute to the pressure, which is thus due to the heavy fermionic baryons. However, computationally this is a very remarkable result. Firstly, we have started with an action in terms of quarks and gluons and expanded about the strong coupling limit. The emergence of baryonic degrees of freedom with a weak attractive interaction [9] is completely dynamical and a result of the calculation. Secondly, the equation of state for any finite lattice spacing shows saturation and thus lattice fermions do not feature a polytropic equation of state, which is also clearly visible in Fig. 10. The fact that our continuum extrapolation is well described by a physically sensible polytrope then appears to be an endorsement of our calculation. It will now be very interesting to investigate the prefactor of the polytropic behaviour, which must depend on the contributing degrees of freedom, their masses and interactions.

Conclusions
We have extended a previously derived three-dimensional effective lattice action for QCD to the order u 5 κ 8 in the cold limit with a combined character and hopping parameter expansion, starting from the full Wilson action. The effective action has a sign problem mild enough to permit controlled simulations of the cold and dense regime for heavy quarks, where our expansion is valid. The additional orders fully confirm our previous results and demonstrate that the systematics in this approach can be monitored and controlled.
In the second part of the paper we have exploited the fact that our effective theory formally corresponds to a spin model with multi-point couplings over various ranges, which we were able to map onto a linked cluster expansion. This permits a completely analytic evaluation of cold and dense thermodynamics which is entirely unaffected by the sign problem. The convergence of the linked cluster expansion is excellent and fully reproduces the numerical simulations of the effective theory in the range of its validity. In this framework we were furthermore able to identify a class of diagrams consisting of chains of arbitrary length representing meson exchange, which can be summed up to all orders in the hopping expansion.
For sufficiently heavy quarks, continuum extrapolations of the thermodynamical functions are possible to provide the equation of state for heavy bulk nuclear matter. Our final result is consistent with a polytropic system of non-relativistic fermions as expected on physical grounds.

Appendix A Effective action in the cold limit
The effective action to order κ 8 and to leading order in 1 Nτ reads where the sums of {i, j, k, l} go over all spatial directions. At this order in 1 Nτ all gauge corrections come from the rescaling of the coupling constants h 1 (u, κ) and h 2 (u, κ). Moving away from this limit will result in gauge corrections that depend on the spatial geometry of the various terms.
Appendix B Generating functional for the W n,m terms In our calculation of the effective theory we need to convert terms on the form: to be functions of the Polyakov loop L = tr W . We can accomplish this using the following generating functional, G[α, β] = tr ln α + βh 1 W , together with the trace-log identity and the expression which is derived from the Calay-Hamilton theorem. It is then straightforward to show that Appendix C The chain resummation C.1 Higher order terms in the effective theory As a prerequisite to our chain resummation, we need to understand the systematics and coefficients of the terms in the effective action that we wish to sum up. To this end, we recall the expansion of the kinetic determinant, (2.9), in terms of forward and backward spatial hopping matrices [9] P αβ ab (x, y) = κ T αδ ac (x, z) i∈{x,ŷ,ẑ} Here T is the static propagator whose Dirac and colour structure neatly separate, T αβ ab (x, y) = δ x,y δ a,b δ α,β + 1 − γ 0 αβ B ab (τ x , τ y )δ x, y , (C.4) and the expression for B can be found in [9]. After the integration over the spatial gauge links, only a few combinations give non-zero contributions to the final result. This is because gauge integrals over many combinations of link variables vanish [25], dU U n U † m = 0, if n + 2m = 0 (mod 3). We define contractions to be those terms for which the gauge integration gives non-vanishing contributions. Thus a contraction identifies the matrix products with such spatial and temporal coordinates that the gauge links overlap, e.g., In the low temperature limit, where N τ → ∞, the contractions consisting of pairs will dominate the result and we restrict our attention to those for the remainder. At this stage we no longer need to distinguish between hops in positive or negative spatial directions. In contrast to temporal hops, which get boosted by the baryon chemical potential, there is no asymmetry between them, and the gauge integration only depends on the number of links in a term. We thus switch to a notation focussing on the dominant pairings, where X and Y symbolise the remainder of the term. Every contracted P, M pair is labelled by an arbitrary number i, and the terms are therefore invariant under relabelling. The six pairings of tr P M P M P M in table 2 are contained in tr 1 1 2 2 3 3 , tr 1 2 3 3 2 1 , tr 1 2 3 1 2 3 . (C.10) The number of equivalent labellings can be read from the notation, as the three terms have 2, 3 and 1 distinct cyclic permutation(s), respectively. In this notation the non-zero terms The combinatorial prefactors 1/g of the trace products are determined by the symmetries of the individual terms. We have 1 g = # of unique cyclic permutations of the traces n 2 !n 4 ! · · · n N ! 2 n 2 4 n 4 · · · N n N .
The numerator is the number of cyclic permutations within all traces that stay different under relabelling. The n i in the denominator is the number of trace factors over i matrices (e.g. the third term in (C.11) has n 2 = 2 and the fifth n 4 = 1), and N is the total number of matrix factors (or the order of κ) of the term. · · · 1 1 · · · W 1,1 ( x), (C.13a) · · · 1 2 2 1 · · · W 2,1 ( x)W 1,1 ( x + i), (C.13b) · · · 1 2 3 3 2 1 · · · W 2,1 ( x)W 2,1 ( x + i)W 1,1 ( x + i + j). (C.13c) The prefactors of terms in this resummation can be calculated from symmetry arguments. Assume that we know the symmetry prefactor 1/g of a term with N open ends. Extending one of these can be done in N/g distinct ways, which all break the previous symmetry. The sum of all such insertions thus have a prefactor of N/g which is N times that of the base diagram. Instead of counting the number of permutations we can think of the number of ways to add n links to N open ends. The total combinatorial factor for a graph with N open ends, an internal symmetry of g and with n link insertions is therefore (C.14)

C.4 Recursive gauge integration for the chain
Finally we will compute the spatial link integrals for the chain and see that we reproduce the substitution in Eq. (5.3). We argued that chain of length n can be represented as an assumption we will finally settle in this section. The expression for the chain has a recursive structure, and it is therefore natural to define the matrices G m such that where the sum over the temporal variables could be trivially evaluated as the W n,m 's are independent of their time argument.
To tie it all together let us consider a generic contribution which has N open ends, where the matrices M i are the rest of the term, comparable to the left hand side of Eq. (5.2). Inserting the expression for G 1 gives We can now attach chains of length n i to each of the N open ends so that the total length of the chain is n, corresponding to the right hand side of Eq. (5.2) C n = n 1 ,n 2 ,...,n N 1 g {n i } tr G n 1 +1 (τ 1 , τ 2 ; x 1 )M 1 G n 2 +1 (τ 3 , τ 4 ; This gives a symmetry factor that depends on the partitioning of the attachments {n i }. We insert the expression for G n from Eq. (C.27), which gives us where x jk is the k'th position of the chain originating from the j'th open end, corresponding to x j from Eq. (C.28). The degrees of freedom are the directions the hops can take. We see that the base diagram, what is explicitly left in the trace, is the same for the term with and without attachments. There is the integral over temporal gauge links which we evaluate by embedding this term onto a skeleton cluster expansion graph with the same geometry. We sum over the subclass of non-overlapping chains, meaning that the coordinates x jk never overlap with each other, nor the positions of the base diagram. This subclass of terms is labelled C * . The integrals over the remaining temporal gauge links therefore factorise and we can sum all partitions of the attachments {n i } into one term. The sum of the symmetry factors of the partitions is exactly what we computed in Eq. (C.14), n 1 ,n 2 ,...,n N 1 g {n i } = 1 g N − 1 + n n , (C.32) which means that the integral over this embedding, C * , is Here Ω is the embedding factor of the base diagram and every factor in the chain brings a lattice embedding of 2d as argued in section 4. Finally we can sum over the total length of the attachments, n, and for brevity we only include the n-dependent factors of the previous expression,