QCD in the heavy dense regime for general $N_c$: On the existence of quarkyonic matter

Lattice QCD with heavy quarks reduces to a three-dimensional effective theory of Polyakov loops, which is amenable to series expansion methods. We analyse the effective theory in the cold and dense regime for a general number of colours, $N_c$. In particular, we investigate the transition from a hadron gas to baryon condensation. For any finite lattice spacing, we find the transition to become stronger, i.e. ultimately first-order, as $N_c$ is made large. Moreover, in the baryon condensed regime, we find the pressure to scale as $p\sim N_c$ through three orders in the hopping expansion. Such a phase differs from a hadron gas with $p\sim N_c^0$, or a quark gluon plasma, $p\sim N_c^2$, and was termed quarkyonic in the literature, since it shows both baryon-like and quark-like aspects. A lattice filling with baryon number shows a rapid and smooth transition from condensing baryons to a crystal of saturated quark matter, due to the Pauli principle, and is consistent with this picture. For continuum physics, the continuum limit needs to be taken before the large $N_c$ limit, which is not yet possible in practice. However, in the controlled range of lattice spacings and $N_c$-values, our results are stable when the limits are approached in this order. We discuss possible implications for physical QCD.


Introduction
The QCD phase diagram is important for many aspects of current nuclear, heavy ion and astro-particle physics, yet it remains largely unknown. This is because lattice QCD at finite baryon chemical potential has a severe sign problem, which prohibits straightforward Monte Carlo simulations. Various workarounds to extend the Monte Carlo method introduce additional approximations and are limited to the high temperature and/or low density region, with baryon chemical potential µ B /T < ∼ 3 [1]. No sign of criticality is found in this region, where the transition from a hadronic gas to a quark gluon plasma proceeds by an analytic crossover [2][3][4][5]. The same conclusion is reached by analytic non-perturbative approaches like Dyson-Schwinger equations [6] or the functional renormalisation group [7]. Despite continuing efforts, a genuine solution to the sign problem, and hence fully nonperturbative access to the cold and dense region of QCD, are missing to date.
This situation has motivated the study of QCD also in unphysical, but controllable parameter regions, where the sign problem can be overcome by either algorithmic or analytic methods. In this work our goal is to bridge two such approaches. We employ an effective lattice theory derived from the standard Wilson action by combined strong coupling and hopping (inverse quark mass) expansions. The effective theory is valid on reasonably fine lattices, as long as the quarks are sufficiently heavy to be described by the next-to-next-to-leading order in the hopping expansion. In this parameter range the theory correctly reproduces the critical heavy quark mass at zero density, where the first-order deconfinement transition changes to a smooth crossover [8], and furthermore allows for an extension to finite baryon chemical potential, including the cold and dense regime around the onset transition to baryon matter.
Here we consider the effective theory for general colour gauge group SU (N c ), in order to establish contact with another effective approach in the continuum, namely QCD at large N c . In particular, we analyse thermodynamic functions around the onset transition to baryon matter in the cold and dense regime, for varying and large N c . This allows us to address, by direct calculation, various conjectures made in [9] regarding the phase diagram and the effective degrees of freedom at large N c . There, the authors argue for the existence of quarkyonic matter, which is characterised by its pressure scaling as p ∼ N c and has both baryon-like and quark-like aspects. Phenomenological consequences of this form of matter in physical QCD have been assessed in [10][11][12][13]. A general, qualitative discussion about the possibilities for the phase diagram in (T, µ, N c )-space as well as references to earlier work can be found in [14].
We begin with a brief review of the effective lattice theory in section 2. This material is not new and can be skipped by readers familiar with it, but is needed as reference point when interpreting the following results. In section 3, a summary of the conjectured phase diagram at large N c is followed by our proper calculations for general N c and the analysis of the results for large N c . Finally, section 4 concludes what is expected for physical QCD.

The effective lattice theory
Consider the partition function of lattice QCD with the standard Wilson action at finite temperature, T = 1/(aN τ ), realised by compact euclidean time with N τ slices and (anti-) periodic boundary conditions for (fermions) bosons. An entirely equivalent formulation in terms of temporal lattice links only is obtained after performing the Gauss integral over the quark fields and integrating the gauge links in spatial directions, (2.1) With the spatial links gone, the effective action depends on the temporal links only via Wilson lines closing through the periodic boundary, For SU (2) and SU (3) the effective action can always be expressed in terms of Polyakov loops, L(x) = TrW (x), whereas for larger N c in general traces of higher powers of W appear as well.
This effective action is unique and exact. However, the integration over spatial links causes long-range interactions of Polyakov loops at all distances and to all powers so that in practice truncations are necessary. For non-perturbative ways to define and determine truncated theories, see [15][16][17][18]. Here, we use an effective theory based on expanding the path integral in a combined character and hopping parameter series. Both expansions result in convergent series within a finite radius of convergence (for an introduction, see [19]). Truncating these at some finite order, the integration over the spatial gauge links can be performed analytically to provide a closed expression for the effective theory. Going via an effective action results in a resummation to all powers with better convergence properties compared to a direct series expansion of thermodynamic observables as in [20][21][22]. Since the Wilson line W (x) contains the length N τ of the temporal lattice extent implicitly, the effective theory is three-dimensional. Note that, in the case of 4d Yang-Mills theory, this representation by a 3d centre-symmetric effective theory is the basis for the Svetitsky-Yaffe conjecture [23,24] concerning the universality of SU (N c ) deconfinement transitions. Including the quark determinant via the hopping expansion introduces centre symmetry breaking terms and additional effective couplings.
Let us briefly summarise the expansions used in order to perform the spatial link integrations. The gauge part of the action is a class function with respect to the product of the links of one plaquette: where V ∈ SU (N c ). Therefore it can be expanded in the characters χ r of the irreducible representations r of SU (N c ) at every plaquette, In this formula, d r denotes the dimension of the representation and a r (β) is the character expansion coefficient divided by the expansion coefficient of the trivial representation. The expansion coefficients can be computed exploiting the orthogonality of the characters: We drop the overall factor of c 0 in equation (2.4), as it cancels in expectation values. For the integration of the spatial links following this expansion, one can use the formulas for those cases where not more than two non-trivial representations share a common link. In earlier publications [25] this was used to derive the effective gauge action for SU (3) to rather high orders in the coefficient of the fundamental character u(β) ≡ a f (β)/d f = a f (β)/N c . The coefficients of the higher dimensional representations can be expressed in terms of the fundamental one, see also section 3.1, and therefore the expansion can be organised according to powers of the fundamental character. The dependence of u on the lattice gauge coupling β = 2N c /g 2 can be specified either as a power series or numerically, It is known to arbitrary precision, and u is always smaller than one for finite β-values.
For the hopping expansion it is useful to split the quark matrix in temporal and spatial hops between nearest neighbours, This is because the static determinant containing all temporal hops (and only those) can be computed exactly [26,27]. We then do the hopping expansion of the kinetic quark determinant using det Q kin = exp(Tr(ln(Q kin ))) . (2.13) This leads to an expansion in powers of S, which is proportional to the hopping parameter, (2.14) The expansion terms are then ordered according to their number of spatial hops while the temporal ones are resummed to all orders. Since the hopping expansion is in inverse quark mass, the effective theory to low orders is valid for heavy quarks only. For a derivation of the effective theory to order O(κ 4 ) in spatial hops, see [27]. In this case the relevant integrals for the fermionic contributions are Generically, the effective action obtained in this way has the following form: The λ i are defined as the effective couplings of the Z(N c )-symmetric terms S s i , whereas the h i multiply the asymmetric terms S a i . In particular, h 1 ,h 1 are the coefficients of L, L * , respectively, and to leading order correspond to the fugacity of the quarks and anti-quarks, Here, is the constituent quark mass in lattice units of a baryon computed to leading order in the hopping expansion [28], while is the effective coupling of a nearest neighbour L x L y interaction. The partition function for SU (3), including just these simplest interactions, reads In this expression the first line represents the pure gauge sector, the second line is the static determinant and the third line the leading correction from spatial quark hops. This partition function has a weak sign problem and can be simulated with either reweighting or complex Langevin methods [8,27]. Since the effective couplings correspond to power series of the expansion parameters, they are themselves small in the range of validity. Hence, the effective theory can also be treated by linked-cluster expansion methods known from statistical physics, with results for thermodynamic observables in quantitative agreement with the numerical ones [29]. In this way, full control over the sign problem is achieved.

The deconfinement transition
The phase diagram of QCD with heavy quarks is depicted schematically in Fig. 1. At zero density, the thermal transition is a first-order deconfinement transition. It is a remnant of the centre symmetry-breaking transition of the SU (3) pure gauge theory, which gets weakened by explicitly centre-breaking finite quark masses ∼ 1/m q , until it ends in a second-order point for some critical mass m c q . In the effective theory, this phase transition appears as spontaneous breaking of the Z(3)-symmetry at some set of critical couplings , which can be determined by numerical simulation. Inversion of the effective couplings then gives predictions for β c (N τ ), κ c (N τ ), which can be compared with the results from full QCD simulations.  For SU (2) and SU (3)-Yang-Mills theory, the simplest effective theory with only a nearest neighbour coupling (first line in equation (2.22)) correctly reproduces the universality of the respective deconfinement transitions, and the predicted β c (N τ ) are within 10% of their true values for N τ = 2, . . . , 16 [25]. For QCD with heavy quarks, the simplest effective theory including the static determinant and κ 2 -corrections, predicts κ c to better than 10% on N τ = 4 [8]. Contrary to full QCD, the effective theory can be simulated at finite chemical potential to determine the location of the critical end point as a function of quark mass [8]. This qualitative behaviour of the deconfinement transition in the heavy quark region is also found by continuum studies using a Polyakov loop model [30] and in the functional renormalisation group approach [31].

The onset transition to finite baryon number
Going out along the chemical potential axis at low temperature, the system crosses the onset transition beyond which the ground state consists of condensed baryon matter. In order to interpret our analysis for general N c , let us first recall the situation for N c = 3 in some detail. The qualitative features are best understood in the strong coupling limit, β = 0, with the static quark determinant only, where (2.22) is reduced to the second line. The partition function then factorises into one-site integrals which can be solved analytically. Since we are interested in low temperatures, where mesonic contributions are exponentially suppressed by their fugacity factors, we simplify the analysis by settinḡ h 1 = 0. For N f = 1 the partition function then reads [27,32] corresponding to a free baryon gas with two species. With one quark flavour only, there are no nucleons and the first prefactor indicates a spin 3/2 quadruplet of ∆-baryons whereas the second term is a spin 0 six-quark state or di-baryon. The quark number density is (2.24) and at zero temperature exhibits a discontinuity when the quark chemical potential equals the constituent mass m. This reflects 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 [33]. Once the baryon chemical potential µ B = 3µ is large enough to make a baryon (m B = 3m = m LO B in the static strong coupling limit), a discontinuous phase transition to a saturated crystal takes place. Note that saturation density here is 2N c quarks per flavour and lattice site and reflects the Pauli principle. This is clearly a discretisation effect that disappears in the continuum limit.
For two flavours the corresponding expression for the free baryon gas reads where we have now distinguished between the h 1 coupling for the u-and d-quarks. In this case we identify in addition the spin 1/2 nucleons as well as many other baryonic multiquark 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 [27]. Again, the onset transition to finite baryon density is a step function from zero to saturation density, which now is 2N c N f quarks per site. The step function behaviour gets immediately smeared out to a smooth crossover, as soon as a finite temperature, N τ < ∞, is switched on. This implies that the firstorder line of the nuclear liquid gas transition is exponentially short as a result of the large quark masses, as expected from nuclear physics Yukawa potentials with meson exchange. Indeed, the interaction energy per baryon, which sets the scale for the critical endpoint of the nuclear liquid gas transition, can be extracted from the dimensionless combination of energy density and baryon number of the system, (2.26) in the limit of zero temperature. In the strong-coupling limit one finds Thus the length of the liquid gas transition in Fig. 1 is a function of quark mass and decreases to zero towards the static limit. Including κ 4 -corrections, a first-order transition ending at some finite temperature is explicitly seen [27]. Including the gauge coupling, and with sufficiently many corrections at hand, also the lattice spacing can be varied and the approach to the continuum can be studied. In [29] it was shown that through orders u 5 κ 8 for a sufficiently heavy quark mass, continuum-like behaviour is obtained immediately after the onset transition, with the qualitative features discussed here.

QCD for large N c
Since in the framework of the effective theory we can work fully analytically, it is straightforward to investigate what happens when the number of colours N c is varied and made large. In particular, we aim to explore some large N c considerations leading to the prediction of quarkyonic matter [9].
There is a lot of interesting literature on QCD at large N c , which we are unable to represent properly. In particular, baryon matter in the combined heavy quark and large N c limits has been considered by a mean field analysis in the continuum [34,35]. Here, our approach is quite different in working on the lattice with large but finite quark masses, for general N c and beyond mean field.
The essential qualitative features of QCD at lage N c were established in the early works [36,37]. The 't Hooft limit is defined by In this case the theory has the following properties: The authors of [9] used these and various other ingredients to draw qualitative conclusions for the QCD phase diagram. Fig. 2 shows their conjectured phase diagram in the large N c limit. From finite temperature perturbation theory it follows that, in the plasma phase, p ∼ N 2 c . With quark loops suppressed, the phase boundary of the deconfinement transition is pure gauge-like and unaffected by chemical potential. It thus forms a horizontal line, staying first-order everywhere. On the other hand, in the hadronic, low density phase, thermodynamics is quantitatively well described by a weakly interacting hadron resonance gas [38,39]. Statistical mechanics then implies that the baryonic contribution to the pressure is exponentially suppressed with baryon mass, so p ∼ N 0 c there. In [9] a similar combination of perturbative (valid for large µ) and statistical mechanics arguments for baryons suggests that, for low temperatures and µ B > m B , the pressure scales as p ∼ N c . The authors termed this phase "quarkyonic", since it shows aspects of both quark matter and baryon matter. In particular it is argued that, for zero temperature, excitations relative to the Fermi sea should be baryon-like for (p − p F ) < ∼ Λ QCD and quark-like for (p − p F ) Λ QCD , implying a shell structure in momentum space as in Fig. 2 (right). Since the Fermi momentum is p F < ∼ Λ QCD at the onset transition and then grows with quark chemical potential, this picture suggests the possibility to smoothly interpolate between baryon matter (right after the onset) and quark matter (at very large densities). We will now address these issues by direct calculations using the effective lattice theory.

The effective lattice theory for general N c
Note that N c = 2 has already been analysed in detail [40], with interesting physics results for two-colour QCD. Our aim here is to go in the other direction and to increase N c . For the gluonic part, the derivation of the effective theory for general and large N c has already been discussed in [41]. Note that the integration rules (2.7) and (2.8) are true for arbitrary N c . Therefore, in cases where only these formulas are relevant, one simply has to replace d r and a r by their appropriate generalisations to N c . Specifically, the character expansion coefficients a r can be obtained via [42] (3.2) In this formula, the α i are a set of N c positive descending integers with α Nc = 0 which label the representation and correspond to Young tableaux. Following [42] one may re-express all coefficients of higher representations in terms of the fundamental representation using double Young tableaux. The characters corresponding to a double young tableau can be determined in terms of the traces of powers of U and U † [43].
In [41] the fermionic contributions were also discussed using the hopping expansion. However, only a subset of spatial hoppings to O(κ 4 ) was considered, and temporal hoppings were included up to O(κ 2Nτ ). As we mentioned in section 2, we work in a scheme where temporal hoppings are resummed to all orders. Nevertheless, equations (2.15) and (2.16) are valid for general N c so the fermionic contributions obtained in this way are legitimate also for general N c . Spatial baryon hoppings contribute at O(κ Nc ) and therefore they are suppressed for large N c . To evaluate the contributions of meson hoppings, integrals of the type are needed. Since these integrals give the same result for U ∈ U (N c ) and U ∈ SU (N c ) one can, for the fermionic contributions in the vacuum and at finite temperature, work with U (N c ) instead of SU (N c ) at large N c . Likewise, when analysing the pure gauge theory for large N c , the same simplification applies [44]. However, at finite chemical potential temporal quark hoppings in the positive direction are boosted by a factor of e aµ and therefore integrals of the type are relevant. These integrals vanish for U (N c ) when b = a, but not for SU (N c ), where they contain baryonic contributions. For the cold and dense regime, we therefore need to calculate integrals over SU (N c ).

Evaluation of the SU (N c ) effective theory in the strong coupling limit
We begin our analysis in the strong coupling limit, u(β = 0) = 0 and thus λ 1 = 0, h 1 = h LO 1 , and employ the effective theory to order O(κ 4 ) in spatial hoppings. We also neglect terms containingh 1 , since they are exponentially suppressed at low temperatures. The theory with these approximations already shows the most salient features of baryon dynamics, and we discuss modifications by the neglected couplings in later sections. Thus we get the free energy density where we have introduced the notation The required integrals are related in the following way, z (22) = z (11) − z (21) , (3.10) (21) .
Therefore, we only need to integrate z 0 , which corresponds to the integration over the static determinant, and z (21) . Note that all integrands are class functions of SU (N c ) group elements, f (W ) = f (V W V −1 ), which are invariant under a change of basis. These functions only depend on the eigenvalues z i of a group element. Furthermore, for our purposes it is sufficient to specialise to functions which factorise in the following way For such functions, the integration over the group can be expressed as [45] SU (Nc) (3.13) Using this formula one obtains (

3.15)
To evaluate the occurring determinants we showed, using the techniques explained in [46], that det 1≤i,j≤N where we have introduced the underline notation for the falling factorials n k = n · (n − 1) · · · (n − k + 1). (3.17) Applying this formula results in det 1≤i,j≤Nc With the free energy density at hand, it is now possible to compute all thermodynamic functions for any desired value of N c within the framework of our hopping expansion. Specifically we use the well known thermodynamic relations for the pressure baryon number density and energy density The derivative of κ with respect to a is computed at constant baryon mass, which is given to first order in κ by (2.20) resulting in a ∂κ ∂a = κ ln(2κ). (3.25)

The onset transition for increasing N c
Of particular interest is the behaviour of the onset transition to finite baryon density, which is shown in Fig. 3 for different choices of N c . We observe a steepening of the transition with increasing N c , which asymptotically ends up in a step function, i.e. a first-order transition, even though we started with a smooth crossover at N c = 3.
Decreasing the values of N τ , i.e. increasing the temperature, flattens the curves with fixed N c , but for asymptotically large N c a step function is obtained for any finite starting value of N τ . Thus, growing N c appears to make the onset transition to baryon matter always first-order. (With increasing temperature one may question the neglect of λ 1 ,h 1 . Their inclusion is discussed in sections 3.6, 3.7.) Note that, in the limit of infinite N c , the transition is between the vacuum and a saturated lattice, similar to what happens in the static strong coupling limit at finite N c . This saturated state is a discretisation artefact and will move towards infinity as the continuum is approached, as we discuss in section 3.8.

Thermodynamic functions for large N c
Since we have an explicit formula for the free energy for general N c , one can easily obtain the asymptotic behaviour of thermodynamic observables for large N c . We study the behaviour of the different orders in the hopping expansion separately. This is necessary because, beyond the onset of baryon condensation, the leading static term represents lattice saturation, which is an unphysical artefact of discretisation. As discussed in section 2.3, correction terms do not contribute to saturation, but modify the shape of the curves as they enter their low and high density asymptotes. These effects will remain after continuum extrapolation and thus are physically significant. The general strategy for the asymptotic analysis is most easily illustrated for the leading order contribution to the pressure at N f = 1 Note that, just like in the SU (3) case in (2.23), the prefactors before h Nc 1 can be understood from spin-degeneracy. Specifically, a colour neutral state consisting of N c fermions is antisymmetric in colour space under particle exchange. The only completely symmetric spin state of N c spin 1/2-particles is that with s = N c /2 . States with this spin and spin components −N c /2 ≤ s 3 ≤ N c /2 are degenerate, explaining the N c + 1 prefactor.
When h 1 < 1 then the term h Nc 1 is strongly suppressed (stronger than N k c can grow for any k) and a Taylor expansion around h Nc 1 = 0 gives For h 1 > 1 the term with the highest power of h Nc 1 determines the asymptotic behaviour and one obtains Order hopping expansion For higher N f and higher orders the terms become more complicated, but the general behaviour stays the same. For N f = 2 our findings on both sides of the onset transition are summarised in Table 1. A clear picture emerges: for h 1 < 1 all terms, due to the static determinant as well as the corrections, come with a factor h Nc 1 to some power. Since the fugacity contains m LO B /(N c T ) in the exponential, this factor will for low temperatures always dominate the powers of N c and result in a stronger exponential suppression before the onset transition. In other words, the curves for all quantities will be squeezed ever more tightly against the chemical potential axis as N c gets large. Since we do not know the hopping expansion of the baryon mass for general N c , we expressed our results in units of the leading order expression (2.20), which is responsible for onset happening at m LO The more interesting situation is h 1 > 1, where we first focus on the baryon number density. As explained in Section 2.3, the leading order contribution in the hopping expansion corresponds to the static determinant only, for which the onset to baryon matter is a first-order step function. This remains true for large N c , with the lattice quark saturation density going as a 3 n sat = 2N f N c , i.e. the baryon density behaves as a 3 n sat B ∼ const.. The most intriguing result of this section is the N c -scaling of the pressure beyond baryon onset, p ∼ N c . Preliminary results based on leading and next-to-leading order were reported in [47]. Stability of this finding through three orders suggests it to hold to any order in the hopping expansion, and thus for all current quark masses. In this case strongly coupled QCD beyond the onset transition to baryon matter satisfies the definition of quarkyonic matter [9]. Note that there is a finite interaction energy per baryon in units of baryon mass also for N c → ∞, as was conjectured in [9]. Its value at leading order κ 2 is determined by d(d + 1)/2, where d refers to the number of spatial dimensions.

The transition region
The results in the previous section were obtained by first fixing h 1 , i.e. fixing the quark chemical potential, and then considering the limit N c → ∞. Right around the onset transition one can also consider h Nc 1 ∼ 1. According to equation (2.18), this means that the quark chemical potential is adjusted such that µ − m ∼ 1/N c , where am = ln(2κ) is again the leading order constituent quark mass. In this regime, the asymptotic behaviour is dominated by the prefactors of the powers of h Nc 1 , which are polynomials in N c , see equations (3.14), (3.15), (3.16) and (3.19). For large N c one then obtains for the pressure in the hopping expansion This indicates that the hopping expansion does not converge for large N c . In Fig. 3 this shows up by the formation of uncontrolled wiggles in the central region for sufficiently large N c (a first indication of this is seen for N c = 9). Thus, we cannot make a statement about the large N c behaviour in this window of parameter space. Fortunately, the width of this region shrinks to zero for large N c and does not affect the observations in the previous sections. (For the same reason, it is not clear to us whether this behaviour is related to a phase transition in N c , as conjectured in [48]).

Inclusion ofh 1 -corrections
For statements at higher temperatures, or lower N τ , we also need to consider what happens whenh 1 is included. In this case terms appear which mix h 1 andh 1 . For h 1 > 1 the large N c analysis is clearly unchanged, since in this case the highest powers of h Nc 1 determine the asymptotic behaviour. For h 1 < 1 µ-independent terms with equal power of h 1 andh 1 become relevant as they are not suppressed when N c → ∞. However, these contributions are of mesonic nature and expected to scale as ∼ N 0 c . This can be shown explicitly for N f = 1 at leading order in the hopping expansion, as the contribution of the static determinant is known to be [49] with T (k) = min(k,2Nc−k)+3 3 and P (k) = (N c + 1 − k)(k + 1). For µ < a ln(2κ) (i. e. h 1 < 1) the µ-dependent terms vanish as N c → ∞. Obviously, this is not the case for the µindependent terms. Their contribution can be evaluated for N c → ∞ by using the following variant of the geometric series (which can be obtained by differentiation): This results for N c → ∞ in the pressure

Gauge corrections and 'tHooft scaling
Our discussion so far has been for the strong coupling limit β = 0, and so does not correspond to the intended 't Hooft limit of large N c . In this section we argue that our observations carry over to the t'Hooft limit once we include gauge corrections, at least to the orders considered.
Including the leading gauge corrections to the fermion determinant, as well as partial resummation of the corresponding diagrams to all orders, proceeds in the same manner as in [27] for any N c . Through O(κ 2 ) the free energy density now reads where u can be computed using (3.2) and h 1 now includes corrections, In taking the 't Hooft limit, i. e. keeping λ H = 2N 2 c /β fixed, one has for λ H > 1 [42,44] Therefore, these gauge corrections only modify the asymptotic behaviour of the thermodynamic functions by a constant ∼ N 0 c . Furthermore, we also consider the leading order contribution from the pure gauge sector to the effective theory (the first line in equation (2.22)) with (3.38) to leading order in the character expansion. The first correction to −f in equation (3.6) due to the inclusion of this term, combined with the centre-symmetric part of the static determinant, reads Employing the previously outlined integration techniques one obtains  Figure 4: Baryon density, including gauge corrections, for growing N c with the 't Hooft coupling held fixed. The qualitative behaviour is the same as in the strong coupling limit.
In the 't Hooft limit λ 1 = 1 , and therefore the asymptotic analysis of this term can be done in the same way as for the strong coupling contributions. The result is dW Tr(W ) n Tr(W † ) n = n! for n ≤ N c (3.43) are relevant. When taking the large N c limit, order by order these contributions to the pressure are µ-independent and scale as ∼ λ k 1 ∼ N 0 c . The last statement hinges on the fact that the N c dependence of λ 1 is solely determined by u. In [41] corrections to λ 1 to O(u 8 ) were computed and, although some corrections do introduce additional N c factors, those cancel order by order when all corrections are summed up. A similar observation, including higher representations, was made in [22] in the context of a strong coupling expansion of pure gauge theory without using an effective theory.
The influence of the gauge corrections is illustrated in Fig. 4 for two different choices of the 't Hooft coupling. Clearly, the small quantitative modifications by the gauge corrections do not alter the qualitative N c -behaviour observed earlier. Of course, in higher orders the situation might be more complicated, as new interactions can arise in the effective theory with non-trivial N c -dependence. Nevertheless, the dominant contribution to the large N c limit of baryon dynamics should always be represented by powers of h Nc like in the leading contributions considered here, since N c occurs in the exponent in these cases.

Approaching the continuum
This situation leaves, however, one caveat. Even if one would be able to include gauge contributions to all orders, the interchange of the N c → ∞ limit and the strong coupling expansion was observed to be "highly suspicious" in the case of QCD in 1+1 dimensions [44]. Our analysis so far was based on taking N c large before a continuum limit. The fact that the density at the onset transition immediately jumps to lattice saturation indeed suggests that the limits should be taken in the opposite order, if one is interested in continuum physics. In this case, the interplay between large N c and the Pauli principle should lead to a finite continuum density, just as for N c = 3.
To get an idea if our results are consistent with this expectation we investigated the behaviour of the baryon density towards the continuum. To set the scale at SU (3) we use the same strategy as in [27], which gives a rough estimate of the parameter space. At first, since heavy quarks have little influence on the running of the coupling, we use the non-perturbative beta-function of pure gauge theory to get a relation between β SU (3) and a/r 0 , where r 0 is the Sommer parameter [50]. Using r 0 = 0.5 fm sets a physical scale for the lattices and the temperature can be adjusted by N τ via T = 1/(aN τ ). To obtain the corresponding β for SU (N c ) we keep λ H = 2N 2 c /β SU (Nc) = 18/β SU (3) fixed. Finally, the leading order expression (2.20) is used to keep the constituent quark mass constant for different a.
The outcome of this procedure is illustrated in Fig. 5, where the lattice spacing is varied for fixed N c = 3, 9. Each continuum extrapolation leads to a finite value for the density, which is smaller or larger, for h 1 < 1, h 1 > 1, respectively, as N c is increased. It is thus apparent that the transition steepens with growing N c also when the continuum is approached first. Similarly, the pressure is shown in Fig. 6 for two different lattice spacings and N c ∈ [3,9]. Since this is far from the large N c limit, we explicitly checked that the first subleading contribution goes as ∼ N 0 c . The figure shows that the full result follows the predicted N c -scaling to hold in a region where the lattice is only about half filled, i.e., not yet dominated by lattice saturation. This behaviour appears to be stable as the lattice is made finer. While for N c = 3 and sufficiently heavy quarks a continuum limit can be explicitly taken [27,29], for large N c this is difficult in practice, because in the problematic transition region, cf. section 3.5, the required length of the hopping series grows exponentially with N c . We are therefore unable to explicitly demonstrate first-order behaviour or the scaling p ∼ N c of the large N c -limit in the continuum.

The phase diagram with growing N c
We have seen that, in the strong coupling limit, the onset transition to baryon matter becomes first-order for any temperature. In the last section we provided evidence for a steepening of the liquid gas transition with N c also in the continuum. While we are unable to take the continuum and large N c limits in this order, we now argue on physical grounds that the onset transition has to become first order when N c → ∞. Within the effective theory (as well as physical QCD), the endpoint of the nuclear liquid gas transition is located where the temperature starts to exceed the binding energy per baryon. In table 1, we found the interaction energy in units of baryon mass to scale as ∼ const. Consequently, the binding energy in N c -independent units scales as ∼ N c as expected, and the critical endpoint of the liquid gas transition moves to ever larger temperatures. On the other hand, the deconfinement transition temperature T d is only sensitive to meson physics, and in the large N c -limit is within ∼ N −2 c of its value at N c = 3 [51]. Hence, in the limit N c → ∞, temperatures in the range 0 < T < T d never exceed the binding energy between baryons and the onset transition must be of first order.
This implies that the critical endpoint of the baryon onset transition increases with growing N c until it hits another discontinuity, as indicated in Fig. 3 (right). We know already from perturbation theory that also the deconfinement transition line "straightens out" with growing N c , as the deconfinement transition becomes less and less sensitive to the quark contributions. Altogether, we then observe how the predicted rectangular phase diagram of Fig. 2 emerges continuously by increasing N c , as indicated in Fig. 7 .

Quarkyonic matter on the lattice?
While lattice saturation is a mere discretisation artefact and may seem uninteresting from a continuum perspective, it adds an intriguing feature here. The lattice saturation density is clearly determined by the quark degrees of freedom. Besides counting degrees of freedom, this follows also from the fact that the saturation density is precisely the same in the case of large isospin chemical potential [27]. Thus, approaching saturation, the lattice is filled with quark matter. On the other hand, the onset transitions at finite baryon as well as isospin chemical potentials are related to the condensation of hadrons, and not quarks. This follows from the fact that the critical chemical potential is different in the baryonic and isospin cases [27], i.e. m B /3 = m π /2.
For increasing baryon chemical potential, a lattice filling up with baryon number is thus consistent with the picture of quarkyonic matter, in the sense that it shows a smooth transition from baryon matter to quark matter in a remarkably narrow range of chemical potentials. As the lattice is made finer, the saturation level in physical units increases and is reached at larger chemical potentials. Eventually, in the continuum limit, the interplay between the attractive baryon interaction and Pauli repulsion will lead to the physical saturation density known from nuclear matter, while the quark matter observed at lattice saturation gets shifted to larger chemical potentials (possibly infinite), with a quarkyonic regime in between.

Conclusions
We have studied the large N c -behaviour of QCD in the cold and dense regime within an effective lattice theory derived by combined strong coupling and hopping expansions, which is valid for sufficiently heavy quarks. At low temperatures and µ B ∼ m B it exhibits a transition to baryon condensation, which is the heavy quark analogue of the nuclear liquid gas transition. By considering the effective theory for general gauge group SU (N c ) we have shown that, in the strong coupling limit and through three consecutive orders in the hopping expansion, the pressure in the baryon condensed phase scales as p ∼ N c and the onset transition becomes first-order for large N c . This behaviour is stable under inclusion of the leading gauge corrections. We have pointed out that the continuum limit has to be taken before the large N c limit, which makes definite conclusions for continuum physics much harder to reach. Nevertheless, we have shown our findings to be stable also in this ordering for the range of lattice spacings and N c -values we were able to consider with our truncated series. Our results are thus consistent with the large N c phase diagram and the definition of quarkyonic matter proposed in [9].
For N c = 3, onset to baryon matter happens at µ c B < m B (because of the binding energy between baryons) whereas h 1 = 1 at some µ B > m B . The onset transition then marks the condensation of baryons, which smoothly turn into quarkyonic matter, whose pressure scales as p ∼ N c and whose effective degrees of freedom can in principle change smoothly from baryon-like to quark-like as a function of chemical potential.
The QCD parameter values realised by nature are m u,d ∼ 2 − 5 MeV, in contrast to the heavy quarks on which our analysis above is based. What can be said about the physical situation? As remarked in section 3, the large N c analysis is independent of the current quark masses, with m B ∼ N c always and Feynman diagrams with quark loops suppressed. Thus, whether there is a deconfinement transition, a chiral transition or a crossover at N c = 3 is immaterial for the forming of a first-order horizontal deconfinement line at large N c . On the other hand, the baryon onset transition is also present in a different effective lattice theory derived from QCD, which is valid in the chiral limit and the strong coupling regime [52]. (Moreover, it is also expected from nuclear physics [53].) If our results generalise to all orders in the hopping expansion and are stable in the proper order of limits, then we expect the N c -scaling of the pressure as well as the evolution of the baryon onset transition with growing N c to look qualitatively the same when starting from the physical situation, i.e. these would be genuine features of SU (N c )-QCD. Whether or not there is also a chiral transition at some larger chemical potential cannot be decided within the current framework, but requires additional investigations in an effective theory including chiral symmetry, such as [52]. matter under extreme conditions" and by the Helmholtz International Center for FAIR within the LOEWE program of the State of Hesse.