Dense two-color QCD towards continuum and chiral limits

Tamer Boz , Pietro Giudice , Simon Hands , and Jon-Ivar Skullerud 1,5,*,† Department of Theoretical Physics, National University of Ireland Maynooth, Maynooth, County Kildare, Ireland Istituto Arcivescovile Paritario Santa Caterina, Piazza Santa Caterina 4, 56127 Pisa, Italy INFN Sezione di Pisa, Largo Pontecorvo 3, 56127 Pisa, Italy Department of Physics, College of Science, Swansea University, Singleton Park, Swansea SA2 8PP, United Kingdom School of Mathematics, Trinity College, Dublin 2, Ireland


I. INTRODUCTION
The structure of strongly interacting matter at high densities and low to moderate temperatures remains an outstanding problem, with applications to compact stars, neutron star mergers, and the next generation of heavy-ion colliders at FAIR and NICA. First-principles studies of this regime are hindered by the sign problem: with chemical potential μ ≠ 0 the Euclidean action becomes complex, and can therefore not be used as a probability weight in Monte Carlo simulations, which are the mainstay of lattice gauge theory, the method of choice for first-principles, nonperturbative quantum field theory. Despite recent progress in alternative sampling approaches such as the density of states method [1], complex Langevin [2] and Lefschetz thimble and related approaches [3,4], we do not as yet have any method that has been shown to yield valid and reliable results for real QCD.
The problem may be circumvented by studying QCDlike theories without a sign problem, such as theories with adjoint fermions in any gauge group, QCD with isospin chemical potential [5], or QCD with gauge groups SU(2) (QC 2 D) [6][7][8] or G 2 (G 2 -QCD) [9,10]. Although these theories all have qualitative features which distinguish them from real QCD at nonzero baryon chemical potentialnotably a gauge-invariant Bose-Einstein condensate (BEC) above an onset chemical potential μ o -they share salient features such as spontaneous chiral symmetry breaking and confinement at T ¼ μ ¼ 0, and may be used as laboratories for strongly interacting theories at high density. Lattice results from these theories may also be used as a check on the approximations made in other approaches which do not suffer from the sign problem, including Polyakov-loop extended Nambu-Jona-Lasinio models [11,12], massive perturbation theory [13,14], quark-meson(-diquark) coupling models [15], the functional renormalization group [16] or Dyson-Schwinger equations [17,18].
In a previous series of papers [6][7][8]19], we have studied the phase structure of QC 2 D with N f ¼ 2 Wilson fermions. The main findings of these studies have been that at high density and low temperature, there is a "quarkyonic" phase [6] where the diquark condensate and quark number density scale with the quark chemical potential μ in the same way as in a system composed of noninteracting fermions disrupted by a Bardeen-Cooper-Schrieffer (BCS) gap. The diquark condensate, signaling superfluidity, vanishes at a critical temperature which appears to be approximately independent of μ above the onset chemical potential μ o ¼ m π =2 [8]. At high temperature, there is a transition to a deconfined quarkgluon plasma, with the pseudocritical temperature T d decreasing with increasing μ [8]. It is as yet unclear whether T d goes to zero at any finite μ and hence whether there is any deconfinement transition at high density and low temperature; early indications of such a transition [19] may have been complicated by lattice artefacts.
These studies have all been carried out with quite heavy quarks (m π =m ρ ¼ 0.8) and on fairly coarse lattices (a ¼ 0.18-0.23 fm). The aim of the current paper is firstly to gain control over lattice artefacts by reducing the lattice spacing at fixed m π =m ρ , and secondly to explore the quark mass dependence by studying a system with lighter quarks at fixed lattice spacing. The latter is of particular significance as it might point to a "BEC region" which can be described using chiral perturbation theory (χPT) incorporating both mesonic and baryonic (diquark) Goldstone degrees of freedom [20].
There have been a number of other lattice studies of dense QC 2 D in recent years, using staggered [21][22][23] and Wilson [24] fermions. These have to a large extent confirmed the picture outlined above, with some additions. Notably, in [21], with a smaller pion mass than in [6][7][8], a BEC region was found where the diquark condensate and quark number density agree with predictions from chiral perturbation theory, followed by a transition to a quarkyonic region at higher μ. Also, the chiral condensate was found to vanish in the chiral limit in both the BEC and quarkyonic regions. The system was found to be confined at low temperature, with deconfinement only setting in at much larger chemical potential [22]. Similar conclusions were found in [24].
It is also worth noting that a recent study of QCD with nonzero isospin chemical potential [5] found a phase diagram very similar to that of [8], namely a pion condensed phase at low T and large μ, with a critical temperature that is nearly independent of μ, and a deconfinement transition line that intersects with the pion condensation transition.
The structure of this paper is as follows. In Sec. II we describe our simulation parameters and determination of the lattice spacing for our new (fine) ensemble. Results from this ensemble are presented in Sec. III. First, in Sec. III A we present results for the superfluid order parameter, the diquark condensate, including its scaling with chemical potential and estimates for the critical temperature. Section III B contains our results for the Polyakov loop and deconfinement transition, while Sec. III C contains results for the quark number density. Section IV contains our results from simulations with lighter quarks. We summarize our findings in Sec. V.

II. SIMULATION DETAILS AND SCALE SETTING
We study QC 2 D with a conventional Wilson action for the gauge fields and two flavors of Wilson fermion. The fermion action is augmented by a gauge-and isosinglet diquark source term which serves the dual purpose of lifting the low-lying eigenvalues of the Dirac operator and allowing a controlled study of diquark condensation. The quark action is where i ¼ 1, 2 is a flavor index and Further details about the action and the Hybrid Monte Carlo algorithm used can be found in [19].
We have studied three ensembles, which in the following we call "coarse," "fine," and "light." The parameters are shown in Table I, together with the values obtained for the pion (pseudoscalar meson) mass m π , ratio of pion to rho (vector meson) mass m π =m ρ and lattice spacing a. The coarse ensemble is the same as was used in [6][7][8]. The parameters for the fine ensemble were chosen to give the same value of m π =m ρ as the coarse ensemble, while those of the light ensemble were chosen to give approximately the same lattice spacing as the coarse ensemble, but with a smaller value of m π =m ρ ≈ 0. 6. Further details about the coarse and light ensemble parameters can be found in [7].
To determine the lattice spacing, we extracted the static quark potential VðrÞ from rectangular Wilson loops Wðr; τÞ by fitting Wðr; τÞ ¼ expð−VðrÞτÞ for τ=a ¼ T min ; N τ − 1.  The lattice spacing was then determined by fitting the static quark potential to the Cornell form and taking the string tension to be ffiffiffi σ p ¼ 440 MeV. The static quark potential for the fine ensemble is shown in Fig. 1 and the best fit values are given in Table II.
We have performed scans in μ at fixed temperature for four temperatures on the fine ensemble, and three temperatures on the light ensemble. The lattice volumes, temperatures and range of chemical potentials are given in Table III, along with those from the coarse ensemble that are used for comparison. For most of these temperatures and chemical potentials we have produced configurations at two values for the diquark source, namely ja ¼ 0.02, 0.03 on the fine ensemble, and ja ¼ 0.02, 0.04 on the light ensemble. A third diquark source value (ja ¼ 0.01 for the fine ensemble, ja ¼ 0.03 for the light ensemble) has been added for selected parameter values to control the systematics of the j → 0 extrapolations. For a few parameters, mostly at the highest and lowest μ-values, only a single j-value has been used; in these cases, only results for the Polyakov loop, which is only weakly dependent on j, will be shown here. In addition, we have performed temperature scans at fixed values of μa ¼ 0, 0.2, 0.3, 0.4, 0.5 on the fine ensemble, with N τ ¼ 18-4 corresponding to T ¼ 63-356 MeV.
We note that by assigning a temperature T > 0 to our ensembles with N τ ≥ N s we deviate from what is commonly done in lattice thermodynamics studies. This is appropriate in the presence of a chemical potential (which modifies the temporal boundary conditions). At weak coupling and low temperature, relativistic quarks form a Fermi surface with Fermi momentum k F ≃ μ q . The resulting ground state is highly degenerate, so that unlike the case when μ ¼ 0 the discrete momenta on the finite volume do not lead to a large energy gap. For this reason in our view one should always consider T as nonzero when studying systems with μ q ≠ 0, even if N s < N τ .
For enhanced separation of scales m π ≪ m ρ , the behavior as μ increases at zero temperature may be analyzed using χPT [20], in this context an effective theory of tightly-bound qq mesons and qq;qq baryons. For μ ¼ 0 the chiral symmetry of continuum QC 2 D is spontaneously broken from SUð2N f Þ to Spð2N f Þ yielding N f ð2N f − 1Þ − 1 Goldstones; for N f ¼ 2 these are the pseudoscalar pion triplet and a scalar diquark/antidiquark pair. As μ is raised there is a second-order onset transition at μ o ¼ m π =2 to a phase where the baryon density n q and superfluid diquark condensate hqqi are both nonzero: Here hqqi 0 is the chiral condensate at μ ¼ 0 and F π the χPT parameter known as the pion decay constant. Note that the onset transition in physical QCD is first order.

III. RESULTS FROM FINE ENSEMBLE
A. Diquark condensation Figure 2 shows the diquark condensate, divided by the square of the chemical potential, as a function of chemical potential, for all temperatures. Physically, hqqi is the density of relativistic quark pairs contributing to the superfluid condensate. In the case of a weakly coupled BCS condensate at the Fermi surface, this should be roughly equal to the momentum-space volume of   (4) predicts hqqi to be μ-independent for μ ≫ μ o and hence the quantity plotted should fall off like μ 2 . For the two lowest temperatures, we find that hqqi=μ 2 is almost independent of both μ and T for μ ≳ μ o . This agrees with what was found previously for the coarse ensemble [7]. For μ < μ o the diquark condensate rises gradually from zero; we take this to be primarily an artefact of the linear extrapolation in j.
To further investigate the diquark source dependence, in Fig. 3 we show hqqi as function of j for different values of the chemical potential. With our current data, we see no evidence of deviation from a linear form at any μ. Specifically, the form hqqi ¼ Aj 1=3 , which should hold exactly at μ ¼ μ o , does not fit the data for any of our μvalues, even with the addition of a constant term. Fits to a more general power-law form, hqqi ¼ A þ Bj α , yield powers ranging from α ∼ 0.9 for μ < μ o to α ∼ 0.6 at μa ¼ 0.4. The apparent linear j-dependence at all μ may be an indication that nonanalytic behavior only sets in at lower j-values than we have here. To improve on the diquark source extrapolation, it may be necessary to carry out a reweighting procedure as outlined in [5], or explore much smaller values of the source ja [21].
At T ¼ 90 MeV (N τ ¼ 16) we see that hqqi is significantly smaller for all values of μ, suggesting that this temperature is near the critical temperature for the superfluid to normal phase transition. This is again in agreement with what was found for the coarse lattice in [7,8]. Finally, at the highest temperature, T ¼ 120 MeV (N τ ¼ 12), the diquark condensate is consistent with zero, suggesting that at this temperature we are in the normal phase.
To study the superfluid to normal phase transition in more detail, we have performed temperature scans at fixed values of chemical potential aμ ¼ 0.3, 0.4, 0.5 and aj ¼ 0.02, 0.03. These are fixed-scale temperature scans; i.e., the temperature is varied by changing N τ , without changing the lattice spacing. This means that although our data for hqqi are not renormalized, this will just contribute an overall factor without changing the shape of the curves.
The results of these scans are shown in Fig. 4. From the data extrapolated to j ¼ 0 using a linear Ansatz, we see evidence of a phase transition at T s ≈ 110 MeV, independent of the chemical potential. In order to determine this transition in a more controlled manner, we find the inflection points for ja ¼ 0.03 and 0.02 using a cubic spline interpolation, and extrapolate these to j ¼ 0. The results are shown in Table IV. This yields a somewhat lower temperature T s ¼ 90-100 MeV, which is consistent with the result quoted in [8], T s ¼ 93ð8Þ MeV. We note that the data are for a single volume and we can therefore not determine the order of the transition, but it is expected to be a second order transition (in the O(2) universality class), and the data are consistent with this.
B. Deconfinement Figure 5 shows the order parameter for deconfinement, the Polyakov loop hLi, for our four different temperatures. The renormalized Polyakov loop L R is given in terms of the bare Polyakov loop L 0 and the temperature T ¼ 1=ðaN τ Þ by Just as in [8], we use two different renormalization schemes, The results in Fig. 5 have been obtained using Scheme A. We see no evidence of any deconfinement transition at our lowest temperature, T ¼ 45 MeV, corresponding to N τ ¼ 32. This suggests that the high-density deconfinement transition found on coarser lattices at temperatures of 40-50 MeV [7,19] is primarily a lattice artefact. At higher temperatures we see that hLi increases rapidly from zero above a chemical potential μ d ðTÞ which is a decreasing function of temperature, in agreement with previous results.
To determine the transition line, we study the variation of the renormalized Polyakov loop with temperature at fixed values of the chemical potential μa ¼ 0.   increased. However, this effect is not large: the inflection point in Scheme B actually appears to increase between μ ¼ 0 and μa ¼ 0.2, and even at our highest chemical potential, μa ¼ 0.5, the inflection point in Scheme B is still consistent with the μ ¼ 0 inflection point in Scheme A. As can be seen from Fig. 7, the transition in Scheme A happens at a considerably lower temperature than in Scheme B, which may be taken as an indication of the width of the crossover and associated uncertainty in the transition temperature. It would be useful to compare these scheme-dependent results with other quantities sensitive to deconfinement, such as the static quark potential [22] or the entropy of a static quark [25].
Our results for the deconfinement transition temperature are shown in Table V and Fig. 8. We find a large discrepancy between the two schemes, and this discrepancy increases with μ; however, the uncertainties are large, in particular in scheme A which has a lower T d . There are indications that, rather than decreasing monotonically toward 0 at large μ, T d approaches a constant value which is above (scheme B) or close to (scheme A) the superfluid transition temperature T s ≈ 90 MeV. If this is confirmed, it means there is no deconfining transition at T ¼ 0, and that the superfluid phase remains confined at all μ.

C. Quark number density
The final quantity to consider is the quark number density n q , which forms the basis for computing bulk thermodynamic properties such as the pressure. In Fig. 9 we show the quark number density extrapolated to zero diquark source (using a linear extrapolation) at the lowest temperature studied, corresponding to N τ ¼ 32. As was the case for the diquark condensate, we do not find significant evidence of a deviation from a linear behavior for the range of j-values we have. The data are plotted in dimensionless form by normalizing by the density n SB for a gas of massless noninteracting fermions. Figure 9 shows results obtained using two choices for n SB : first using the continuum result and second, using a form n lat SB evaluated as a mode sum using the action (1), (2) with κ ¼ 0.125 and U μ ≡ 1 on a finite lattice as described in [19]. Both forms were explored in [7], with the conclusion that while n lat SB offers a better correction for UV artifacts at large μa, it is prone to significant IR artifacts for μ ≳ μ o , so that n cont SB is preferred   (4) 210 (17) 214 (9) in this regime. In the current work we have addressed this issue by evaluating n lat SB on a much larger spatial volume than that used for simulation, specifically 96 3 × N τ (with 128 3 × N τ checked to see that this was sufficient). This removes the IR artifacts but retains the UV corrections.
We find that n q =n SB rises rapidly above the onset transition, then descends to reach a plateau for μ ≳ 500 MeV. This is in qualitative agreement with the predictions of χPT shown in Fig. 3 of [19]. Within the fairly large errors resulting from the j → 0 extrapolation the results obtained in the neighborhood of onset with continuum and lattice free fermions are compatible, but for μ ≳ 700 MeV the n cont SB curve continues to rise while the lattice normalization yields a plateau with n q =n SB ≲ 1; for the reasons given in the previous paragraph this is the normalization we prefer and will use henceforth. Also shown in Fig. 9 is the ratio evaluated using free massive Wilson fermions with κ ¼ 0.120; in this case the value of n q =n SB on the plateau is consistent with unity. Due to the difficulties inherent in assigning a bare quark mass for interacting Wilson fermions, we therefore do not draw any physical conclusions from the plateau height at this stage.
An important observation is that there is no longer evidence for a regime at high μ where n q =n SB increases above its value on the plateau (cf. Fig. 11 of [7]; were the high-μ behavior in that plot physical, then a corresponding rise would be expected on the fine lattice at μ ≈ 700 MeV). In conclusion, at the lowest temperature studied the equation of state looks "quarkyonic", i.e., with n q ðμÞ ≈ n SB , all the way along the μ-axis, with no evidence of a qualitative change associated with deconfinement. This is consistent with our results for the Polyakov loop, which show no sign of a deconfinement transition at low temperature, and the high-μ increase seen in [7] is therefore most likely a lattice artefact.
As shown in Fig. 10, this behavior persists for the lowest three temperatures, which according to the results in Sec. III A are all in the superfluid region (or near the transition temperature, for N τ ¼ 16). The data suggest the plateau value of n q =n SB falls with increasing T, although uncertainties following the j → 0 extrapolation are significant. At the highest temperature, T ¼ 119 MeV (N τ ¼ 12), which falls in the quark-gluon plasma phase, we see indications of n q =n SB monotonically increasing with μ, in qualitative agreement with the findings of [7].  To make this comparison quantitative, we show the results from the fine and coarse ensemble together, in units of the onset chemical potential μ o ¼ m π =2 (allowing also for a comparison of results for the different quark masses), in figure 11. At low temperature (upper panel) there is quantitative agreement between the two ensembles for μ ≲ 2μ o . For larger μ the rise in n q =n SB seen for the coarse ensemble (which might have signaled a transition to a different state of matter) is absent for the fine ensemble, and instead we see that n q remains close to n SB throughout. This is consistent with our results for the Polyakov loop, which show no sign of a deconfinement transition at low temperature, and the high-μ increase seen in [7] is therefore most likely a lattice artefact. The same pattern is repeated at high temperature (lower panel), where now we see n q approach n SB from below. However, part of the difference between the two ensembles may in this case be due to the slightly different temperatures (120 vs 130 MeV for the fine and coarse ensembles respectively).
Next we discuss the pressure, which as outlined in [7,19] can be obtained in the limit of low T via an integral of the form p ¼ R μ n q ðμ 0 Þdμ 0 . In order to arrive at the dimensionless ratio p=p SB starting from our data there are several different quadrature schemes available: here we use introduced as "scheme II" in [7] (μ 0 is the lowest available value in the dataset). Reassuringly, with n lat SB now defined on a large spatial volume to eliminate IR artifacts we find results compatible with those computed using scheme I, which was not the case in [7]. The results for data extrapolated to j → 0, at low temperatures where there are enough μ-points to control the numerical integration, are shown in Fig. 12.
The main result is that p=p SB increases sharply after onset, reaching a plateau at μ=μ o ≈ 2. It appears to approach the plateau from below, in contrast to the χPT prediction that the SB limit is approached from above [19]. While, as in the discussion of n q =n SB , it is premature to assign a precise value to the height of the plateau, it is worth recalling that e.g., in the Van der Waals equation of state the ideal gas pressure receives a downward correction due to attractive forces between particles.

IV. RESULTS FROM LIGHT ENSEMBLE
We now study the effect that reducing the quark mass may have on the phase structure and equation of state. While our parameters are still very far from the chiral limit, this may give us an idea of which, if any, qualitative changes may occur as we approach this limit. Figure 13 shows the diquark condensate from the light ensemble for our three different temperatures, extrapolated  to j ¼ 0. On the face of it, these results are qualitatively different from our results with heavier quarks, in that the diquark condensate no longer scales like μ 2 in the region just above the onset transition. This might be taken as an indication that a BEC window is opening up where, rather than following the BCS scaling hqqi ∝ μ 2 , the condensate behaves according to the predictions (4) of zero-temperature chiral perturbation theory.
To facilitate this comparison, we have included in Fig. 13 the χPT curve (4), with an arbitrary prefactor, and we see that the data for μ o < μ < 400 MeV make contact with this curve.
This interpretation is, however, complicated by our lack of control of the j → 0 extrapolation. For nearly all points, we only have two values of the diquark source (ja ¼ 0.02, 0.04) available and have used a simple linear extrapolation. For μa ¼ 0.2, 0.25 and 0.3 on the 12 3 × 24 lattice we also have data for ja ¼ 0.03, and have also used a quadratic extrapolation as well as a power lawþ constant form, hqqiðjÞ ¼ A þ Bj α . While the quadratic extrapolation gives results roughly consistent with the linear form, the power law form gives a result which is consistent with BCS scaling at all μ. It should be noted that near the onset transition, we expect a power-law scaling with α ¼ 1=3.
We also find that hqqi ≠ 0 also for μ < μ o , and we take this to be an indication that the linear extrapolation breaks down in this regime. Figure 14 shows the diquark condensate as a function of j for those chemical potentials where we have 3 j-values at our disposal. We see that the data are consistent with a linear behavior as was also found for the fine ensemble, but may also be described with a power-law þ constant form. A better control over the diquark source extrapolation, for example along the lines of [5], is required to determine whether there is indeed a BEC window for these parameters.
The hqqi results for T ¼ 87 MeVðN τ ¼ 12Þ are almost identical to those for T ¼ 43 MeV, except for μ ≳ 600 MeV. However, at the highest temperature, T ¼ 130 MeV (N τ ¼ 8) we see that the diquark condensate is much smaller, and vanishes at high μ. It is worth noting, however, that hqqi does not vanish at all μ, as was the case at comparable temperatures for the coarse and fine ensembles. This may be taken as a first indication of a μ-dependent superfluid to normal transition temperature. A controlled j → 0 extrapolation combined with temperature scans at fixed μ would be required to draw a firm conclusion.
Looking now at the deconfinement transition, Fig. 15 shows the Polyakov loop as function of chemical potential for our three temperatures. In the absence of a μ ¼ 0 temperature scan to fix the renormalization constant, we show the unrenormalized Polyakov loop. This does not affect the shape of each of the curves, only the relative magnitude of the data at different temperatures. We see the same picture as before, where hLi ≃ 0 at low temperature (except for μ ≳ 800 MeV, corresponding to μa ≳ 0.8, which in the light of previous results we attribute to a lattice artefact), but increasing with μ at higher temperature. At T ¼ 130 MeV (N τ ¼ 8), we note that hLi ≠ 0 also at μ ¼ 0, suggesting that this temperature is near the transition region for this ensemble.
In Fig. 16 we show the quark number density n q =n SB calculated using the same procedure outlined in Sec. III C. At low temperature, this has a plateau for μ > μ o , in line with our previous results, while at high temperature we see a steeper increase with μ, as seen before, which is much more pronounced than that for the fine lattice shown in Fig. 10. On a coarse lattice the impact of UVartefacts in this regime cannot be excluded. However, as can be seen in the lower panel of Fig. 16, the quark number density at the lowest temperature is also consistent with the χPT prediction (5), as indicated by the j → 0 extrapolated curve tending to a nonzero constant as μ → μ oþ . It therefore remains an open question whether we for these parameters have a window characterized by a Bose-Einstein condensate of tightly bound diquarks.
In Fig. 11 we compare our results for the quark number density with those obtained on the coarse ensemble in [7]. We see that the smaller quark mass has a large quantitative effect, giving a much lower number density. Because μ o a is now smaller, the increase in n=n SB at large μ, which is likely to be a lattice artefact, now only appears at larger μ=μ o . At high temperature, however, we do not see any plateau in n q =n SB , as was the case for the fine ensemble, but instead a monotonic increase with μ.
Finally, we note that results for p=p SB for the light ensemble are included in Fig. 12. The notable feature is that the rise from zero following onset is now much steeper, and an approximate plateau established by μ=μ o ≲ 1.5. As before, we are reluctant to overinterpret the disparity in plateau height between fine and light ensembles.

V. CONCLUSIONS AND OUTLOOK
We have performed lattice simulations of two-color QCD at different lattice spacings and different quark masses, as a step toward the continuum and chiral (or light-quark) limits of this theory. Our investigation of the phase structure of the theory in the ðT; μÞ plane has yielded results in qualitative agreement with earlier studies [6][7][8]; most notably it has confirmed the quarkyonic phase at low T and intermediate to high μ. In this phase, quarks are confined but the bulk thermodynamics as well as the diquark condensate behave as if the system consists of a Fermi sphere of weakly interacting quarks. This phase is found to extend to larger μ as the lattice spacing or the quark mass is reduced. A striking result is that at the lowest temperatures explored the quark number density is found to be very close to the Stefan-Boltzmann value. Because we have not yet determined the physical quark mass, it is not possible to state with confidence whether the ratio n q =n SB is greater or less than one in this regime. However, Fig. 11 does suggest that the ratio falls as the quark mass decreases. The superfluid phase transition temperature is found to be independent of μ at least for our heavier quark mass, corresponding to m π =m ρ ¼ 0.8. The transition temperature T s ≈ 90 MeV is in quantitative agreement with that found on a coarser lattice in [8]. This qualitative behavior also agrees with what has been found for QCD at high isospin density [5]. For lighter quarks, however, there are indications that T s may increase with μ.
Although we find a deconfinement transition at high temperature, the deconfinement transition previously seen at low T and high μ appears to be a lattice artefact-at T ∼ 45 MeV the Polyakov loop appears to increase from zero at μa ∼ 0.8 for all lattice parameters. Simulations using staggered fermions at a relatively high temperature T > 100 MeV [21,22] have found a deconfinement transition at μ ≈ 750 MeV characterized by the string tension dropping to zero. This is not in contradiction with the absence of any deconfinement transition at T ≈ 45 MeV.
The location of the high-temperature deconfinement transition is not yet clear. Our results suggest a broad crossover with a transition temperature T d that decreases with increasing μ. The analysis is complicated by lack of precise data for the renormalized Polyakov loop at low temperature, as well as a small but nontrivial diquark source dependence. There are indications that T d may tend to a constant at large μ and that T d ≳ T s for all μ.
There remain significant uncertainties relating to the extrapolation to zero diquark source j, with our data unable to distinguish between a linear j-dependence and the nonanalytic behavior expected at least in the vicinity of the onset transition. This may be mitigated by adopting the reweighting method introduced in [5] and implemented in the context of QC 2 D in [24]; this will be left for future investigations.
Another source of uncertainty is the use of relatively small volumes, with N s < N τ for our lowest temperatures.
Additional simulations on larger volumes would be required to reliably determine if our assignment of a nonzero temperature to these ensembles is correct. We note that finite volume effects were studied previously in Ref. [7]; there it was found that the results for the various physical quantities, and in particular the Polyakov loop, agreed between the two volumes considered, suggesting that finite volume effects are indeed small.
A reliable continuum extrapolation would require significantly finer lattices than those used in this study. The first step toward further reducing the lattice spacing errors would be to employ an improved fermion action. This, together with algorithmic improvements, would also allow us to study lighter quark masses which are out of bounds with the unimproved Wilson action used here.