Heavy dense QCD and nuclear matter from an effective lattice theory

A three-dimensional effective lattice theory of Polyakov loops is derived from QCD by expansions in the fundamental character of the gauge action, u, and the hopping parameter, \kappa, whose action is correct to \kappa^n u^m with n+m=4. At finite baryon density, the effective theory has a sign problem which meets all criteria to be simulated by complex Langevin as well as by Monte Carlo on small volumes. The theory is valid for the thermodynamics of heavy quarks, where its predictions agree with simulations of full QCD at zero and imaginary chemical potential. In its region of convergence, it is moreover amenable to perturbative calculations in the small effective couplings. In this work we study the challenging cold and dense regime. We find unambiguous evidence for the nuclear liquid gas transition once the baryon chemical potential approaches the baryon mass, and calculate the nuclear equation of state. In particular, we find a negative binding energy per nucleon causing the condensation, whose absolute value decreases exponentially as mesons get heavier. For decreasing meson mass, we observe a first order liquid gas transition with an endpoint at some finite temperature, as well as gap between the onset of isospin and baryon condensation.


Introduction
The phase diagram of QCD at finite temperature and baryon density is still largely unknown today, because lattice QCD suffers from a severe sign problem when chemical potential for baryon number is non-vanishing. Several methods have been devised to circumvent this obstacle (see e.g. [1] and references therein), but all of them introduce additional approximations that are valid for small quark chemical potentials only, µ/T 1. In order to reach higher chemical potentials and/or low temperatures, methods are required which at least potentially may solve this problem. Among these are Complex Langevin Dynamics (CLD) [2,3], transformation of the degrees of freedom into so-called dual variables as demonstrated in scalar models [4,5], and the formulation of quantum field theories on a Lefshetz thimble [6]. In particular, CLD has recently been applied to full QCD in a previously inaccessible parameter range [7]. However, even if an approach should finally succeed in solving the sign problem, it will remain very hard to study the regime of cold and dense matter. This is because, in order to avoid the limiting artifact of saturation at finite lattice spacing, very fine lattices are required for high density, which implies in turn very large temporal lattice extents near T = 0.
In this work we further elaborate on an effective theory approach [8][9][10][11], where analytic strong coupling and hopping expansion methods are used to derive an effective lattice action whose numerical simulation is feasible also in the cold and dense regime. The sign problem can be handled by complex Langevin simulations in a controlled way, and in certain parameter ranges even Monte Carlo simulations are possible. Moreover, the effective action resembles a three-dimensional spin model, such that the numerical effort is vastly smaller than for full lattice QCD simulations. At the present stage of the project, simulations can still be run on time scales of days on university PC clusters. The drawback is that the effective action is only valid in parameter ranges where the expansion converges, which is currently restricted to the heavy mass region and the confined phase. Even there, the effective theory is unsuitable for long range correlation functions, but it gives accurate results for bulk thermodynamic quantities and phase transitions [12]. In particular, it has already provided predictions with better than 10% accuracy for the critical couplings of SU(2), SU(3) Yang-Mills [8], the critical quark masses where the deconfinement transition changes to a crossover [9] and the tricritical point of the deconfinement transition at imaginary chemical potential [13]. A similar approach is used in [14][15][16][17] with staggered fermions. There, the chiral regime can be studied directly but the strong coupling series is much harder to compute and no continuum extrapolations are possible so far.
The paper is organised as follows. In section 2 we summarise the derivation of the effective action in the pure gauge sector and give a systematic calculation of the fermion determinant. In section 3 we analyse the effective action by analytic methods to leading and next-to-leading order in the small effective couplings. Section 4 is devoted to a systematic study of the validity of complex Langevin simulations. Finally, section 5 contains our physics results for the cold and dense regime of QCD with heavy quarks. Readers not interested in the technical aspects of the derivation and simulation may skip sections 2, 4 and read sections 3, 5 only.

The effective action
Starting point is a (3+1)-dimensional lattice with Wilson's gauge and fermion actions for N f flavours, which after Grassmann integration may be written as with elementary plaquettes U p , the quark hopping matrix for the flavour f , and U ab −ν (x) = U †ab ν (x −ν). The effective action is then defined by integrating out the spatial link variables In the first line we split into a part which is Z(N c ) centre symmetric and a part with symmetry breaking terms. For the present work it is more convenient to split the action into a purely gluonic part and a fermionic part due to the determinant, which contains both symmetric and symmetry breaking contributions. All terms depend only on temporal Wilson lines W x or their traces, the Polyakov loops, The effective action features an infinite tower of interaction terms between loops to all powers and at all distances, where S x i denote i-point-interactions. These are completely determined in terms of Wilson lines and the parameters of the original theory. Note that, without truncations, the effective action is unique and exact. Nonperturbative determinations of the effective action [18][19][20][21][22] can in principle be applied at all parameter values. In practice they necessarily imply truncation and modelling, which may have to be different in different parameter regimes. In our approach we compute the effective action in a combined strong coupling and hopping parameter expansion, which orders terms according to their leading powers in β, κ. By summing up all temporal windings we make sure that we have the complete dependence on chemical potential, or fugacity, in each order of the hopping parameter expansion.

Pure gauge theory
For the Yang-Mills part, it is advantageous to perform a character expansion where the factor c 0 (β) can be neglected as it is independent of gauge links and cancels in expectation values. In earlier publications [8,9,23], we have shown how to compute the effective gauge theory up to rather high orders in the fundamental character expansion coefficient u(β) ≡ a f (β) = β 18 + . . . . In leading order we have a chain of N τ fundamental plaquettes winding around the temporal direction and closing via periodic boundary conditions. Therefore the leading order is a two-point interaction, where higher order corrections of λ(u, N τ ) as well as a discussion of higher order interaction terms can be found in [8]. In the leading order expression of eq. (2.6) we already see that λ(u, N τ ) is suppressed for large N τ , since u < 1, see also [9] for a further discussion. In this work we aim at temperatures T ≤ 10 MeV with lattice parameters β < ∼ 6, N τ ≥ 100, where λ < ∼ 10 −25 can be safely neglected.

Static quark determinant
The quark determinant is expanded in a hopping expansion. In order to keep the complete dependence on chemical potential it is useful to split the quark matrix in positive and negative temporal and spatial hops, The static determinant is then given by neglecting the spatial parts, with propagation in the temporal direction only. Calculating the space and spin determinant we get (2.9) Note that this includes all windings of Wilson lines around the temporal direction and thus the full fugacity dependence. A well-known relation valid for SU(3) then allows us to reformulate this in terms of Polyakov loops, with the abbreviation 11) and the constituent quark mass am = − ln(2κ) = am B 3 , to leading order of eq. (2.49). When det[Q stat ] is exponentiated, the parameter c also constitutes the effective onepoint coupling constant of S f 1 to leading order [9], h 1 = c,h 1 =c . (2.12)

Kinetic quark determinant
In order to compute a systematic hopping expansion about the static limit, we define the kinetic quark determinant which we then split into parts describing quarks moving in positive and negative spatial directions, P = k P k and M = k M k . The reason for this is that the trace occurring in eq. (2.13) is also a trace in coordinate space. This means that only closed loops contribute and hence we need the same number of P s and Ms in the expansion of the logarithm.

Static quark propagator
We now compute the static quark propagator (1−T ) −1 appearing in eq. (2.13). Since (1 + γ µ )(1 − γ µ ) = 0, hops in forward and backward time direction do not mix and the full static quark propagator is given by In order to compute the positive static quark propagator, we use the series expansion The inverse is then given by is a temporal Wilson line from τ 1 to τ 2 and we have suppressed its spatial location. If τ 1 = τ 2 , the Wilson line winds around the lattice, W (τ 1 , τ 1 ) = W . The contribution in negative time direction (Q − stat ) −1 τ 1 τ 2 can then be obtained from (Q + stat ) −1 τ 1 τ 2 by the following replacements and reads Finally we split the temporal quark propagator in spin space as well as in propagation in positive and negative temporal direction according to

Gauge integrals for the leading fermionic action
Next we compute the leading strong coupling contribution to the effective action by performing the group integrations. We will arrange the fermionic part of the effective action as Since it is not known how to analytically perform the gauge integral over links in the exponent, we have expanded it in a Taylor series. After the integration we shall see that it is possible to resum some terms back into an exponential. At the order κ 4 there are zero-point contributions (or vacuum graphs) from closed hops around a plaquette. In a strong coupling series these only contribute after being dressed with a plaquette, ∼ κ 4 u, and thus are neglected here. The one-point contributions of the Polyakov loops constitute the static determinant and have been computed already.

Two-point interaction
Dealing with more than one trace, as in i Tr P i M i 2 , it will be necessary to explicitly display spatial coordinates, i.e. (2.26) Here we have to consider three different possibilities: The two nearest-neighbour contributions may share 0, 1 or 2 sites, where only the last one contributes to the two-point interaction. To the order κ 4 it is then The first contribution to the two-point interaction is of order κ 2 : Here we have used the expressions eq. (2.24) for B, evaluated the trace over spin and coordinate space and introduced the coupling The group integrations have been computed via Note that this enforces the spatial link variables to be at the same temporal location and yields a factor N τ rather than N 2 τ from the two temporal traces. From now on we will skip the last step, where one has to insert the definitions of A and B and perform the temporal sums. Explicit expressions for all types of terms appearing in the following can be found in the appendix.
The next correction to the two-point interaction is of order κ 4 : In this calculation it can happen that there is a spatial link which is occupied by four matrices and we need the group integral (see e.g. [24]) The next contribution of order κ 4 comes from eq. (2.26), which is a two-point interaction in the case that x = y and i = j: Higher corrections to the two-point interaction start with O(κ 6 ).

Three-point interaction
The three-point interaction starts at O(κ 4 ); The different contributions are evaluated to be where the sum is only over terms where the two traces share one spatial point.

Four-point interaction
There are only two four-point interactions to order κ 4 : After integration the first contribution vanishes in the strong coupling limit and only gives a non-zero contribution if a plaquette is inserted into the fermionic loop: Since we only calculate the action to order κ m u n with m + n = 4 we neglect this term. The second contribution is where the sum is only over terms where the traces share no common spatial point.

Resummations
In order to include as many terms as possible and improve convergence we perform a resummation. Note that in order to perform the gauge integration, we had to expand the exponential of hopping matrices, e.g., After the integration it is possible to resum many of the resulting terms back into an exponential, Nc u,i Tr Bu,uTr B u+î,u+î + O(κ 6 ) . (2.44) Inspection of higher order terms indicates that this should always be possible. Note that terms that have been resummed, like the higher orders in eq. (2.44), have to be excluded in the appropriate higher order to avoid double counting.

Leading gauge corrections to the strong coupling limit
Leaving the strong coupling limit, i.e. for β = 0, the gauge action has to be included when performing the group integration. As a consequence the effective coupling constants depend on the gauge coupling also. The leading gauge corrections are of order N τ κ 2 u coming from attaching plaquettes to the Wilson line, cf. figure 1, and This can also be exponentiated by summing over multiple attached, disconnected plaquettes at different locations and we see that in this way the Polyakov line receives mass corrections due to interactions. Note that this generates overcounting in higher orders, but in our opinion the resummation effects of this procedure more than compensates for this additional care. Let us finally also give the gauge correction for the prefactor of the leading order of S f This correction does not appear to exponentiate.

Effective action for the cold and dense regime
The terms evaluated in the last sections and displayed in the appendix can now be added up to provide the complete effective action. Fortunately, simplifications occur because some terms have the same structure. Moreover, in this work we focus on the cold and dense regime and mostly simulate with N τ > 100, for which λ < ∼ 10 −25 , and terms that are of subleading order in N τ as well as terms proportional toh 1 are neglected, sinceh 1 → 0 as T → 0. For N f = 1 we then simulate the simplified action (2.48) For N f = 2 some care has to be taken when introducing the determinant for the second flavour, which introduces mixing terms that are not present in the above expression.

Hadron masses in strong coupling and hopping expansion
In order to interpret the results in the following sections, it is convenient to also have the leading order of the meson and baryon masses, To the orders given here, these expressions are the same for N f = 1, 2 degenerate masses. From the second equation we extract the running of the hopping parameter in the strong coupling limit for later use, 3 Analytic analysis of the effective theory 3.1 NLO perturbation theory for N f = 1 A lot of insight about the behaviour of the effective theory can be gained by studying the static strong coupling limit, where the partition function factorises into a product of one-link integrals which can be solved analytically. For the case of N f = 1 we previously observed the onset transition as a step function from zero density to lattice saturation [10]. Here we extend this analysis beyond the static strong coupling limit by using perturbation theory in the small couplings λ, h 2 , permitting a clear understanding how the nuclear liquid gas transition is driven by interactions.
To this end we consider the partition function with nearest-neighbour interaction between a Polyakov loop and its conjugate, as well as between two Polyakov loops, i.e. including the couplings λ, h 1 , h 2 . Here we are interested in the cold and dense regime. Near the zero temperature limit and for µ > 0, the anti-quark contributions vanish exponentially becauseh 1 → 0 for T → 0 and the simplified partition function is Note that the coupling h 1 parametrises (µ−m) and moreover approaches one around the onset transition. Therefore it cannot serve as an expansion parameter. On the other hand, there are physically interesting parameter regimes where λ, h 2 are sufficiently small to allow for a power series expansion. The leading orders for the partition function and pressure read In the cold and dense regime we are working with N τ ≥ 10 for which λ(β = 6.0, N τ ) < 10 −5 plays no quantitative role, so we neglect it from here on. The static strong coupling limit is obtained for λ = h 2 = 0 and has already been discussed in [10]. In this case the partition function factorises into one-link partition functions z 0 , i.e. it represents a non-interacting system. We identify z 0 to consist of baryons, a spin 3/2 quadruplet and a spin 0 baryon made of six quarks. Note that the Pauli principle for N f = 1 does not admit spin 1/2 doublets. The quark number density and the energy density then follow as where we have made use of eq. (2.50).

The nuclear liquid gas transition for
With these formulae at hand, it is easy to analyse the physics of the cold and dense regime. Let us begin with the static strong coupling limit. At high density, the lattice is populated until it is saturated with six quarks per lattice site according to the Pauli principle, lim Note that the dominating contribution to z 0 is a bosonic baryon. However, it is a composite of quarks such that the Pauli principle, built into the partition function in the original QCD action, is still contained in z 0 . Another limit of interest is that of zero temperature. In this case we have lim T →0 Thus we find the so-called silver blaze property, i.e. the thermodynamic functions stay zero as the chemical potential is raised until it crosses the constituent quark mass. Then it is possible to excite baryons and the onset phase transition to nuclear matter takes place. In the static strong coupling limit, this transition is a step function from zero to saturation density. This step function gets immediately smeared out to a smooth crossover as soon as a finite temperature (N τ < ∞) or coupling h 2 is switched on, cf. figure 2.
We can unambiguously identify this transition as baryon condensation by also looking at the energy density. Away from the static limit, there are non-vanishing attractive quark-quark (and hence baryon-baryon) interactions parametrised by h 2 . These are identified by the quantity which gives the energy per baryon minus its rest mass in units of m B . For temperatures approaching zero, this is the binding energy per baryon. In perturbation theory, the result is where we have used the leading order of eq. (2.49) to express the hopping parameter by the meson mass. This result beautifully illustrates several interesting physics points. In the non-interacting static limit with κ = h 2 = 0, there is no binding energy and hence no true phase transition for the onset to nuclear matter. For finite κ we see from the behaviour of z 3 , z 0 that for µ < m and T → 0 the binding energy is also zero, another manifestation of the silver blaze phenomenon. On the other hand, for µ > m, T → 0 it is explicitly negative and its absolute value increases with decreasing meson mass. This is in complete accord with expectations from nuclear physics models based on meson exchange. The binding energy as a function of chemical potential is shown in figure 3 (left), the asymptotes towards larger chemical potential are due to lattice saturation. Plotting against the number density, we obtain the equation of state as qualitatively expected for nuclear matter, figure 3 (right). In particular, the binding energy per baryon gets more negative as the quarks get lighter, until we see a minimum forming. Note that all curves eventually should turn upwards again, but for finite lattice spacing they are limited by the saturation density. With the explicit mass dependence in the binding energy and without a continuum extrapolation, quantitative predictions for physical QCD cannot be made until the physical flavour content and masses can be controlled. Nevertheless, it is interesting to keep in mind the physical binding energy per nucleon, ǫ ≈ 0.017 and the nuclear saturation density, n B0 /m 3 proton ≈ 0.016.

The static strong coupling limit for N f = 2 at finite baryon density
For β = 0, the partition function consists of the static determinant factors only We again consider the zero temperature limit at µ > 0, for which the anti-quark contributions vanish. After the gauge integration the result reads All exponents of h n u h m d come in multiples of three, n + m = mod 3. Just as in the one-flavour case (with h d = 0), this has the form of a free baryon gas where the prefactors give the degeneracy of the spin multiplets. Note that for N f = 2 we also find the standard spin 1/2 nucleons and many more combinations. To illustrate the prefactors, consider the example h 2 u h d . There is the spin 1/2 doublet, the proton, as well as a spin 3/2 quadruplet, the ∆ + , i.e. six states altogether. The states corresponding to h 2 d h u are the neutron and the ∆ 0 , while h 3 u , h 3 d are the ∆ ++ , ∆ − quadruplets, respectively. It continues with six-quark states. For example, h 4 u h 2 d has 10 allowed spin-flavour combinations, corresponding to three spin 1 triplets and one spin 0 singlet. These are consistent with an interpretation as di-baryon states built of ∆ ++ ∆ 0 or pp. Thus, eq. (3.11) contains all baryonic spin-flavour multiplets that are consistent with the Pauli principle, i.e. up to a maximum of 12 constituent quarks.
The quark density reads (3.12) In the high density limit numerator and denominator are dominated by the term with the highest power and we obtain lim µ→∞ (a 3 n) = 2 · 2 · N c ≡ 4(a 3 n B,sat ) .

(3.13)
This is the saturation density with two spin, two flavour and three colour states per lattice site and fermion. In the zero temperature limit we have again the silver blaze property followed by a transition to lattice saturation lim T →0 (3.14) 3.4 The static strong coupling limit for N f = 2 at finite isospin density Finite isospin density is realised for µ I = µ u = −µ d [25]. Choosing µ u > 0, the zero temperature limit impliesh u , h d → 0 for T → 0. Omitting the corresponding terms from eq. (3.11) and performing the gauge integration we find the expression With isospin chemical potential, d-anti-quarks are now playing the same role as uquarks and the partition function is a free gas of baryons, anti-baryons and mesons. Differentiating with respect to isospin chemical potential gives the isospin density, . Also in this case, we find saturation in the high density limit, lim µ→∞ (a 3 n I ) = 2 · 2 · N c ≡ 4(a 3 n I,sat ) . (3.17) Just as in the case of finite baryon density, the high density expression is dominated by a bosonic composite state which "knows" that it consists of constituent quarks, of which only a finite number can be packed on one lattice site. The saturation level is hence identical to that for large baryon chemical potential. Similarly, in the zero temperature limit we find again the silver blaze property followed by a non-analytic transition to isospin condensation, lim T →0 Note that for static quarks, m B /3 = m π /2 and the onset transition to nuclear or isospin matter fall on top of each other as a function of quark chemical potential. We shall see in our numerical investigations that a gap between them opens up as expected when interactions between the hadrons are switched on.

Simulation of the effective theory by complex Langevin
The effective theory specified in the last sections has a sign problem. With less degrees of freedom and the theory being only three-dimensional, the sign problem is milder than in the original theory such that Monte Carlo methods are feasible at finite temperatures and chemical potentials µ/T < ∼ 3 [9]. If, however, one is interested in cold dense matter in the zero temperature limit, the sign problem becomes strong and Monte Carlo methods are restricted to small volumes. Fortunately, the effective theory is amenable to simulations using complex Langevin algorithms (for an introductory review, see [26]) and the onset transition to nuclear matter could be demonstrated explicitly for very heavy quarks [10]. In this section we discuss the validity of complex Langevin for the effective theory. We will only sketch the general method here, as there is an abundant literature on this subject [26][27][28][29][30]. The basic idea is to introduce a fictitious Langevin time θ, in which a field theoretical system with a generic field φ evolves according to the Langevin equation where η(x, θ) denotes Gaussian noise. In the case of a complex action, the field variables have to be complexified too, φ → φ r + iφ i . In our case, the degrees of freedom of the effective theory are the temporal Wilson lines We may further simplify this by taking the trace of the Wilson lines and parametrising the resulting Polyakov loops in terms of two angles, bringing them into a diagonal form [31] L(θ, φ) = e iθ + e iφ + e −i(θ+φ) .  Hence the integration measure we use in our simulation is the reduced Haar measure [dφ] e 2V . This means instead of an integration over SU(3) matrices we have 2 complex degrees of freedom on every spatial lattice point. Furthermore, having only diagonal matrices their inversion is trivial. With these ingredients eq.(4.1) was solved numerically using stepsizes of around ǫ = 10 −3 and applying the adaptive stepsize technique proposed in [32] to avoid numerical instabilities.

Criteria for correctness
It is well known that the complex Langevin algorithm is not a general solution to the complex action problem since it converges to the wrong limit in some cases, including some parameter ranges for QCD [26,33]. The failure can be attributed to insufficient localisation of the probability distribution in the complex field space, and a set of criteria was developed to check whether this localisation is sufficient in a given simulation [29]. A necessary condition is that the expectation value of all observables O[φ] vanishes after a Langevin operatorL has been applied to them, .
While, strictly speaking, this test is necessary on all observables of the theory, in practice only a select few can be tested. Note that in the framework of our effective theory, all observables are expressed as functions of Polyakov loops and one might hope that its proper behaviour propagates to more complicated functions of it. In figure 4 we show the expectation value of the Polyakov loop as a function of the step size of the Langevin algorithm for the effective theory to order κ 2 (left) and κ 4 (right). In both cases the criterion is satisfied in the limit of vanishing stepsize.

The logarithm of the static determinant
Another problem related to the distribution in the complexified field space has recently been pointed out for all partition functions containing a complex determinant [34]. Its contribution to the effective action is ∼ log det, and the logarithm develops a cut along the negative real axis, i.e. it is multi-valued. This may cause a problem whenever the calculation of the drift term for the Langevin time requires a derivative to be taken across the cut. In [34] it was found for a random matrix model that these crossings lead to incorrect predictions for observables if they happen frequently in a Monte Carlo history. Here we can see another benefit of the effective theory compared to a Langevin simulation of full QCD. In the effective theory, only the static determinant features this problem, while the corrections to the effective action in the hopping expansion are exponentials of polynomials. We have monitored the static determinant during the Langevin evolution, an example is shown in figure 5 at baryon density slightly below (left) and above (right) the onset transition to nuclear matter. Note that the static determinant is dominated by the Polyakov loop. One observes the expected distortion from the centre symmetric distribution of the vacuum state to the distribution preferring the real centre sector, and this distortion is amplified exponentially with chemical potential. For simulation purposes, the crucial observation is that there are nearly no crossings of the negative real axis, in accord with the satisfied convergence criterion above. We have monitored such scatter plots over a wide range of parameter values. Occasionally crossings of the negative axis do occur, but the observed frequency was < 10 −4 in all cases.  possible, i.e. on small volumes. In [10] we have shown a successful comparison for very small hopping parameters κ ∼ 10 −4 . Figure 6 shows that this test is also passed for significantly larger values κ ∼ 0.01. We conclude that complex Langevin simulations of the effective theory constructed here are fully controlled for the entire coupling range investigated, 0 < β < 6 and 0 < κ < 0.12. This is an algorithmic advantage over Langevin simulations in full QCD, where gauge cooling techniques [35] are required to control the field distribution and even then simulations at small lattice couplings are ruled out [7].

Convergence region of the hopping series
An important task is to find the region of validity of the effective theory. By this we mean the region, determined by a self-consistent test, where the truncated effective theory is a good approximation to the full theory. As criteria we choose the difference between expectation values of observables, calculated from the κ 2 and the κ 4 action, O κ 2 , O κ 4 . These can be evaluated as a function of the expansion parameter Nc , and the convergence region is where the difference is smaller than the desired accuracy. Since we are interested in the onset of baryon number, we choose the density in lattice units a 3 n as an observable and compute it at a fixed value of the coupling h 1 = 0.8. As can be seen in figure 7, the static limit is only a valid approximation in the κ → 0 limit. Note that the resummed action offers a slightly better convergence. Therefore, we will use this version for our simulations. The expansion parameter already shows that the region of convergence is limited in the direction of low temperatures and light quarks, i.e. one can reach lower quark masses at larger temperatures.

Setting a scale and units
Setting a scale and performing continuum limits along lines of constant physics is a computationally very demanding task. Rigorously speaking, this is truly possible only at or near the physical point. On the other hand, the effective theory considered here is only valid for very heavy quarks, due to the truncated hopping series. While it exhibits most qualitative features of physical QCD, its spectrum is still far from the experimentally observed one. For this reason we do not attempt to accurately fix our hadron masses. (In the mass ranges considered this would anyway demand heavy quark effective theories [36]). Instead we only provide a very rough guide where we are in parameter space.
Our procedure is as follows: heavy quarks have little influence on the running of the coupling. Thus we use the non-perturbative beta-function of pure gauge theory for the lattice spacing in units of the Sommer parameter, a(β)/r 0 [37]. Taking r 0 = 0.5 fm sets a physical scale for our lattices, while N τ tunes temperature via T = (aN τ ) −1 . In a very rough approximation we then use the strong coupling expressions eq. (2.49) for the hadron masses.

The nuclear liquid gas transition in heavy dense QCD
In our previous work [10] we performed a continuum extrapolation for the transition to cold nuclear matter based on the κ 2 action. In figure 7 we repeat this calculations including the κ 4 corrections. This allows us to simulate smaller lattice spacings a = 0.08 fm without leaving the region of convergence, since reducing a while keeping m B /T and T fixed means going to higher κ and N τ . Nevertheless, the extrapolation suffers from considerable uncertainties, resulting in large errors in the high density phase. This can be seen in fig. 8 (left), where we show the two best fits for our data at µ B /m B = 1 at several lattice spacings. This is the chemical potential where different extrapolation fits differ the most. The systematic truncation error for our κ 4 data is estimated as the difference to the data obtained from the κ 2 action and included in the error bars in the figure. This data was then fitted to get a value for a → 0. For each value of the chemical potential we tried several fits (linear and quadratic) with one to three degrees of freedom. For the best fits we always achieved χ 2 red < 2 as long as µ B /m B < 1.0014. For the continuum result we quote the average of the two best fits, the error was estimated as difference between those two fits. We note that the results at κ 4 are somewhat higher than our κ 2 -results in [10]. This is because inclusion of κ 4 is the first order allowing for a realistic estimate of the truncation error, and thus permits inclusion of data with smaller lattice spacing. This results in the continuum extrapolated baryon number density in figure  8 (right), where we display the results for N f = 1, 2 for a temperature T = 10 MeV. In the low density region the "silver blaze" property, i.e. the independence of the thermodynamic functions of chemical potential can be seen. The growing uncertainties in the high density region are caused by the unphysical saturation on the lattice which limits the density to 2N f N c quarks per lattice site, while in the continuum no such saturation exists. As expected, the onset of nuclear matter happens at a critical value µ c B < m B , due to the nuclear binding energy. The location of the onset suggests a very small binding energy ∼ 10 −3 m B for the heavy quarks considered here, in accord with our perturbative analysis, section 3.2. This explains why the onset transition is a smooth crossover rather than the first-order transition expected for light quarks. The endpoint of the nuclear liquid gas transition sits at a temperature of the order of the binding energy and is not visible for very heavy quarks. In accord with expectation, the onset with two flavours is steeper than with one flavour.
It is now straightforward to compute the other thermodynamic functions and from them the equation of state. Figure 9 shows the pressure as a function of baryon chemical potential as well as a function of baryon density, whereas the binding energy per nucleon is shown in figure 10. Note that in all plots the error bars include the systematic uncertainty of both, the truncation of the effective theory as well as the continuum extrapolation. The plot of the binding energy is particularly intriguing. For small density it is zero, another manifestation of the silver blaze property, until it turns negative, thus causing the condensation of nuclear matter. At larger density, lattice saturation is reached before the expected upturn of the curve. Nevertheless, the shape of the curve suggests that the minimum has been reached near the right border. Its numerical value of the order of 10 −3 is consistent with that observed from the location of the onset transition in figure 8 (right).

Nuclear liquid gas transition for light quarks
As in our previous work [10], the accessible quark masses in the convergence region of the effective theory are too high to realise the expected first order transition for the onset of nuclear matter. Finite size scaling analyses reveal the transition to be a smooth crossover, in accord with the interplay between accessible temperatures and the values of the binding energies. Of course it is highly interesting to see whether the effective theory includes the expected physics features when the quark mass is lowered. We now consider κ = 0.12, corresponding to a small quark mass, and very low temperatures parametrised by N τ ∼ O(10 3 ). We stress that this choice of parameters is far outside the convergence region of our κ 4 -action, cf. figure 7. In other words, there is no reason to expect the results to accurately represent QCD and an attempt at a continuum extrapolation makes no sense. Nevertheless, this is an interesting check of the qualitative features of the effective theory. Figure 11 shows distributions of the Polyakov loop in the onset transition region for three choices of N τ , corresponding to increasing temperatures from left to right. We clearly observe the coexistence of two phases at the lowest temperatures, which indicates a first order transition between them. As the temperature is raised (N τ is lowered), the two-state signal weakens and merges to a single gaussian distribution, signalling a weakening and eventual disappearance of the first-order transition. This picture is corroborated by a finite size analysis of the quark number susceptibility in figure 12. First-order and crossover transition are clearly distinguished by diverging and finite susceptibility as a function of volume. Thus we conclude, while our κ 4 -action used in this work is not quantitatively reliable in this parameter range, it displays all the qualitative features expected for the nuclear liquid gas transition: a first-order transition from the vacuum to nuclear matter which weakens with temperature until it vanishes in a critical endpoint. We therefore expect higher orders in the effective action to only correct the quantitative details of this transition.

Isospin vs. baryon chemical potential
Let us finally consider the situation in the two-flavour theory with finite isospin chemical potential, µ I = µ u = −µ d . In section 3.4 we have discussed the situation in the static strong coupling limit, where the onset transition for pion condensation at µ I = m π /2 happens at the same chemical potential as the one for baryon condensation at µ B = µ B /3. With interactions included, this gets modified in two ways. Firstly, we have m π /2 < m B /3 in this case, and secondly the onset gets shifted to smaller chemical potentials by the non-vanishing binding energy. The first effect also leads to the expected gap opening between the onset to pion condensation vs. that to baryon condensation [38], when plotted against quark chemical potential, as shown in figure 13.

Conclusions
In this work we further elaborated the construction of an effective three-dimensional lattice theory for QCD thermodynamics. It is formulated entirely in terms of Polyakov loops and calculated from the 4d Wilson action as a strong coupling and hopping series which is now complete to order κ n u m , (n + m) = 4. In the static strong coupling limit, the effective theory can be solved exactly, providing the complete spin-flavour structure of the hadron spectrum as well as an onset transition from zero density to lattice saturation. The interacting effective theory has a sign problem that can be handled by complex Langevin simulations with fully satisfied convergence criteria. Moreover, the sign problem is mild enough that on small volumes Monte Carlo simulations are feasible, even at real chemical potential. The couplings of the effective theory are sufficiently small to also permit a perturbative evaluation, which agrees with numerical results in wide regions of the parameter space. Altogether this allows for a controlled and very efficient evaluation of thermodynamic functions and critical couplings.
Working in the heavy quark region near the static limit, where continuum extrapolations of thermodynamic functions are feasible, we have explicitly demonstrated the onset transition to cold nuclear matter and calculated the nuclear equation of state for the first time directly from QCD. In particular, we find a negative binding energy per nucleon as the expected reason for baryon condensation. In accord with expectations from models of nuclear interactions, the binding energy is governed by exponentials of the meson mass and suppressed for heavy quarks. Decreasing the quark mass beyond the convergence region of our expansion, we indeed observe the nuclear onset transition to emerge as a first order liquid gas transition with an endpoint at some small temperature. In this parameter range also the expected gap opens up between the onset of pion condensation in the case of finite isospin chemical potential and the nuclear onset at finite baryon density.
In summary, the effective lattice theory described in this work contains all the qualitative physics expected for cold nuclear matter. It remains to be seen whether high enough orders of the hopping expansion can be generated in the future in order to reach physical quark mass values. However, since the hopping convergence is much faster at high temperatures, the current effective theory might already be useful to describe the finite temperature phase structure of QCD with light quarks. Work in this direction is in progress.

A Wilson line contributions to the effective action
In this appendix we list final expressions for all types of terms appearing in the kinetic determinant derived in section 2.5. Tr Tr Tr Tr The TrP i M i P j M j and TrP i M j P j M i contributions are the same, just with different directions and a factor of 1 2 in front.