Lattice study of static quark-antiquark interactions in dense quark matter

In this paper we study the interactions among a static quark-antiquark pair in the presence of dense two-color quark matter with lattice simulation. To this end we compute Polyakov line correlation functions and determine the renormalized color averaged, color singlet and color triplet grand potentials. The color singlet grand potential allows us to elucidate the number of quarks induced by a static quark antiquark source, as well as the internal energy of such a pair in dense quark matter. We furthermore determine the screening length, which in the confinement phase is synonymous with the string breaking distance. The screening length is a decreasing function of baryon density, due to the possibility to break the interquark string via a scalar diquark condensate at high density. We also study the large distance properties of the color singlet grand potential i a dense medium and find that it is well described by a simple Debye screening formula, parameterized by a Debye mass and an effective coupling constant. The latter is of order of unity, i.e. even at large density two-color quark matter is a strongly correlated system.


Introduction
Knowledge of the properties of QCD at large baryon density is needed to interpret the results of heavy ion collisions experiments. In particular, this is the case at the future experiments of NICA (JINR, Dubna) and FAIR (Darmstadt, Germany), which are designed to study the region of high baryon density. Input from the theory side is hence urgently needed. An understanding of the properties of matter in the corresponding region of the QCD phase diagram is also extremely important in astrophysics, for example, for a correct description of the fusion of neutron stars.
In general, lattice QCD is a powerful tool to study the non-perturbative properties of strongly interacting matter from first principles. By virtue of lattice simulations vital insight into QCD at finite temperature [1], nonzero magnetic field [2], isospin chemical potential [3,4] and chiral chemical potential [5][6][7] has been obtained.
An interesting area of finite temperature lattice simulations is the study of the interaction between a quark-antiquark pair and the interaction of the pair with the QCD medium (see e.g. [8][9][10][11][12][13]).
The in-medium properties of QCD are prominently encoded in the correlation function of Polyakov loops (i.e. Polyakov lines with maximum temporal extent). As the Polyakov JHEP05(2019)171 loop constitutes a quasi order parameter of the strong interactions, its correlator is greatly affected by the phase transitions which take place in QCD. One pertinent example is the confinement/deconfinement phase transition which considerably modifies the value of the correlator.
Furthermore the Polyakov loop correlator is directly related to the free energy of the in-medium quark-antiquark pair. In the confinement phase the free energy extracted is known to be a linear increasing function at intermediate distances. At zero temperature and distances ∼ 1.2 fm the free energy asymptotes to a plateau due to the string breaking phenomenon. On the other hand in the deconfinement phase at large distances the free energy also flattens off, the reason being a screening of the interactions between the quark and antiquark due to liberated colored medium degrees of freedom. The question of whether or how the screening properties of QCD may be captured by an analogous and equally simple Debye screening formula in analogy with the Abelian theory is a ongoing field of research. The properties of the correlation function of Polyakov loops at finite temperature in QCD have been thoroughly studied in lattice simulations [8][9][10][11][12][13]. More recently the Polyakov loop correlator on the lattice has been compared to effective field theory predictions, both in a perturbative setting in pNRQCD, as well as perturbatively matched EQCD [14]. For analytic studies of the Polyakov loop see e.g. [15,16].
While it is an interesting proposition to carry out similar studies of the Polyakov loop at finite Baryon density in QCD, the usual methods of lattice QCD unfortunately break down because of the so-called sign problem. Instead approaches, such as the Taylor expansion or analytical continuation (both in quark chemical potential) allow one to obtain useful results at small values of the chemical potential (see e.g. [17][18][19] ). There are a lot of analytical attempts to study properties of dense matter (see e.g. [20][21][22]. However, most of them are not based on first principles and it is not clear how to reliably estimate the systematic uncertainty of different models.
Instead of pursuing the question of finite density physics in QCD directly, we here turn to the study of theories, which are similar to QCD but are not plagued by the sign problem. We believe that in particular the study of dense two-color QCD [23,24] allows us to learn about the properties of regular QCD at nonzero chemical potential. Other candidate theories not further pursued here are e.g. QCD at nonzero isospin chemical potential [25][26][27]. Of course we cannot expect to obtain quantitative predictions from such a strategy, while vital qualitative insight may be gained.
Two-color QCD at finite chemical potential has been studied with lattice simulations quite intensively before, see, e.g. [28][29][30][31][32][33] and references therein. Mostly these papers are aiming at the study of the phase diagram of two-color QCD in the region of small and moderate baryon densities.
The phase structure of two-color QCD at large baryon densities was studied in our previous paper [34], where lattice simulations were carried out at a relatively small lattice spacing a = 0.044 fm. Compared to previous works, it allows us to extent the range of accessible values of the baryon density, up to quark chemical potential µ > 2000 MeV, avoiding strong lattice artifacts. The main result of the paper [34] is the observation of the confinement/deconfinement transition at finite density and low temperature. In view JHEP05(2019)171 of this finding we are interested in studying the properties of the Polyakov loop correlation function in cold dense quark matter and to shed light on how they are affected by the confinement/deconfinement transition, which takes place at finite density. This is the central question addressed in this paper.
The manuscript is organized as follows. In the next section 2 we describe our lattice set-up. In section 3 we give an overview of the present status of the cold dense two-color QCD phase diagram. Section 4 is devoted to the calculation of the renormalized Polyakov loop correlation functions, as well as of the color averaged, color singlet and color triplet grand potentials. In addition, in the same section we determine the renormalized Polyakov loop and the grand potential of a single quark/antiquark. Using the color singlet grand potential, the quark number and internal energy induced by a static quark-antiquark pair are obtained in section 5. We consider the string breaking phenomenon in dense quark matter in section 6 and Debye screening in section 7. In the last section 8 we discuss our results and conclude.

Simulation details
In our lattice study we used the tree level improved Symanzik gauge action [35,36]. For the fermionic degrees of freedom we used staggered fermions with an action of the form whereψ, ψ are staggered fermion fields, a is the lattice spacing, m is the bare quark mass, and η ν (x) are the standard staggered phase factors: η 1 (x) = 1, η ν (x) = (−1) x 1 +...+x ν−1 , ν = 2, 3, 4. The chemical potential µ is introduced into the Dirac operator (2.2) through the multiplication of the links along and opposite to the temporal direction by factors e ±µa respectively. This way of introducing the chemical potential makes it possible to avoid additional divergences and to reproduce well known continuum results [37]. In addition to the standard staggered fermion action we add a diquark source term [28] to equation (2.1). The diquark source term explicitly violates U V (1) and allows to observe diquark condensation even on finite lattices, because this term effectively chooses one vacuum from the family of U V (1)-symmetric vacua. Typically one carries out numerical simulations at a few nonzero values of the parameter λ and then extrapolates to zero λ. Notice, however, that this paper is aimed at studying the region of large baryon density where lattice simulations are numerically very expensive. For this reason, in this paper we have chosen a different strategy. Most of our lattice simulations are conducted at a single fixed value λ = 0.00075. In order to check the λ-dependence of our results for chemical JHEP05(2019)171 potentials aµ = 0.0, 0.1, 0.2, 0.3, 0.4 we carry out additional lattice simulations for values of λ = 0.0005, 0.001.
Integrating out the fermion fields, the partition function for the theory with the action S = S G + S F can be written in the form which corresponds to N f = 4 dynamical fermions in the continuum limit. Note that the pfaffian P f is strictly positive, such that one can use Hybrid Monte-Carlo methods to study this system. First lattice studies of the theory with partition function (2.3) have been carried out in the following papers [30,38,39]. In the present study we are going to investigate instead a theory with the partition function which corresponds to N f = 2 dynamical fermions in the continuum limit. It is known that the symmetries of the staggered fermion action are different from those of two-color QCD with fundamental quarks [28]. In particular, the symmetry breaking pattern of QC 2 D with fundamental quarks is SU Notice, however, that in the naive continuum limit for the staggered action with the diquark source term, the action factorizes into two copies of N f = 2 fundamental fermions [32]. In addition, for sufficiently small lattice spacing a the β-function of the theory (2.4) measured in [32] corresponds to the β-function of QC 2 D with two fundamental flavors. For these reasons we believe, that the partition function (2.4) after the rooting procedure corresponds to QC 2 D with N f = 2 fundamental fermions with a correct continuum chiral symmetry breaking pattern.
The results presented in this paper have been obtained in lattice simulations performed on a 32 4 lattice for a set of the chemical potentials in the region aµ ∈ (0, 0.5). At zero density we performed scale setting using the QCD Sommer scale r 0 = 0.468(4) fm [40]. In this case the string tension associated to µ q = 0 amounts to √ σ 0 = 476(5) MeV at a = 0.044 fm. Numerical simulations in the region of large baryon density require considerable computer resources. For this reason, for the present paper we performed our study at a pion mass of m π = 740(40) MeV, where the cost is manageable. We will preferentially choose a smaller pion mass in future simulations.
To calculate Wilson loops we have employed the following smearing scheme: one step of HYP smearing [41] for temporal links with the smearing parameters according to the HYP2 parameter set [42] followed by 100 steps of APE smearing [43] (for spatial links only) with a smearing parameter α AP E = 2/3. We found that the HYP2 parameter set provides a better signal-to-noise ratio for Wilson loops and correlators of the Polyakov loops than the HYP1 set. As for the calculations of Polyakov loop, color-averaged (2.5), colorsinglet (2.6) and color-triplet (2.7) correlators one step of HYP2 smearing for temporal JHEP05(2019)171 μ m π/2 1 GeV hadronic phase deconfinement dense quark matter BEC BCS Figure 1. Schematic phase diagram of dense two-color QCD at low temperatures.
links was performed, but in the case of singlet and triplet correlators the Coulomb gauge without residual gauge fixing was fixed at first. The formulas used in the calculation are The reason why we performed smearing is that, due to the large time extension (L t = 32) correlators of the Polyakov loops by default are very noisy. Thus one has to introduce some smoothing technique to extract the signal. An analogous smearing scheme was applied in the papers [44,45].
3 The phase diagram of dense two-color QCD at low temperatures Let us explore the tentative phase structure of two-color QCD as basis for our study of interquark interactions. Based on symmetry arguments it is possible to build a chiral perturbation theory (ChPT) for sufficiently small chemical potential [23][24][25]. This ChPT can be used to predict the phase transitions at sufficiently small values of the chemical potential. In particular, it was predicted that for small values of chemical potential (µ < m π /2) the system is in the hadronic phase. In this phase the system exhibits confinement and chiral symmetry is broken. At µ = m π /2 = 370(20) MeV (aµ 0.08 in lattice units) there is a second order phase transition to a phase where scalar diquarks form a Bose-Einstein condensate (BEC phase). The order parameter of this transition is the diquark condensate q T (Cγ 5 ) × τ 2 × σ 2 q , where Cγ 5 is the matrix which acts on Dirac indices and τ 2 , σ 2 are Pauli matrices which act on flavor and color indices of the quark field q. In the massless limit there is no chiral symmetry breaking, if the diquarks are condensed. However, for massive quarks the chiral condensate is not zero. Instead it is proportional to the quark mass and decreases with increasing chemical potential. Let us note that dense QC 2 D in the hadronic phase and the BEC phase was intensively studied within lattice simulations in a series of papers [28,32,33,38,39,46] where reasonable agreement with ChPT was observed.

JHEP05(2019)171
In the ChPT the interactions between different degrees of freedom are accounted for by perturbation theory, so they are assumed to be weak and in addition the baryon density in the region of application is small. Together with the fact that in two-color QCD the diquarks are baryons this implies that the system is similar to a dilute baryon gas at µ above µ ≥ m π /2 but below the values corresponding to large density.
Enhancing the baryon density further, we proceed to dense matter, where the interactions between baryons cannot be accounted for within perturbation theory. This transition can be seen through the deviation of different observables from the predictions of ChPT. In our paper [32] this deviation was observed in the diquark condensate, the chiral condensate and the baryon density.
At sufficiently large baryon density (µ ∼ 1000 MeV, aµ ∼ 0.22) some observables of the system under study can be described using Bardeen-Cooper-Schrieffer theory (BCS phase). 1 In particular, the baryon density is well described by the density of non-interacting fermions which occupy a Fermi sphere of radius r F = µ. In addition, the diquark condensate, which plays the role of a condensate of Cooper pairs, is proportional to the Fermi surface. In lattice simulation the BCS phase was observed in the following papers [31,32,47,48] and we found that the transition from the BEC to the BCS phase is smooth [32].
It is worth to note that the Cooper pairs condensate in the BCS phase leads to the mass gap in the fermion spectrum. The mass gap is important parameter which determines a lot of properties of dense quark-gluon matter. Due to the nonperturbative nature of QCD it is not clear how the mass gap depends on the chemical potential. One can expect that in the region of ultrahigh density the strong coupling constant becomes sufficiently small and the gap can be estimated as ∆(µ) ∼ µg −5 exp (−3π 2 / √ 2g) [49]. An important property of this approximation is that at sufficiently large density the ∆(µ) is rising function of the chemical potential. As the result starting from some value of the chemical potential ∆(µ) Λ QCD , i.e. dynamical quarks become heavy and the system in some properties resembles quenched QCD. Most probably in our simulations we have not reached the region ∆(µ) Λ QCD , but we might reach the region where ∆(µ) ∼ Λ QCD . We are going to consider this issue below.
In addition to the transition to the BCS phase at µ ∼ 1000 MeV (aµ ∼ 0.22) there is the confinement/deconfinement transition in dense two-color QCD [34]. This transition manifests itself in a rise of the Polyakov loop and vanishing of the string tension. It was also found that after deconfinement is achieved, we observe a monotonous decrease of the spatial string tension σ s which ends up vanishing at µ q ≥ 2000 MeV (aµ ≥ 0.45). It should be noted that the results of this study suggest that the confinement/deconfinement transition is rather smooth. The Polyakov loop results do not allow us to locate the transition region from the confinement to the deconfinement phase. For this reason we consider here the transition region to be around µ = 1000 MeV. This value was found in our previous study [34], where it was determined by the condition that the string tension, extracted from the Wilson loops, becomes zero within the uncertainty of the calculation. Thus throughout the paper we use the term "the confinement/deconfinement transition region" in the sense of vanishing of the string tension extracted from the Wilson loops.

JHEP05(2019)171
4 The grand potential of a static quark-antiquark pair in dense quark matter In this section we are going to study the grand potential Ωq q (r, µ) of a static quark-antiquark pair placed within a distance of r into the dense medium. It can be represented in terms of the correlator of Polyakov loops whereTr = 1 2 Tr and the Polyakov loop is given as the trace of a product of gauge links in The quantity c(µ) denotes a divergent renormalization constant, which is related to the self-energy of a quark or antiquark source. In the limit r → ∞, the correlation between the Polyakov lines becomes negligible and the grand potential Ω ∞ (µ) is given by the squared expectation value of the volume-averaged To find the grand potentials from formulae (4.1) and (4.2) one has to determine the renormalization constant c(µ).
In pure gauge theory the expectation value of the Polyakov line which is defined as is the order parameter of the confinement/deconfinement transition. In particular, L ren (µ) vanishes in the confined phase, whereas it is non-zero in the deconfined phase. After inclusion of dynamical quarks in the simulations, the expectation value of the Polyakov line is no longer an order parameter. However, one can interpret the Ω ∞ (µ)/2 as the grand potential of one quark or one antiquark in dense quark matter. Thus one may expect that in the confined phase Ω ∞ (µ) is much larger than that in the deconfined phase. Below we will also need the color-singlet grand potential Ω 1 (r, µ), which is defined as Notice that contrary to the color averaged grand potential Ωq q (r, µ), the singlet one Ω 1 (r, µ) is not gauge invariant. So, in order to calculate Ω 1 (r, µ) we have to fix the gauge and we choose here conventionally the Coulomb gauge. Both the color averaged and the color singlet grand potentials are calculated up to renormalization constants. Now let us define the relative normalization of these observables. It is clear that at sufficiently large spatial separation between quarks the relative orientation of charges in color space is not important due to screening. For this reason the authors of [8][9][10] chose a relative normalization of the color averaged and color singlet free energies, such that they are identical at large distances. In our paper we are going to use the same relative normalization between the Ωq q (r, µ) and Ω 1 (r, µ).

JHEP05(2019)171
For two-colors, the color averaged grand potential Ωq q (r, µ) can be represented [50] through the color singlet Ω 1 (r, µ) and the color triplet grand potential Ω 3 (r, µ) as Let us consider short distances (rµ 1), i.e. distances where the Debye screening can be neglected. In this limit the running of the coupling constant is determined by the scale ∼ 1/r, and the influence of the chemical potential on the running coupling can be neglected. The perturbative one-gluon exchange expression for the color singlet and the triplet grand potentials at short distances (rµ 1) can be written as We have already discussed relative normalization between the grand potentials Ωq q (r, µ) and Ω 1 (r, µ). Thus to renormalize the grand potentials it is sufficient to renormalize one of them. Let us consider Ω 1 (r, µ). To do this we use the procedure proposed in [8][9][10], adopted here for the calculation at finite density. The grand potential at finite temperature and chemical potential is defined as where function U 1 (r, T, µ) is the internal energy, the function S 1 (r, T, µ) is the entropy and the function N 1 (r, T, µ) is the quark number density of a static color singlet quarkantiquark pair. From the above discussion it is clear that at short distances r the grand potential does not depend neither on the chemical potential µ nor on the temperature T . This implies that at short distances the entropy S 1 = −∂Ω 1 /∂T and the quark number density N 1 = −∂Ω 1 /∂µ are zero. It is also clear that at short distances the internal energy equals to the interaction potential in a quark-antiquark pair at zero temperature and density. So, at short distances the grand potential Ω 1 (r, T, µ) coincides with the zero temperature and density potential V (r) which is calculated in appendix A. Similarly to papers [8][9][10] we fix the renormalization constant c (µ) through the matching condition for Ω 1 (r, µ) at short distances to the short distance behavior of the interaction potential V (r). The renormalization for the grand potential Ωq q (r, µ) can be fixed using matching at large distances, r where the color averaged and the color singlet grand potentials are expected to be identical. Evidently, this procedure allows us to get rid of the divergent self-energy contributions and uniquely fixes the renormalization constants c(µ) and c (µ).
Having gone through the renormalization procedure we are ready to present the results of the calculation of the renormalized grand potentials Ωq q (r, µ), Ω 1 (r, µ), Ω 3 (r, µ). In figure 2 we plot the renormalized Ω 1 (r, µ) as a function of distance for different values of the chemical potential. In figure 3 we plot the grand potential Ωq q (r, µ). To get an idea how the Ωq q (r, µ), Ω 1 (r, µ), The color singlet grand potential as a function of distance for few values of the chemical potential under study. The black curve is the potential of a static quark-antiquark pair at zero density and temperature. Note the absence of a Coulombic small distance regime, due to smearing.  The black curve is the potential of the static quark-antiquark pair at zero density and temperature.Note the absence of a Coulombic small distance regime, due to smearing.
can be extracted from the Polyakov loop correlator at large distances. In the calculation we take Ω 1 (∞, µ) = Ω 1 (L s /2, µ) and calculate the renormalized Polyakov loop applying formula (4.3). Notice that for the calculation it is important that the grand potential extracted from the Polyakov loop correlator goes to a plateau value. In the confined, phase the plateau in the grand potential is due to string breaking, which takes place at large distance. Due to the relatively small spatial lattice size we can observe the string breaking only for sufficiently large chemical potential (µ > 440 MeV). For this reason the calculation of the Ω 1 (∞, µ) and L ren (µ) based on the renormalized correlator of the Polyakov loops can be carried out for µ > 440 MeV. The results for Ω 1 (∞, µ) and L ren (µ) are shown in  we conduct one step of the HYP smearing. The Polyakov loop is renormalized according to where L ren (µ = 1030 MeV) is the Polyakov loop measured in the previuos approach, 2 and L bare (µ) is the bare Polyakov loop measured on the lattice. Similarly to the renormalization of the correlators of the Polyakov loops, the approach based on (4.8) gets rid of infinite ultraviolet divergence and uniquely fixes the renormalization.
Having calculated the renormalized Polyakov line, we can find the Ω 1 (∞, µ) using formula (4.3). The results for Ω 1 (∞, µ) and L ren (µ) are shown in figure 5 and in figure 6 by blue circles. From these figures one sees that both approaches to the calculation of the Ω 1 (∞, µ) and L ren (µ) are in agreement with each other.
Here a few comments are in order: the measurement of the Polyakov loop correlation functions in this section was carried out at λ = 0.00075. As was discussed above, in order to check the λ-dependence of our results we carried out a similar study at aµ = 0.0, 0.1, 0. The confinement/deconfinement transition at finite temperature manifests itself in an increasing value of the Polyakov loop and its rise may become quite rapid in the transition region [10]. A similar behaviour can be seen from figure 6, where the confinement/deconfinement transition is observed through the rise of the Polyakov line. At the same time from this figure we don't see any specific region in the chemical potential where the rise of the Polyakov line is dramatically different as compared to other regions. This observation corroborates our previous finding that the finite density confinement/deconfinement transition in two-color QCD is rather smooth. 3

JHEP05(2019)171
Let us put these finding into further context by considering the ideas proposed in [25], where the authors proposed an interesting scenario for the confinement/deconfinement transition in dense matter. Although the paper is devoted to QCD at finite isospin density, this theory is very similar to two-color QCD at finite baryon density. So the results of [25] are very likely to be applicable to two-color QCD as well. The idea can be summarized as follows. At sufficiently large density the mass gap ∆(µ) becomes much larger than the confinement scale Λ conf , which differs from Λ QCD at zero density. The critical temperature at which the diquark condensate melts is of order of ∆(µ). The critical temperature at which confinement/deconfinement transition takes place is of order of Λ conf . So, at high density the confinement/deconfinement transition takes place at much smaller temperature than melting of the diquark condensate. The confinement/deconfinement transition on the (T, µ) phase diagram is almost a horizontal line which might have small slope (see the tentative phase diagram in figure 1 of paper [25]). The lattice study presented in this paper has been carried out at small temperature through the variation of the chemical potential, i.e. on the horizontal line on the (T, µ) phase diagram. The confinement/deconfinement transition observed in our study might result from the cross of these two lines. This scenario might also explain large width of the transition region and its smoothness. While the currently available data is consistent with parts of the above argument we are yet unable to unambigously confirm or reject the underlying hypotheses.
Let us now pay attention to the region µ > 2000 MeV. In this region the Polyakov loop/grand potential reaches a maximum/minimum and then drops/rises. Below it will be shown that the region µ > 2000 MeV differs from the region µ < 2000 MeV not only for the Polyakov loop but also the grand potential. In turn also derived observables, such as the screening length R sc , the Debye mass and effective coupling constant show a distinctive behavior.
At the moment we do not fully understand the physics, which is responsible for this behavior. One possibility is that the value of the chemical potential µ ∼ 2000 MeV is exceptional since there is nonzero spatial string tension in the region µ < 2000 MeV whereas the spatial string tension is zero for µ > 2000 MeV. This might imply that the point µ ∼ 2000 MeV separates systems with and without magnetic screening.
Another hypothesis which describes the region µ > 2000 MeV may be the following. We have already discussed that at the ultrahigh density the ∆(µ) Λ QCD , dynamical quarks become heavy and the system is similar to quenched QCD. The last property implies that at the ultrahigh density two-color QCD returns to confinement phase, where the Polyakov line is small. At zero density the Polyakov line is also small, but it starts to rise at nonzero density. These facts allow us to state that at some value of the chemical potential the Polyakov line reaches a maximum. One can expect that to the left of this maximum ∆(µ) < Λ QCD and the quarks can be considered as light degrees of freedom. To the right of this maximum ∆(µ) > Λ QCD and the quarks are heavy. So, it is reasonable to assume that the maximum at µ ∼ 2000 MeV takes place in the region where ∆(µ) ∼ Λ QCD . This hypothesis is supported by our data for the Debye mass (see below). However our data do not allow us to choose if one of these hypotheses is correct and further study is required to settle this question.

JHEP05(2019)171
From figure 5 one sees that the potential Ωq q (r, µ) changes its sign at µ ∼ 1300 MeV, which at first sight may be unexpected. However, let us recall that Ωq q (r, µ) is not a grand potential of the whole system. On the contrary, it is the difference of the grand potential of dense quark matter with a static quark-antiquark pair and dense quark matter without a static quark-antiquark pair. So, Ωq q (r, µ) in our context is the additional grand potential due to the introduction of the quark-antiquark pair to quark matter. Now figure 5 can be interpreted as follows: introducing a static quark-antiquark pair increases the grand potential of the system for µ < 1300 MeV and decreases the grand potential of the system for µ > 1300 MeV. An explanation of this fact will be presented in the next section.
The authors of [31] studied QC 2 D with N f = 2 quarks within lattice simulation with dynamical Wilson fermions. In particular, they measured the Polyakov loop as a function of the chemical potential and observed the following behavior: it remains zero up to aµ ∼ 0.75 and then quickly rises. In the region aµ > 1, due to the saturation the quark degrees of freedom, quarks are no longer dynamical and the theory becomes quenched QCD and exhibits confinement, i.e. the Polyakov loop goes to zero for aµ > 1.
Further measurement of the string tension carried out in [51] did not confirm the presence of a confinement/deconfinement transition and the decrease of the string tension with chemical potential. Although the behavior of the Polyakov line in figure 6 seems similar to that obtained in [31], we believe that this behavior can be explained as a lattice artifact of Wilson fermions. In order to explain the results of [31], let us recall that in Wilson fermions one has one light quark and 15 heavy quark species with masses ∼ 1/a. If the chemical potential is aµ ∼ 1 the heavy quarks are not suppressed any longer and additional color degrees of freedom are released to the system under study. We believe that this mechanism is responsible for the rise of the Polyakov line observed in [31]. Notice that this mechanism does not work for staggered quarks as used in our paper, since there are no heavy species. In addition we observe a considerable rise of the Polyakov loop already at aµ ∼ 0.2. The decrease of the Polyakov line in figure 6 hence cannot be attributed to the saturation since it starts at aµ ∼ 0.4, which is rather far away from the saturation in staggered fermions [32].

The quark number and internal energy induced by a static quark-antiquark pair
The authors of [11] calculated the free energy of a static quark-antiquark pair in QCD at finite temperature. In addition they calculated the entropy of the QCD medium in the presence of a static-quark antiquark pair. In dense quark matter there is in addition a contribution of the entropy to the grand potential (4.7). However, this contribution is not important at low temperature. What becomes important in dense quark matter is the quark number induced by the static quark-antiquark pair N (r). This quantity can be calculated as follows N (r, µ) = − ∂Ω(r, µ) ∂µ .  Notice that the N (r, µ) is the quark number which arises in the system due to the introduction of the quark-antiquark pair to the dense matter. So, it is the difference of the quark number with and without static quark-antiquark pair.
In this paper we calculate the N 1 (r, µ) for the color singlet grand potential. Similar observables can be calculated for the color averaged and color triplet grand potentials. To find the N 1 (r, µ) we determine the derivative of the grand potential over chemical potential through the finite difference approximation. The results of the calculation of the N 1 (r, µ) for a few values of the chemical potential are shown in figure 7. We found that the smaller the chemical potential, the larger uncertainty of the calculation in the N 1 (r, µ). For this reason we show only for N 1 (r, µ) for sufficiently large chemical potentials. From figure 7 one sees that N 1 (r, µ) is rising from zero at short distances to some plateau value N 1 (∞, µ), which is an important observable, since it is proportional to the derivative of the Polyakov loop of a single quark/antiquark in a dense medium over the chemical potential. One might expect that at the critical chemical potential where the confinement/deconfinement phase transition takes place there is an inflection point of the Polyakov loop. For this reason the confinement/deconfinement phase transition might manifest itself in the maximum of N 1 (∞, µ). For the same reason the authors of [11] observed a maximum of the entropy at the critical temperature. Our results for N 1 (∞, µ) are shown in figure 8. Unfortunately due to large uncertainties of the calculation we are unable to locate this maximum. If we ignore the entropy contribution to the grand potential, which is small in cold dense matter, one can calculate also the internal U (r, µ) energy using the formula U (r, µ) = Ω(r, µ) + µN (r, µ). (5.2) In this paper we calculate the U 1 (r, µ) for the color singlet grand potential, the result of which is shown in figure 9.

String breaking in dense quark matter
Let us consider again figure 2 and figure 3. We have already mentioned that in dense QC 2 D the confinement/deconfinement transition takes place at µ ∼ 1000 MeV. Despite this fact from figure 3 we see that Ωq q (r, µ) reaches the plateau already at µ = 447 MeV. This happens because of the string breaking phenomenon, which for µ = 447 MeV takes place at r ∼ 0.5 fm. Of course string breaking occurs also for smaller chemical potentials, but we do not observe it, as it takes place beyond our spatial lattice size. one also sees that the larger the chemical potential, the smaller the distance at which the string breaking takes place. Let us consider this phenomenon more quantitatively. 4 To study the string breaking phenomenon we introduce the screening length R sc which can be calculated from the solution of the equation [8] V µ=0 (R sc ) = Ωq q (∞, µ), (6.1) where V µ=0 (r) is the static potential at zero density. For Ωq q (∞, µ) we take the grand potential calculated from the renormalized Polyakov loop measured on the lattice(see figure 5). The results of this calculation are shown in figure 10. This plot tells us that the larger the chemical potential the smaller the string breaking distance. In order to understand this behaviour, let us recall that in three-color QCD the string breaking phenomenon can be explained by the possibility to break the string between static quarks by a quark-antiquark pair created from vacuum. If the length of the string is larger than the critical one it becomes energetically favorable to break the string and form two heavy-light meson instead of increasing the length of the string.
In dense two-color QCD in addition to the possibility to break the string by quarkantiquark pairs it becomes possible to break the string by two quarks. As the result of 4 At sufficiently large values of the chemical potential the screening properties of the dense medium are described by the Debye screening phenomenon rather than by the string breaking. In this section we consider the chemical potential values corresponding to nonzero string tension extracted from the Wilson loops [34].

JHEP05(2019)171
this phenomenon, after the string has been broken, one ends up with a heavy-light meson and one heavy-light diquark. Due to confinement, the two quarks have to be extracted from some hadron. The two-color baryon -the scalar diquark is a good candidate for such a hadron. Indeed at nonzero µ the scalar diquark is a lightest state in the system. At µ > m π /2 there is condesation of the scalar diquarks, so the two quarks can be extracted from the diquark condensate. This picture is supported at large µ in the BCS phase, where one has a Fermi sphere with radius µ. Evidently one cannot break the string by taking two quarks deep inside the Fermi sphere, since in that case, the quarks which break the string due to the interactions have to move from one point of the Fermi sphere to some other point inside the Fermi sphere. However, all points inside the Fermi sphere are occupied. So, the only possibility to break the string is to take two quarks close to the Fermi surface. In the confined phase, quarks on the Fermi surface are condensed as diquarks. Thus we again confirm the picture that two quarks, which break the string between a quark-antiquark pair, can be taken from the available diquarks.
It is clear that the QCD string breaks when the energy accumulated in the string σr becomes larger than double binding energy E b of light quark in the field of static color source. This binding energy non-trivially depends on the mass gap in the fermion spectrum ∆(µ) and on dynamics of QCD. One can expect that for the ∆(µ) . So, at sufficiently large density the string breaking phenomenon is determined by the mass gap in the fermion spectrum.
Further, let us consider the following model of string breaking: if one diquark penetrates inside the string, it breaks the string with some probability ω. It is clear that the ω(∆) is a decreasing function of the mass gap ∆(µ). Notice that if the ∆(µ) Λ QCD the ω(∆) is determined by QCD dynamics and weakly depends on ∆(µ) and µ. For the density of diquarks n, the string length R and the transverse area S the number of diquarks inside the string is n × R × S. If the string breaking events are independent, the total probability to break the string P ω(∆) × n × R × S. The condition for the string breaking is P ∼ 1. From last statement we conclude that R sc ∼ 1/(nω(∆)). From the derived formula one sees that there is a competition between the baryon density which rises with the chemical potential and ω(∆) which decreases with µ. From figure 10 it is seen that up to the µ ∼ 2000 MeV the R sc is decreasing function of the chemical potential and up to the deconfinement phase, where there is string breaking, the density wins. So, one can expect that up to µ ∼ 1000 MeV the quarks are light and ∆(µ) < Λ QCD .
If one increases the chemical potential then at some density R sc becomes so small that the string cannot be created, i.e. at the instant of creation it will be immediately broken by the two-color baryons -diquarks. This is our hypothesis of the deconfinement mechanism in two-color dense quark matter. It is not clear how to find unambiguously the distance at which the string ceases to be stable. From the interaction potential at zero density (see figure 2) it is found that this happens in the region r ∈ (0.2, 0.3) fm. Using figure 10 one can infer that the interactions in this interval are screened, as the chemical potential lie within µ ∈ (900, 1300) and which agrees with the position of the confinement/deconfinement transition. We believe that this fact confirms our hypothesis.

JHEP05(2019)171
For a chemical potential larger than the µ ∼ 1300 MeV(aµ ∼ 0.3) the R sc is smaller than 0.2 fm. At such small distances the entropy S = −∂Ω/∂T , the quark number density N = −∂Ω 1 /∂µ are small and the grand potential is mainly determined by the interaction potential at zero temperature and density. At the same time the renormalized interaction potential (see figure 2) is negative at distances r < 0.2 fm. For this reason the grand potential of one quark in dense quark matter becomes negative, which was observed in the previous section.
In addition to the R sc in figure 10 we plot the average heavy quarkonia J/Ψ, χ c , ψ radii which where estimated in appendix B within a simple potential model. It is clear that if the screening length is close to the heavy quarkonium radius this state is considerably modified by dense quark matter. From figure 10 one sees that the heaviest state the ψ due to its rather large radius should be considerably modified at nonzero density before the transition to BEC phase. The χ c meson will instead be modified in the BEC phase. Finally we predict that the J/Ψ meson will be modified in dense quark matter but below the deconfinement transition. Notice, however, that if the radius of a charmonium equals to the R sc at some density n 0 , the dissociation of this charmonium takes place at densities larger than n 0 .
A more detailed study of quarkonium dissociation in two-color dense quark matter will be presented in a future study. In particular the question of the presence of an imaginary part in the interquark potential at finite density, which may further destabilize the bound states will be carefully investigated.

Debye screening in dense quark matter
In the region µ > 900 MeV the system under study transitions from the confined to the deconfined phase. In the deconfined phase the contribution of the string is markedly reduced and one may attempt to describe R sc in a dense quark-gluon plasma via an analogy with the Abelian theory, i.e. purely Coulombic Debye screening.
The scale of the Debye screening in perturbation theory is denoted by the Debye mass, which to one-loop order (for the N c = 2) reads To describe or results for the R sc it is reasonable to assume that the screening length is inversely proportional to m D (µ). For this reason we fit our data by the formula where the A is some factor. We fit our data in the region µ ∈ (900, 1800) MeV and use a two-loop approximation for the running of the coupling constant α s (µ) (see formula (B.2) with N f = N c = 2). The fit describes our data well (χ 2 /dof 0.8) and the best fit parameters are A = 1.4 ± 0.4, Λ = 140 ± 80 MeV. In the region µ > 2000 MeV the R sc goes to plateau and the data cannot be described by the formula (7.2). Now let us study how the Debye screening phenomenon manifests itself in the large distance behavior (rµ 1) of the grand potential. In this case the dominant scale is the chemical potential, i.e. the running coupling constant depends only on µ: g(r, µ) = g(µ). For sufficiently large density one can apply perturbation theory to calculate grand potentials. Perturbatively the grand potential Ωq q (r, µ) is determined by two-gluon exchange and it is rapidly decreasing with distance function. Contrary to Ωq q (r, µ) the color singlet grand potential Ω 1 (r, µ) is determined by one-gluon exchange. In this paper we consider only Ω 1 (r, µ), whose leading order contribution has the form where m D is the Debye mass given by the expression (7.1). It tells us that due to Debye screening at sufficiently large distance the expression (Ω 1 (∞, µ) − Ω 1 (r, µ))r is an exponentially decreasing function of the distance. We plot (Ω 1 (∞, µ) − Ω 1 (r, µ))r in logarithmic scale in figure 11. From this figure the exponential decrease at large distance is seen starting from the µ = 850 MeV what confirms Debye screening phenomenon in deconfined dense quark matter. The deviation from a purely Coulombic Debye-like behavior at intermediate distances may be related to the remnants of the string, which is not perfectly screened.

JHEP05(2019)171
Further we fit our data in the deconfinement phase for Ω 1 (r, µ) at sufficiently large r by the formula (7.3). The results for m D /µ and α s (µ) as a function of the chemical potential are shown in figure 12 and figure 13. From figure 12 it is seen that the dependence of the Debye mass on the chemical potential is m D ∼ µ. Due to large uncertainties of the calculation, we are not able to resolve the running of the coupling constant with µ. For the same reason we are not able to observe the running of the α s , as shown in figure 13. The running coupling is constant within the uncertainty of the calculation up to the µ < 2000 MeV and it starts to drop in the region µ > 2000 MeV. The ratio m D /µ as well as the ratio 1/(R sc µ) ∼ m D /µ starts to drop in the region µ > 2000 MeV.
We have already mentioned that the maximum of the Polyakov line might be explained by the hypothesis that in the region µ < 2000 MeV the quarks are light and in the region µ > 2000 MeV they become heavy. The drop of the ratio m D /µ might be explained by the same hypothesis as follows. The one-loop formula (7.1) was derived assuming that the quarks are light. However, if in the region µ > 2000 MeV they become heavy due to the mass gap, the Debye mass is smaller than what is given by (7.1). Moreover, one can expect that the larger the chemical potential the larger the mass gap and the smaller the ratio m D /µ. higher order radiative corrections. In the deconfined phase at finite temperature and zero density one also obtains a large coupling constant (see e.g. [9,10]). Despite the large coupling constant, it turns out that the one-loop formula (7.1) works quite well. To show this, we compute the ratio m D /(µ √ α s ). At the leading order approximation this ratio is 4/π. In figure 14 we plot this ratio and find that within the uncertainties of the calculation formula (7.1) works quite well for m D and α s extracted from the color singlet grand potential.

Conclusion and discussion
In this paper we continued our study of two-color QCD at finite density and low temperature based on lattice simulations. Our simulations were performed on 32 4 lattices with rooted staggered fermions at a relatively small lattice spacing a=0.044 fm, which allowed us to study two-color QCD at very large baryon densities (up to quark chemical potential µ > 2000 MeV) while avoiding strong lattice artifacts. The aim of the present paper was the study of the interaction between a static quarkantiquark pair in two-color dense quark matter. To this end we performed computations of the Polyakov loop correlation functions and calculated the color averaged, color singlet and color triplet grand potentials. To handle appropriately the divergent self-energy contribution to the Polyakov loop correlation functions, we conduct renormalization through a matching of the color singlet grand potential to the static interaction potential of quarkantiquark pair at short distances. Having determined the renormalized grand potentials, we calculated the renormalized grand potential of a single quark/antiquark and average

JHEP05(2019)171
Polyakov loop. In addition we calculated the quark number induced by a static quark antiquark pair and its internal energy.
The confinement/deconfinement transition at finite density manifests itself in an increasing value of the Polyakov loop. The finite density transition does not show a region of rapid rise of the Polyakov loop contrary to the finite temperature case. For this reason we conclude that the transition from confinement to deconfinement at finite density is smooth.
The numerical results obtained in this paper have further been considered in light of the proposed scenario for the confinement/deconfinement transition in QCD-like theories presented in [25]. At sufficiently large density the confinement/deconfinement transition is almost a horizontal line with small slope which is located far below the temperatures at which the diquark condensate melts (see the tentative phase diagram in the paper [25]). The lattice study presented in this paper has been carried out at small temperature through the variation of the chemical potential, i.e. on the horizontal line without slope. The confinement/deconfinement transition observed in our study might result from the cross of these two lines. This scenario might also explain large width of the transition region and its smoothness. While the currently available data is consistent with parts of the above argument we are yet unable to unambiguously confirm or reject the underlying hypotheses.
In addition we calculated the screening length R sc which is defined as where V µ=0 (r) is the static potential at zero density and the Ω(∞, µ) is the grand potential of a static quark-antiquark pair at infinite distance. In the confined phase, the screening length is determined by the string breaking length, whereas in the deconfined phase R sc is determined by the Debye screening phenomenon. The result of the calculation of the screening length shows that consistent with intuition, the larger the chemical potential the smaller R sc , the string breaking distance. We believe that the decrease of the string breaking distance with density can be attributed to a further string breaking mechanism in dense matter. In dense two-color QCD, in addition to the possibility to break the string by a quark-antiquark pair, it becomes possible to break the string by two quarks which can be extracted from a two-color baryon -the scalar diquark. As the result of this phenomenon, after the string breaking one end up with one heavy-light meson and one heavy-light diquark. Lattice studies show [32] that in the region µ > m π /2 the scalar diquark condensate increases with the chemical potential, i.e. it becomes easier to find two quarks and to break the string.
If one increases the chemical potential then at some density R sc becomes so small that the string cannot be created at all. Once created it will be immediately broken by the two-color baryons -the scalar diquarks. This is our hypothesis of the deconfinement mechanism in two-color dense quark matter.
The behavior of the string breaking distance in dense matter and the deconfinement mechanism are not specific only for two-color QCD. We believe that a similar process can be realized in SU(3) QCD with the difference that one has to replace two-quark baryon in SU(2) by three-quark baryon in SU(3). In particular, one can expect that the screening length which has the same definition as in two-color QCD is decreasing function of the JHEP05(2019)171 chemical potential. In turn, the larger the density the smaller the string breaking distance. For three colors this behavior can be explained as follows: at nonzero chemical potential one has a nonzero baryon density in the system. Baryons which form this density can break the string, splitting it into one quark and one diquark. After the string breaking one has one heavy-light meson and heavy-light baryon. Finally the larger the the chemical potential, the larger the number of baryons which can break the string, i.e. the string breaking distance is a decreasing function of the chemical potential.
Notice that in three-color QCD one might also have a similar mechanism of deconfinement at finite density, as we proposed above for two-color. In particular, deconfinement takes place at the density at which the R sc is so small that the string can not be created.
In the previous section we considered the large distance behavior of the color singlet grand potential in the deconfined phase. In analogy with Debye screening in the Abelian theory and using leading order perturbation theory, we attempt to quantitatively describe the observed behavior and find good agreement with the lattice data.
We calculated the Debye mass and the coupling constant for various chemical potentials. The coupling constant extracted in this way takes on values α s ∼ 1, which tells us that despite the large baryon density, the system remains strongly coupled. It was also found that despite the large coupling constant the one-loop formula for the Debye mass works well at large distances within the uncertainty of the calculation.
In this paper we found that the region µ < 2000 MeV physically differs from the region µ > 2000 MeV. This manifests itself in different behavior of the following observables: the Polyakov line, the grand potential, the screening length R sc , the Debye mass and effective coupling constant.
While we do not yet fully understand the physics, which is responsible for this behavior, one possibility is that the value of the chemical potential µ ∼ 2000 MeV is exceptional since it divides the region µ < 2000 MeV with a spatial string tension from that at µ > 2000 MeV where it vanishes. This may imply that the point µ ∼ 2000 MeV separates systems with and without magnetic screening.
The other hypothesis about the region µ > 2000 MeV is the following. At the chemical potential µ ∼ 2000 MeV the mass gap ∆(µ) ∼ Λ QCD . For this reason to the left of the µ ∼ 2000 MeV the quarks are light and to the right of this value the quarks are heavy. The theory with heavy dynamical quarks is similar to quenched QCD what explains the difference between regions µ < 2000 MeV and µ > 2000 MeV. Unfortunately our data do not allow us to choose if one of these hypotheses is correct and further study is required to settle this question.
Finally we are going to discuss lattice artifacts which result from the saturation effect. It is known that at large values of the chemical potential aµ ∼ 1 a saturation effect starts to be seen. The essence of this effect is that all lattice sites are filled with fermionic degrees of freedom, and it is not possible to put more fermions on the lattice ("Pauli blocking"). It is known that the saturation effect is accompanied by the decoupling of the gluons from fermions. Thus, effectively due to saturation, the system becomes simply gluodynamics, which is confined at low temperatures. From this consideration it is clear that in order to study the properties of quark matter at large baryon density one should have sufficiently small lattice spacing such that the properties are not spoiled by this kind of artificial JHEP05(2019)171 confinement at large values of the chemical potential. We believe that because of the saturation effect the deconfinement in dense SU(2) matter has not been seen before.
The results of the study presented in this paper are obtained for the chemical potentials µ < 2200 MeV (aµ ≤ 0.5). We believe that our results are not spoiled by saturation for the following reasons. First, for the µ > 2000 MeV (up to µ ∼ 2500 MeV [34]) the spatial string tension is vanishing. Second, we do not see a respective rise of the timelike string tension. Moreover, the static potential for µ > 2000 MeV (up to µ ∼ 2500 MeV [34]) is well described by Debye screening potential. So, the properties of the system in the range µ > 2000 MeV are very different from those of plain gluodynamics at small temperatures.
Notice also that in our previous study of dense two-color QCD [32] we found that the onset of the saturation effects are seen at aµ ∼ 0.7-0.8. This was deduced through the decrease of the diquark condensate for aµ > 0.7 (while it is rising with µ for the aµ < 0.7). The rise of the diquark condensate in the continuum is related to the rise of the Fermi surface. The decrease of the diquark condensate on the lattice is evidently related to the onset of the saturation effect what can be seen as follows. Due to finite number of the fermion states in the lattice Brillouin zone there is a chemical potential from which the rise of the chemical potential does not lead to the rise of the Fermi surface on the lattice. Notice that at this value of the chemical potential not all fermion states on the lattice are filled and the saturation takes place at larger values of the chemical potential.
Finally the deviation of the lattice measured baryon density from the baryon density calculated for free fermions is 10 % for aµ = 0.45( 2000 MeV) and 20% for aµ = 0.50( 2250 MeV). We argue that such a deviation even if it could be attributed to saturation cannot lead to considerable modification of physics. Notice also that such a deviation may also be explained by other mechanism than saturation, e.g. the finite lattice spacing which is present for any aµ.
Taking into account all what is written above we believe that in the region under consideration in this paper, (aµ < 0.5) our results are not spoiled by eventual saturation effects. Notice that the strict proof of last statement requires additional lattice simulations at smaller lattice spacing, which are planned in the future.
A The interaction potential of static quark-antiquark pair at zero density In appendix A we are going to calculate the interaction potential of static quark-antiquark pair V (r) at zero density through lattice measurement of Wilson loops. In this paper we use the V (r) in order to renormalize the correlation functions of the Polyakov lines.
For the calculation of Wilson loops we have employed one step of HYP smearing [41] for temporal links with the smearing parameters according to the HYP2 parameter set [42], followed by 100 steps of APE smearing [43] for spatial links only with the smearing parameter α AP E = 2/3. The similar smearing scheme was applied in the paper [45] for the extraction of V (r) from the Wilson loops. In the case of spatial Wilson loops (see below) the smearing scheme was adopted respectively to consider one of the spatial directions as a "temporal direction". 5

JHEP05(2019)171
Having measured Wilson loops for all distances between static charges one can calculate the interaction potential V (r). Notice, however, that the interaction potential is determined up to a renormalization constant. We find this renormalization constant for the interaction potential at µ = 0 using the following procedure. It is known that for the distances r > 0.5 fm [52] the interaction potential is well described by the linear confinement potential corrected by the so-called Lüscher term which describes string fluctuations V (r) = σr − π 12r . (A.1) In order to calculate the interaction potential unambiguously we fit our lattice data by the potential V (r) = σr − π/12r + C. The fit is good and the constant is equal to C = 726 ± 13 MeV. The renormalization constant C determines the energy shift of the lattice potential as compared to (A.1). The renormalized lattice potential can be obtained through the shifting of the lattice potential down by C and it reproduces the potential (A.1) for r > 0.5 fm.
At the end of this section we would like to mention that the potential V (r) calculated as it was described above contains lattice artifacts at small distances r/a ≤ 3. These artifacts result from the HYP smearing which modifies the potential at small distances and the violation of the rotation invariance by our lattice at short distances. However, the measurements of the correlators of the Polyakov lines and the Wilson loops are carried out with the same HYP smearing and on the same lattice. For this reason, the interaction potential determined in this section is appropriate for the renormalization of the Polyakov lines correlators.

B Quarkonia properties in two-color QCD
In appendix B we are going to estimate the two-color quarkonia masses and their sizes using potential model. The background of all potential models is the interaction potential. In this paper we are going to use the Cornell potential [53]  Notice that the coefficient 3/4 is front of the Coulomb term is due to N f = 2 colors in our system. In the calculation we use the string tension σ = (476 MeV) 2 which was calculated in this paper. The effective coupling constant α s is extracted from the Cornell potential fit of our lattice data in the region 3 r/a 12. Thus we obtained the following value α s = 0.31. Notice that this value a larger that that in SU(3) theory α s = 0.21 [10]. We believe that the difference between two-and three-color QCD can be attributed to different RG running of the coupling constants. In particular, the two loops running of the coupling constant is given by the formula We use the non-relativistic Schrödinger equation with the potential B.1 in order to estimate the quarkonia masses and sizes. In the calculation we take the charm quark mass m c = 1850 MeV from the paper [53]. Having solved the equation numerically, one obtains the following values of the states energies and mean squared distance between quarks r 2 :