Strong coupling methods in QCD thermodynamics

For a long time, strong coupling expansions have not been applied systematically in lattice QCD thermodynamics, in view of the succes of numerical Monte Carlo studies. The persistent sign problem at finite baryo-chemical potential, however, has motivated investigations using these methods, either by themselves or combined with numerical evaluations, as a route to finite density physics. This article reviews the strategies, by which a number of qualitative insights have been attained, notably the emergence of the hadron resonance gas or the identification of the onset transition to baryon matter in specific regions of the QCD parameter space. For the simpler case of Yang-Mills theory, the deconfinement transition can be determined quantitatively even in the scaling region, showing possible prospects for continuum physics.


Introduction
Quantum Chromodynamics (QCD) at finite temperatures and baryon densities is the theoretical foundation for an understanding of observed phenomena in a wide range of fields, like the thermal history of the early universe, heavyion collision experiments or the composition and properties of neutron stars. Since the interactions of QCD are strong on scales up to several GeV, ordinary weak coupling perturbation theory fails for temperatures and baryo-chemical potentials T, µ B < ∼ 2 GeV. The non-perturbative treatment by numerical Monte Carlo simulations of lattice QCD is to date still restricted to vanishing or small baryo-chemical potentials, because of a severe sign problem at finite baryon density (for an introduction, see [1]). The continued absence of a genuine algorithmic solution to this problem has motivated renewed interest in strong coupling and hopping expansion methods over the last decade. Contrary to weak coupling expansions, which typically result in asymptotic series, strong coupling expansions in the Euclidean framework are known to yield convergent series with well-defined radii of convergence. In the early days of lattice gauge theory they were used to get analytical results for some physical quantities of interest, such as glueball masses [2,3] or the energy density of lattice Yang-Mills theories [4,5]. These calculations were restricted to zero temperature, with the exception of some mean field analyses of phase transitions in the strong coupling limit [6,7,8,9], and a series for the temperature-dependent string tension [10]. More recent attempts to explore strong coupling methods as a systematic calculational tool for thermodynamics started with an investigation of Yang-Mills theory [11]. While Monte Carlo methods appear superior where they work, analytic results provide additional understanding for the interpretation of data. Moreover a two-stage approach, with an analytic derivation of effective lattice theories and their subsequent analytic or numerical solution, appears to be reasonably promising for finite baryon density. It is currently the only way to investigate the onset transition to baryon matter in lattice QCD directly, albeit not yet for physical parameter values. The purpose of this article is to review these developments with a focus on the concepts and the results, whereas technical details are left to the references.

A proof of concept: the SU (3) spin model
Before approaching the difficulty and complexity of QCD, it is useful to illustrate the basic strategies at the example of a simpler theory that is fully understood and has served several times as a testing ground for methods to deal with the sign problem. We consider the SU (3) spin model with the action S = − The fields L(x) = TrW (x) are complex scalars representing the trace of SU (3)matrices W (x). Writing η = κ exp(µ),η = η(−µ), the model closely resembles an effective action of lattice QCD with static quarks in the strong coupling region, to be discussed in Sec. 4. For µ = 0, it features spontaneous breaking of its Z(3) center symmetry along a line τ c (κ) of first-order transitions with an endpoint, as well as a sign problem at µ = 0. The model can be reformulated free of a sign problem in terms of a flux representation [12], allowing for simulations [13] by means of a worm algorithm. Likewise, a complex Langevin algorithm has been successfully applied [14,15,16].
Here, we are interested in a solution by series expansion methods, which are more familiar in condensed matter physics. Specifically, consider a linked cluster expansion of the free energy density (for an introduction, see [17]), which has been evaluated through 14 orders in τ , and for each order up to convergence in κ [18], The equation of state can be parametrised by the combined derivatives ∆S, corresponding to the interaction measure in QCD. It can be straightforwardly determined in Monte Carlo simulations, and all thermodynamic functions follow by integration or differentiation. Results corresponding to the three highest orders are shown in Fig. 1 (left). Excellent convergence and agreement with Monte Carlo results is observed until the phase transition, representing the radius of convergence, is approached. A finite polynomial is unable to describe a non-analytic phase transition, but is sensitive to it by loss of convergence. A marked improvement near the transition can be obtained by infiniteorder resummations through Padé approximants, defined as rational functions, The coefficients a i , b i are uniquely determined for L + M ≤ N , if N represents the highest available order of the expansion. In this way the [L, M ] approximant reproduces the known series up to and including O(x L+M ). Padé approximants are able to show singular behaviour and can model scaling properties near second-order phase transitions. The [6,6] approximant to ∆S is shown also in Fig. 1 (left). It accurately reproduces the simulated equation of state all the way to the phase transition, which it indicates by a singularity. At µ = 0 the phase transition is known to be first-order for small κ, ending in a critical point at some κ c , which is of the same type as standard liquid-gas transitions. Therefore the susceptibility and the specific heat, when approaching from any direction not parallel to the scaling axes, must both diverge with a common critical exponent, χ ∼ C ∼ |τ − τ c | − , = γ/βδ. This is reflected in a simple pole of the Dlog-Padé's, e.g., A difficulty is that Padé's also show poles for first-order transitions, where they indicate the end of the metastability range, as well as spurious poles, which renders the analysis quite intricate. In order to detect a true critical point one has to require convergence of the poles predicted by different Padé's, as well as convergence of he poles in both observables, D χ , D C . In this way, a critical point can be located within some scatter specifying a systematic error. Once the series are long enough to achieve this, an application to finite chemical potential, real or imaginary, poses no additional problem. Hence, the phase structure and the order of the phase tansition can be fully determined, as shown in Fig. 1 (right), in quantitative agreement with numerical results. This example provides further motivation to explore such methods also in QCD, where no algorithmic solutions to the sign problem exist yet.
3 Strong coupling and hopping expansions in QCD at finite T

Yang-Mills theory
Consider the Yang-Mills partition function with the standard Wilson action with plaquette variables U p and the lattice coupling β = 2Nc g 2 . From the lattice action it is obvious that the exponential can be expanded in powers of β, where the expansion point β = 0 corresponds to the strong (infinite) coupling limit. However, in the continuum limit β → ∞ and convergence is expected to be slow. A powerful resummation including all orders in β is immediately achieved by expanding instead in the group characters χ r (U p ) = TrD r (U p ), which are traces of the representation matrices D r of the plaquette variables. Using the formalism of moments and cumulants [5,19], it can be shown that this expansion exponentiates, such that a series for the (formal) free energy densityf ≡ − 1 Ω ln Z is obtained in terms of clusters C of graphs Φ, quantity in Eq. (2 . 2) is customarily called a free energy density, because the path integral corresponds to a partition function if one formally identifies β with 1/T . Here we are interested in a physical temperature T = 1/(aN t ), realised by compactifying the temporal extension of the lattice. The physical free energy is then obtained by subtracting the divergent vacuum free energy, and the pressure is P = −f . Group integrals are evaluated using the formulae Due to the latter the contributing graphs X i have to be objects with a closed surface.
The graph contributing to the leading order of the free energy density is a tube of length N t with a cross-section of one single plaquette, Fig. 1. Summing over all such graphs on the lattice, their contribution is Thus, the strong coupling limit at β = 0 (and thus T = 0) has zero free enery density or pressure, as expected. Moreover, we see that the pressure rises very slowly with T as the leading order starts at a high power in the expansion parameter only.
Corrections to the free energy density through order u 8 with respect to the leading order term for various N t have been calculated for SU(2) 1) and for SU(3). 5) Next, let us turn to screening masses. The one with the quantum numbers of the lowest lying glueball is defined by the spatial correlation function of plaquettes, .
(2 . 6) At zero temperature the exponential decay is the same as for correlations in the time direction, and thus determined by the glueball masses. At finite temperature, the LO graph for the difference to the vacuum mass is shown in Fig. 1 and gives The finite T effect is to lower the screening mass in the confined phase and is suppressed by high orders in the strong coupling. This explains why in the confined phase it is close to the vacuum glueball mass as observed by Monte Carlo. 6) Fig. 2 Leading order graphs for the free energy density (left) and screening mass (right).
Here Ω = V · N τ is the lattice volume, d r and a r (β) are the dimension and expansion coefficient of representation r, and c 0 is the expansion coefficient of the trivial representation. The combinatorial factor a(C) equals 1 for clusters C consisting of only one so-called polymer X i , which represent closed surfaces of plaquettes. The coefficients of higher representations can be expressed in terms of u(β) ≡ a f (β) = β/18+β 2 /216+. . ., the coefficient of the fundamental representation. It is a known function over the entire range of lattice gauge couplings with u(β) ∈ [0, 1), and constitutes the expansion parameter for the character expansion.
For thermodynamics the physical temperature T = 1/(aN τ ) is realised by compactifying the temporal extension of the lattice. The physical free energy is then obtained by subtracting the divergent vacuum contribution, and the pressure is P = −f . Group integrals are evaluated using the formulae the latter enforces contributing graphs X i to be objects with a closed surface. The leading order graph is a tube of length N τ , Fig. 2. Summing over all such graphs on the lattice, their contribution to the pressure is Thus, the strong coupling limit at β = 0 (and thus T = 0) has zero free enery density or pressure, as expected. Moreover, we recognize a qualitative feature known from full lattice simulations, where it is difficult to extract, and hadron resonance gas descriptions, namely the exponentially slow rise of the pressure with T . In order to get more quantitative, corrections to the free energy density through order u 8 with respect to the leading order term for various N τ have been calculated for SU(2) [11] and for SU(3) [20,21]. A direct comparison with Monte Carlo simulations is easiest for the energy density, which is shown in Fig. 3 for the example of SU (2). On the left, the series is shown order by order in the character coefficient u(β), and evaluated as a function of β. The series loses convergence as it approaches the deconfinement transition, whose critical coupling β c bounds the radius of convergence in this case. On the right, a resummation of the series is achieved by Padé approximants. This clearly improves the convergence and comparsion between different approximants provides an estimate for the systematic error associated with the truncated series. Good agreement with the Monte Carlo data is observed in the controlled region, which is however restricted to relatively small β-values.
At the finite temperature phase transition the free energy f (N τ , u(β c )) is continuous, with a discontinuous first or second derivative, depending on the order of the transition. For SU (2), we have a second-order transition, for which the 'heat capacity' diverges as C(u) ∼ (u c − u) α , with a critical exponent characteristic of the universality class. Its logarithmic derivative has a simple pole with residue α and is therefore well-suited for an analysis by Padé approximants. The averaged predictions of the approximants in Fig. 3 are β c = 1.65 (35) and α = 0.052 (19) [11], compared to the Monte Carlo value β c = 1.880(3) [22] and the 3d Ising exponent α = 0.12. Another interesting qualitative result can be obtained for thermal screening masses, which are obtained from the exponential decay of spatial correlation functions. The correlator of plaquettes is given by At zero temperature (infinite N τ ) the exponential decay is the same as for correlations in the time direction, and thus determined by the glueball masses.
At finite temperature (finite N τ ) the LO graph for the difference to the vacuum mass is shown in Fig. 2 (right) and gives Thus finite temperature lowers the screening mass in the confined phase, but the effect is suppressed by its high order in the expansion parameter. This explains the insensitiviy of screening masses to T observed in simulations of the confined phase, where they stay very close to the vacuum hadron masses up to the phase transition [23]. Such an understanding is not attainable from simulations alone.

QCD with heavy quarks
Including fermions for a treatment of full QCD requires an additional expansion. The Wilson fermion determinant can be formulated as a power series in the hopping parameter κ f = (2am f + 8) −1 for every flavour f . Defining the usual hopping matrix M in Dirac, colour and coordinate space, Grassmann integration over fermion fields gives to leading order in κ f per flavour The Kronecker deltas in M ensure that only closed loops contribute, and l counts their length. Again, for finite temperature effects we only need the difference between finite and infinite N τ , and hence are interested in the loops winding through the temporal boundary. One can formalise this by considering a hopping expansion in the spatial directions only, while the hops in the Euclidean time direction are taken into account fully. Thus the leading order heavy fermionic contributions to the effective action can be written as in Eq. (16), which only holds for finite T , and the dots represent loops that wind more than once. Now one can calculate again the free energy density and pressure, this time performing character expansions in both U and L. Expanding all terms up to O(κ 3Nτ ) and doing the group integrals one finds for two flavours u, d [21] Next we identify the hadron masses to leading order in the hopping expansion, Baryons : The pressure can then be rewritten as which is the expression corresponding to an ideal gas of hadrons. We have thus derived from first principles that QCD reduces to a hadron resonance gas in the strong coupling and heavy mass regime. Note that this also holds for Yang-Mills theory, which can be represented as a glueball gas [11,21]. It is then plausible that this feature is generic and independent of quark mass, which explains the success of hadron resonance gas descriptions of the lattice QCD equation of state in the confined phase, where the physical gauge coupling is fairly strong [24,25].

Effective lattice theories from strong coupling methods
In the previous section we have seen that strong coupling methods lead to interesting qualitative insights in QCD thermodynamics, but they are also limited quantitatively if we are interested in continuum physics. This can be improved significantly by a procedure in two stages: 1) construct an effective lattice theory by integrating out part of the degrees of freedom, and 2) solve the effective theory by series expansion or otherwise. Two types of effective degrees of freedom arise naturally, depending on the integration order, In the first case, fermions as well as all spatial link variables are integrated over, leaving a theory of temporal links only, which on a periodic lattice can be expressed by Polyakov loops. In the second case, all gauge links are integrated, leaving a fermionic effective theory in terms of mesons and baryons, because of gauge invariance. Both representations are perfectly equivalent to QCD. Truncations involved in doing the integrations reduce this equivalence to specific parameter regions, where the corresponding approximations hold.

Effective lattice theory for QCD with heavy quarks
The representation of QCD in terms of Polyakov loops was first developed to characterise the deconfinement phase transition in pure gauge theories [6,27], and has later been considered perturbatively and non-perturbatively in the continuum [28,29,30,31,32] and on the lattice [33,34,35,36,37]. Here our interest is in the derivation of the effective theories by strong coupling methods. In this case their form is uniquely determined and can be systematically worked out and improved. Starting point is Wilson's lattice formulation, as before. However, instead of the free energy we now compute the effective action defined in the first of Eq. (22), This is done again by combining a character expansion of the exponentiated gauge action with a hopping expansion of the exponentiated quark action. The gauge integrals can then be done term by term for the truncated expansion to leave an effective action depending on temporal links only. These combine to temporal Wilson lines W (x), which close through the periodic boundary and implicitly contain the dependence on the Euclidean time extent N τ . As a result, the effective partition function is three-dimensional, resembling a continuous spin model [38,39], Here the first line corresponds to pure gauge theory, the second line to the static quark determinant and the third line to the leading contributions of the kinetic quark determinant. The couplings of the effective theory are functions of the original lattice QCD parameters (for higher order expressions see [39]), and am = − ln(2κ) = am B /3 denotes the leading-order constituent quark mass in a baryon (not to be confused with the current quark mass m q in the QCD action). The higher the order to which the effective action is considered, the more couplings appear which are increasingly non-local, connecting loops W x to all powers and over all lattice separations. Note that the effective couplings are (resummed) power series in the original expansion parameters u and κ, and hence can in turn be treated as small expansion parameters themselves. Any truncated effective theory can thus be treated by series expansion methods as discussed in the simple example of Sec.2. On the other hand, the effective theory has a much reduced sign problem, and can be simulated by even a choice of different algorithms. One current line of work is to design flux representations for this type of actions that can cure the sign problem and be simulated efficiently [40,41].

Deconfinement transition in Yang-Mills theory
To appreciate the advantage of going through an effective theory, let us consider SU (3) Yang-Mills theory, so that h 1 =h 1 = h 2 = 0 and Eq. (24) is reduced to the first line. One can now compute the free energy and the susceptiblitiy of the Polyakov loop as a power series in λ, and improve the results by Padé approximants as in Sec. 3.1. These show a singularity at a critical λ c , indicating the phase transition. Note however, that for SU (3) we have a first-order phase transition, so that the singularity corresponds to the end of the metastability region of the disordered phase, i.e., the true critical coupling is slightly overestimated. With a series through O(λ 8 ) the estimates for λ c shown in Fig. 4 (left) are obtained, which agree with a Monte Carlo determination within 2%. This critical coupling can be converted into the critical Yang-Mills coupling β c (N τ ) by inverting the first of Eq. (26). Comparison with the full 4d Monte Carlo result [22] shows agreement to better than 10% for the entire range of N τ ∈ [2,16]. This is a dramatic improvement compared to the direct calculations in 4d QCD in Sec 3.1, which were only feasible for small N τ . This result also constitutes a completely analytic calculation of the deconfinement transition of lattice Yang-Mills theory. In fact, the β c -values for the larger N τ are in the scaling region and a continuum extrapolation of the critical temperature is possible, as shown in Fig. 5 (left). The resulting T c in the continuum also agrees within 10% with the known value from full 4d simulations. This calculation demonstrates that, in some cases, it is quite feasible to obtain continuum results from strong coupling methods.

Deconfinement transition in QCD with heavy quarks
Once quarks are included via the hopping expansion, a fully analytic evaluation becomes much more complex due to the additional couplings, but the effective theory can still be simulated efficiently. Moreover, while the partition function Eq. (24) still has a sign problem, it is considerably milder than that of 4d QCD, due to the fact that many degrees of freedom have been integrated out already. In [39] a flux representation of the effective theory including the static determinant was simulated, and the resulting phase diagram for N f = 2, N τ = 6 is shown in Fig. 5 (right). Since the center symmetry of the Yang-Mills theory is explicitly broken by the quark determinant, the first-order transition weakens with decreasing quark mass until it vanishes at some critical value m c q . The transition also weakens with chemical potential, so that for m q > m c q one finds a first-order deconfinement transition coming from the T -axis and terminating in a critical endpoint, whose location depends on the quark mass. A continuum extrapolation would require larger N τ -values, which in turn demand higher orders for the hopping expansion to converge. Note that recent 4d QCD simulations using a hopping expanded determinant [44], as well as full QCD simulations [45] permit control over the systematics at µ = 0.

Onset transition to baryon matter
With the same methods the cold and dense regime, where the sign problem is most severe, can also be studied. In particular, the onset transition to baryon matter has been seen explicitly to various orders in the combined expansions [46,43,42]. Once more valuable physical insight can be gained by an analytic calculation to leading order within the effective theory. In the strong coupling (β = 0) limit with a static quark determinant only, the partition function  Fig. 6 Left: Onset of baryon number (crossover) for heavy quarks and different lattice spacings [42]. Right: First-order onset transition for light quarks and low T (large Nτ ) [43].
factorises into one-site integrals which can be solved exactly. In the zero temperature limit, mesonic contributions are exponentially suppressed by chemical potential and for N f = 1 we have This corresponds to a free baryon gas with two species. With one quark flavour only, there are no nucleons and the prefactor of the first term indicates a spin 3/2 quadruplet of ∆'s, whereas the second term corresponds to a spin 0 six quark state or di-baryon. The quark number density is now easily evaluated Thus at T = 0 we find a discontinous transition when the quark chemical potential equals the constituent mass m. This reflects two expected physical phenomena: the "silver blaze" property of QCD, i.e. the fact that the baryon number stays zero for small µ even though the partition function explicitly depends on it [47]. Once the baryon chemical potential µ B = 3µ is large enough to make a baryon (m B = 3m in the static strong coupling limit), a discontinuous phase transition to a saturated baryon crystal takes place. The saturation density is 2N c quarks per flavour and lattice site and reflects the Pauli principle. This is clearly a discretisation effect that has to disappear in the continuum limit.
In the case of two flavours the corresponding expression for the free gas of baryons reads with now two h 1 couplings for the u-and d-quarks. In this case we easily identify in addition the spin 1/2 nucleons as well has many other baryonic multi-quark states with their correct spin degeneracy. A similar result is obtained for mesons if we instead consider an isospin chemical potential in the low temperature limit [43]. Remarkably, the entire spin-flavour-structure of the QCD bound states is obtained in this simple static strong coupling limit. When corrections are added the step function is smoothed out. Fig. 6 (left) shows the baryon density computed through orders O(κ 8 u 5 ) and plotted in continuum units. As the continuum is approached, the artefact of lattice saturation moves to infinity, as expected. Fig. 6 (left) shows a crossover, whereas for light quarks simulations show a first-order transition for large N τ , Fig. 6 (right), and a crossover for smaller N τ . This can again be understood with analytic insight from the hopping expansion. The dimensionless quantity here evaluated at β = 0 in the simpler case N f = 1 [43], gives an interaction energy per baryon in units of the baryon mass when T → 0. For µ < m, it evaluates to zero in accordance with the silver blaze property. For µ > m it is nonzero and implies an attractive interaction, consistent with the baryon gas 'condensing' to a liquid (or crystal on the lattice). Starting as ∼ κ 2 , it decreases with growing quark mass to zero in the static limit, as one would also expect from Yukawa potentials in nuclear physics. Hence, the end point of the liquid gas transition, T c (m), shrinks to zero with increasing quark mass, which explains the difference between the left and right of Fig. 6.

Large N c and quarkyonic matter
Interesting conjectures concerning the QCD phase structure were made some years ago, based on large N c arguments [48]. In particular, the phase diagram in the limit of large N c , with fixed 't Hooft coupling, was argued to look as in Fig. 7  plasma phase, where the pressure scales as p ∼ N 2 c (perturbation theory), from the hadron gas phase, where it scales as p ∼ N 0 c (ideal gas of hadrons). In [48] it is argued that at finite density there should then be a third phase with p ∼ N c , which was termed quarkyonic since it shows aspects of both baryon and quark matter. In particular, at low temperatures the fermi sphere in momentum space is argued to be composed of a baryonic shell of thickness ∼ Λ QCD , and quark matter inside.
With the analytic expansion tools and an effective lattice theory at hand, one can now adress these issues by direct calculations starting from QCD. For large N c , the baryon mass grows as m B ∼ N c , so the size of the constituent quark mass should not matter in that limit. This suggests to investigate the cold and dense region for large N c by direct calculation with the effective theory for heavy QCD. First, the effective theory of the previous sections was derived for a general number of colours [50] (see also [51]). Second, the onset transition to baryon matter has been analysed for general N c in [49], where no large N c approximations are implied. Fig. 8 (left) shows the baryon onset transition to steepen with growing N c , such that it always become first-order in the large N c limit when the t'Hooft coupling and N τ are held fixed. This is independent of the value of N τ , which parametrises temperature. This strenghtening of the onset transition implies its critical endpoint T c (N c ) to move to ever larger temperatures with growing N c , as indicated in Fig. 7 (right). Together with the straightening of the deconfinement line with N c , which follows from the behaviour of Feynman diagrams at large N c , we thus understand how the conjectured phase diagram in Fig. 7 (left) emerges smoothly from the phase diagram with N c = 3 and heavy quarks. Furthermore, right after the onset transition h 1 > ∼ 1, through three orders in the hopping expansion at β = 0, the coefficients of the pressure are proportional to ∼ N c once N c is large, which is suggestive of a property to all orders. This N c -scaling is stable under gauge corrections and, with a leading order correction in N c , holds all the way down to N c = 3 − 9, Fig. 8 (right). Note also, that a lattice filling with baryon number smoothly changes from baryon matter (at the onset of condensation) to quark matter (at saturation) as a function of µ B , which is consistent with the momenum space picture of quarkyonic matter.
It thus seems that lattice QCD with heavy quarks shows the defining features ascribed to quarkyonic matter. On the other hand, for the moderate densities right after the baryon onset, one would expect no discernible difference from ordinary baryon matter. It remains an open question what happens at larger densities, in particular for the case of light quarks, when there may or may not be an additional chiral transition after baryon onset.

Effective lattice theory for chiral and light quarks
The light quark regime has also been addressed already a long time ago by effective theories at strong coupling. In particular, it was attempted to describe dense baryonic matter in Hamiltonian approaches, with constructions based on the strong coupling limit and partial chiral symmetry on the lattice [52,53,54]. Here, we stick to the Euclidean approach, starting from the partition function, to apply systematic expansions. However, being interested in the chiral phase transition, we now consder the lattice action with staggered fermions, x+µ χ x,l , where ( , x, p) label lattice links, sites and plaquettes and η µ (x) are the staggered phases. In the continuum limit, this represents N f = 4 QCD if no rooting is applied. Expanding in β, one obtains with plaquette (anti-plaquette) occupation numbers {n p ,n p }. In the strong coupling limit, β = 0, the gauge action vanishes and the link integration reduces to one-link integrals [55]. Subsequent integration over the Grassmann variables then yields the partition function as a sum over graphs in terms of mesons and baryons [56]. Leading O(β)-corrections have been derived in [57].
Recently, after an intricate reorganisation of the involved degrees of freedom, an all-order dual formulation in terms of monomers, dimers, world lines and world sheets has been given [58,59]. Its contributions through O(β 2 ) can be simplified to The integer valued occupation variables satisfy constraints at sites and bonds, which result from the Grassman and gauge integrations. Note that this formulation in general also contains negative weights, but the remaining sign problem is mostly mild enough to be handled by reweighting techniques. A particular advantage of this formulation is the feasibility to simulate the chiral limit directly, as well as any finite quark mass. On the other hand, gauge corrections are more difficult to include. Moreover, in the present formulation the expansion of the gauge sector in β is slower to converge than the previously discussed character expansion, to which it might be extended in the future. Nevertheless, this formulation offers yet another road to studying QCD at finite density in a complementary parameter regime compared to the previous sections. Early mean field [61] and Monte Carlo [62] studies based on a polymer representation have been restricted to the strong coupling limit, β = 0. Modern simulations are done using a worm algorithm [63], which can be extended to include gauge corrections. Fig. 9 (left) shows the phase diagram in the chiral limit, m = 0, both for β = 0 and with leading gauge corrections O(β) included [57]. In the chiral limit there always is a non-analytic chiral phase transition, with a tricritical point where the first-order transition at finite density meets the second-order line. In the strong coupling limit, this tricritical point coincides with the end point of the nuclear liquid gas transition.
When gauge corrections are switched on, these start splitting up, but surprisingly the first-order lines of the chiral and nuclear transitions are still indistinguishably close. Fig. 9 (right) shows the strong coupling limit, but now with finite quark mass switched on. The second-order transition line changes to crossover, as expected. A very intersesting feature is the decrease with mass of the endpoint T c (m) of the nuclear liquid gas transition. Obtained from a different discretisation and in a very different parameter range, this behaviour agrees with the finding for heavy quarks discussed in the last sections. Note also that the endpoint of the chiral transition quickly moves to µ B > ∼ 3T , which is consistent with all results reported for QCD at the physical point [64,65].

Conclusions
In view of the sign problem of QCD at finite baryon density, which continues to evade purely algorithmic solutions, systematic studies of strong coupling methods, either directly or through development of effective theories, have provided a number of genuinely new results for QCD thermodynamics: -The behaviour of screening masses for temperatures up to the crossover -A derivation of the hadron resonance gas in the strong coupling regime -The deconfinement transition for heavy quarks at all baryon chemical potentials -The chiral phase transition in the massless limit in the region of strong, but finite, couplings -The onset transition to baryon matter and its equation of state, both for heavy quarks on fine lattices and for light quarks on coarse lattices -A smooth connection to thermodynamics at large N c These results were so far obtained in parameter regions away from the physical point in the continuum. Nevertheless, all the physics expected up to nuclear densities has been seen qualitatively on the lattice. It is now a matter of pushing these approaches to higher orders, and to evaluate them with refined algorithms. Given existing automatisation techniques to achieve high-order calculations for spin models and their sign problems, known in the condesed matter literature, the author believes this path is worth to be pursued further.