Dense two-color QCD from Dyson-Schwinger equations
Romain Contant Markus Q. Huber






contant.romain.phy@gmail.com
Institute of Physics, University of Graz, NAWI Graz, Universitätsplatz 5, 8010 Graz, Austria


markus.huber@physik.jlug.de
Institut für Theoretische Physik, Justus-Liebig–Universität Giessen, Heinrich-Buff-Ring 16, 35392 Giessen, Germany




We investigate quantum chromodynamics with two colors at nonvanishing density using Dyson-Schwinger equations.
Lattice methods do not have a complex action problem in this theory.
Thus, we can benchmark our results and the effect of truncations directly by comparing with the corresponding lattice results.
We do so for the gluon propagator, the chiral condensate, and the quark number density and test variations of the employed truncation to improve the agreement.
Finally, we compare the effect of a truncation on the chiral and confinement/deconfinement transitions in the phase diagrams of QCD and QCD with the gauge groups SU(2) and G_2.


12.38.Aw, 14.70.Dj, 12.38.Lg



Dense two-color QCD from Dyson-Schwinger equations
    Markus Q. Huber
    December 30, 2023
==================================================








§ INTRODUCTION


Strongly interacting matter possesses a rich phase structure.
At vanishing chemical potential, there is a crossover from a hadron dominated phase at low temperatures to the so-called quark-gluon plasma phase at high temperatures.
This is well established, for example, by lattice simulations <cit.>.
Switching on a chemical potential, though, the picture becomes less clear.
The otherwise successful method of Monte Carlo lattice simulations faces a complex action problem that currently restricts calculations to a baryon chemical potential of μ_B≲ 2T <cit.>.

In that regime, no critical endpoint has been found by lattice simulations and it is up to other methods to explore the regime beyond that.
For cold and dense strongly interacting matter, a rich phase structure is expected, but the details are still not clear besides that at high enough densities a color superconducting phase will appear <cit.>.

Thus, to explore high densities, alternative approaches to lattice simulations are required.
One such method is functional equations; for general reviews, see, e.g., <cit.>, and for their application to the phase diagram of quantum chromodynamics (QCD), e.g., <cit.>.

They provide access to the full phase diagram but two things need to be dealt with. First, from a technical point of view, nonvanishing density and temperature complicate the equations due to the loss of manifest Lorentz invariance. Second, being exact equations, they can only be solved once approximations are made because they form an infinitely large system of equations.
In the vacuum, several studies indicate that the system has favorable convergence properties <cit.>.
However, in the medium, the corresponding level of truncation is not reached yet, and once it is, it is not guaranteed that the convergence properties remain the same.
Thus, for the time being, calculations in the medium operate at a more basic level involving phenomenological models for interactions.

At vanishing chemical potential, 'intermediate' (gauge dependent) results, like propagators and vertices, as well as 'final' (gauge invariant) results, like condensates, of functional equations can be compared against results from lattice calculations.
Thus, truncations can be assessed directly.
At nonvanishing density, the possibilities for such comparisons are very limited.
One way is to switch to a theory which does not suffer from a complex action problem.
This is the path we follow here.
We will solve Dyson-Schwinger equations (DSEs) of QC_2D <cit.>, which is QCD with the gauge group SU(2) instead of SU(3).
For an even number of quark families it has no complex action problem.
Lattice calculations were performed for this theory <cit.>, and we will compare results for the gluon propagator, the quark number density, and the chiral condensate.
QC_2D was also investigated with various continuum methods, e.g., <cit.>.
Furthermore, it has attracted attention <cit.> as a potential theory for a composite Higgs in the context of technicolor <cit.>.

A second alternative to QCD is QCD with the gauge group G_2 for which lattice simulations do not have a sign problem either <cit.>.
We will also consider this theory shortly in Sec. <ref> where we compare the phase diagrams of QCD and the two QCD-like theories with respect to chiral symmetry and confinement.

The spectrum of QC_2D differs from real QCD as it does not have fermionic baryons.
Instead, it contains color-neutral diquarks which can condense and lead to the existence of a transition from a Bose-Einstein condensate (BEC) to a Bardeen-Cooper-Schrieffer (BCS) phase.
Diquarks can be taken into account with the Nambu-Gor'kov formalism <cit.>.
We will not include diquarks at this point as this would increase the complexity of the calculations.
We will, however, come back to them in Sec. <ref> and shortly touch upon their impact on our calculations.
Calculations with diquarks and a similar setup can be found in Ref. <cit.>.

In the next section we detail our setup.
The results are presented in Sec. <ref>, and we summarize and provide an outlook in Sec. <ref>.






§ SETUP



We solve the gluon and quark propagator DSEs using a model for the quark-gluon vertex and lattice input for the quenched gluon propagators.
The equations are shown in Figs. <ref> and <ref>.
As is made explicit in fig:AADSEApprox, the quark loop is split off in the gluon propagator DSE.
We approximate the rest by quenched lattice results and calculate the quark loop explicitly.
This approximation, developed in a series of works <cit.>, neglects only indirect quark contributions, but the obvious advantage is that neither gluonic vertices are required nor two-loop contributions need to be calculated,
both of which are quantitatively relevant <cit.>.
This particular truncation was already used for calculations of the N_f=2 <cit.>, N_f=2+1 <cit.>, and N_f=2+1+1 <cit.> QCD phase diagrams, of baryon number fluctuations <cit.> and of mesons at nonvanishing chemical potential <cit.>.





The quark propagator DSE reads


    S^-1(p)   =Z_2 S_0^-1(p) +Z_1F C_F g^2_q γ_μ S(q) Γ_ν(p-q;-p,q) D_μν(p-q)


with the momenta p=(p⃗,p_4) and q=(q⃗,q_4).
The quark and gluon propagators have Matsubara frequencies p_4=(2n_p+1)π T and p_4=2n_pπ T, respectively, with n_p∈ℤ.
The integration and Matsubara summation are abbreviated as _q ≡ T ∑_n_q ∈ ℤ∫ d^3 q⃗  /   (2π)^3.
C_F=(N_c^2-1)/(2N_c) is the quadratic Casimir of the fundamental representation.
S(p) is the full quark propagator parametrized as

    S^-1(p)    = i p⃗γ⃗ A(p) + i(p_4+i μ) γ_4 C(p)   + B(p) +  i p_4 γ_4p⃗γ⃗ D(p).

In our calculations, we set D(p)=0.
At vanishing density, we tested its quantitative influence and found at most a change of the order of 0.0001 on the other dressing functions <cit.>.
S_0^-1(p)=i p⃗γ⃗ + i (p_4+i μ) γ_4 + Z_m m is the inverse bare quark propagator, where Z_m and Z_2 are the quark mass and wave function renormalization constants, respectively.
m is the renormalized current quark mass at the renormalization point.
Finally, Z_1F is the quark-gluon vertex renormalization constant for which we use Z_1F=Z_2 Z_1/Z_3.
Z_1 is the renormalization constant of the ghost-gluon vertex which is finite in Landau gauge <cit.> and we choose it as 1.

The quark-gluon vertex is approximated by the following model <cit.>:

    Γ_ν(q;p,l) =Z_3 γ_μΓ_mod(x)      ×(A(p^2) + A(l^2)/2δ_μ, i+C(p^2) + C(l^2)/2δ_μ, 4),
       Γ_mod(x) = d_1/(x+d_2)     + x/Λ^2 + x(α(μ)β_0/4 πln(x/Λ^2 + 1))^2 δ.

p and l are the antiquark and quark momenta, respectively, and q is the gluon momentum.
x depends on the equation in which the vertex model is used.
This is necessary to maintain multiplicative renormalizability with this model <cit.>.
x is (p^2+l^2) in the gluon propagator DSE and q^2 in the quark propagator DSE.
The function Γ_mod(x) contains the logarithmic running in the perturbative regime.
Its strength in the IR is determined by the value of the parameter d_1.
Λ and α(μ) are fit parameters taken from the gluon propagator fit; see below.
The parameter d_2 is fixed at 0.5 GeV^2.
δ is the one-loop anomalous dimension of the ghost propagator given by =-9N_c/(44N_c-8N_f) and β_0=(11N_c - 2N_f)/3 is the lowest coefficient of the QCD beta function.
The model has twice the anomalous dimension of the quark-gluon vertex which is necessary to obtain the correct anomalous dimensions of the propagators <cit.> without including higher loop contributions <cit.>.
An implicit temperature and chemical potential dependence enters via the quark propagator dressing functions A(p^2) and C(p^2).

The gluon propagator DSE reads


    D_μν^-1(p)=[D_μν^YM(p)]^-1-g^2/2 N_f Z_1F_q  [γ_μ  S(p+q)  Γ_ν(-p;-q,p+q)   S(q)].


D_μν^YM contains all gluonic diagrams.
In the medium, the propagator splits into parts transverse and longitudinal with respect to the heat bath,

    D_μν(p)=   D_L,μν(p)+D_T,μν(p)
     =   P^L_μν(p)Z_L(p^2)/p^2 + P^T_μν(p)Z_T(p^2)/p^2,

with

    P^T_μν(p)   =(1-δ_μ4)(1-δ_ν4)(δ_μν-p_μ p_ν/p⃗^2),
    
     P^L_μν(p)   =P_μν-P^T_μν(p),
    
     P_μν   =δ_μν-p_μ p_ν/p^2.

For the gluonic part, the dressing functions Z^YM_T and Z^YM_L are fitted by <cit.>

    Z^YM_ T/L (p^2)= x / (x+1)^2 ( (  c/Λ^2 / x+a_ T/L ) ^ b_ T/L +
     x( α (μ )β _ 0 / 4πln(x+1) ) ^γ)

with x=p^2/Λ^2.
We only fit the lowest Matsubara frequency and use for higher Matsubara frequencies Z^YM_T/L(p⃗^2, p_4^2)=Z^YM_T/L(p⃗^2+p_4^2,0), which is a good approximation according to lattice results <cit.>.
The parameters are c = 11.5 GeV^2 and Λ = 1.4 GeV.
α(μ)=g^2/4π=0.3 is used throughout all calculations.
The fits also set the scale in our calculations.
γ=(-13N_c+4N_f)/(22N_c-4N_f) is the one-loop anomalous dimension of the gluon propagator.

Temperature dependence enters via the fit parameters a_T/L and b_T/L.
They are fitted to the lattice results of Refs. <cit.>.
To have a smooth behavior and access to all temperatures, the parameters are fitted in terms of t = T/T_c <cit.>,

    a_T    = {[    0.46 ≥ t  :  1.41 + 0.43 t; 1 ≥ t ≥ 0.46 :  1.52 + 0.20 t;        t ≥ 1 :   3.60 - 1.88t ]. 
    
    b_T   ={[    0.49 ≥ t  :  2.20 + 0.07 t; 1 ≥ t ≥ 0.49 :  2.43 - 0.40 t;        t ≥ 1 :   2.32 - 0.29t ]. 
    
    a_L   ={[              0.53 ≥ t  :            1.41 - 2.09 t;           1 ≥ t ≥ 0.53 : 0.89 - 1.51 t + 0.77 t^2;                  t ≥ 1 :            -8.16 + 8.31t ]. 
    
    b_L   ={[                0.52 ≥ t  :               2.20 - 1.82t;             1 ≥ t ≥ 0.52 : 1.22 + 0.10 t - 0.05 t ^ 2;                    t ≥ 1 :              -1.48 + 2.75t ].







§ RESULTS



In this section we present our results for various quantities.
One of the main findings is that a straightforward extension of the μ=0 setup does not lead to a satisfactory comparison with lattice results at nonvanishing chemical potential.
We argue that the quark-gluon vertex is to be blamed for that and present calculations that corroborate this.

For our calculations several setups for the parameters are used which are listed in tab:parameters.
We explain the specific choices below where they are used for the first time.
Results along the chemical potential axis were calculated at a temperature of 47 MeV.







 §.§ Gluon propagator





The truncation employed in this work takes into account gluonic effects by unquenching the gluon propagator explicitly.
The effect of finite density on the gluon propagator can thus be studied.
We start with setup I of tab:parameters where the IR strength of the vertex d_1 is fixed such that we get close to the zero chemical potential gluon propagators observed in Ref. <cit.>.
The pion mass in the lattice calculations was quite heavy with m_π=717 MeV.
With the Gell-Mann–Oakes–Renner relation we can at least estimate the corresponding quark mass for which we use 188 MeV.

When solving the gluon propagator DSE, a renormalization of quadratic divergences is needed if a hard UV cutoff is employed.
Very often, a generalization of the Brown-Pennington projector <cit.> is employed.
Here, in order to compare with lattice results, for setup I we use a mass counter term C_sub/p^2 which is fixed by a second renormalization condition <cit.>.
Various other methods exist, but they have mostly been employed only in vacuum calculations see, e.g., <cit.>.

The resulting longitudinal gluon propagators are shown in fig:glP.
As one can see in the comparison to lattice results, the agreement down to 1 GeV is satisfactory.

Interesting quantities to monitor the reaction of the propagators on external parameters like temperature or chemical potential are the longitudinal and transverse screening masses, m_L and m_T, respectively.
They are defined as

    m_L/T=1/√(D_L/T(0)),

where D_L/T(p^2) is the scalar part of the corresponding propagators.
It should be noted that these screening masses are in general not the same as pole masses, which not even need to exist.
Of particular interest is the longitudinal screening mass which shows a distinct temperature dependence at zero chemical potential <cit.>.
With the truncation employed here, the transverse screening mass is completely determined by eq:ZTL and independent of μ, because the quark-loop contribution vanishes at zero momentum in this case.
The transverse screening mass is shown in fig:screen_mass_T.
The constant transverse screening mass appears to be an acceptable approximation when compared with the lattice data <cit.> at least for chemical potential below 800 MeV.
However, the observed increase for larger chemical potential depends on the lattice spacing and is thus at least affected by the discretization if not a complete lattice artifact <cit.>.
Taking the indications of the small dependence of the gluon propagator on chemical potential as working assumption offers the interesting possibility of approximating the gluon propagator as independent of chemical potential.
We will explore this possibility below.





We calculate the longitudinal screening mass in two setups.
We vary the quark masses and fix d_1 such that the chiral crossover transition at zero chemical potential is at 210 MeV.
This corresponds to setup II in tab:parameters.
The resulting longitudinal screening mass is shown in the left plot of fig:screen_mass_L.
For the lowest quark mass, m=1.2 MeV, we find a first order transition between μ=250 MeV and 300 MeV.
Increasing the quark mass moves the transition to higher values of the chemical potential and smooths it to a crossover.
However, the lattice results do not show a transition in the screening mass in this region at all.
The silverblaze point, where a first order transition occurs, is located in this interval, but it is not reflected in the screening masses <cit.>.

A physically better choice for fixing d_1 is to use a condition from the chemical potential axis.
Lattice results indicate that the longitudinal screening mass also remains largely unaffected by chemical potential.
Only at high chemical potential a small increase is seen, which, however, might be affected by the same problems of discretization artifacts and statistics as the transverse screening mass.
Thus, with setup III of tab:parameters we want to see if we can push the transition to higher chemical potential by varying the quark masses again and fixing d_1 such that the increase starts at μ=650 MeV, which is a conservative estimate for the beginning of this increase.
It should be noted that the onset of an increase is clearer for the transverse screening mass.
The quark masses are now considerably higher and consequently the setup would lead to a totally different crossover line between the hadronic regime and the quark-gluon plasma.
However, there is anyway no reason to expect that the vertex is independent of temperature and chemical potential, so having a different d_1 at low temperatures is to be expected.
Also, for real QCD, it was seen that fixing d_1 via the pion in the vacuum or via the transition temperature at μ=0 yields different values <cit.>.

The results for the longitudinal screening mass for setup III are shown in the right plot of fig:screen_mass_L.
Increasing the quark mass lowers the steepness of the increase after the transition.
In addition, the screening mass below the transition flattens.

As mentioned above, it is not settled if the rise of the screening masses is genuine.
However, from our results we see that variations of the quark-gluon vertex strength with temperature and/or chemical potential allow to describe quite different scenarios for the longitudinal screening mass while an extension of the truncation is required for the transverse screening mass.





 §.§ Chiral condensate






The chiral condensate is a standard quantity to test for the breaking of chiral symmetry and thus to distinguish corresponding phases.
It can be calculated from the quark propagator as

    <  ψψ>  = -N_c Z_2 Z_m _q[S(q)]    =-N_c Z_2 Z_m ∑_q_4∫d^3q/(2π)^34 T B(q)/A^2(q)q⃗^2 + C^2(q)q_4^2 + B^2(q).

Since this expression is UV divergent, it needs to be renormalized which is done by subtracting a quark condensate with a heavier renormalized mass m_s from a condensate with a light renormalized mass m_l which we take as the quark mass m in this work:

    Δ_l,h(μ,T) = -<  ψψ>_l + m_l/m_s<  ψψ>_h.

It is also convenient to normalize this expression by the vacuum value

    Δ(μ,T)=Δ_l,h(μ,T)/Δ_l,h(0,0).


We show results for the chiral condensate and the different quark masses of setup II in the left plot of fig:chiral_cond.
We observe the transitions at the same points and of the same types as for the longitudinal screening masses.
As we argued in Sec. <ref>, it is quite likely that the gluon propagator does not depend strongly on the chemical potential.
Thus, we also tested the impact of the dynamic gluon propagator by comparing two calculations with a dynamic gluon propagator and with the gluon propagator fixed at μ=0.
In the right plot of fig:chiral_cond, one can see that the gluon propagator indeed influences the position of the transition.
For the future it will thus be important either to improve the setup for the gluon propagator or to confirm its weak dependence on the chemical potential in a certain range so that using the propagator from μ=0 is a valid approximation.





 §.§ Quark number density




The quark number density can be calculated from the quark propagator by

    n(μ,T)=-∂Ω/∂μ=-N_c  N_f  Z_2 _q [γ_4  S(q)],

where Ω is the grand-canonical potential of QCD given by

    Ω=-T/Vlog𝒵(μ,T).

V is the volume of the system and 𝒵 the QCD partition function.
Eq. <ref> needs to be regularized which is done by subtracting a temperature independent term, see Refs. <cit.> for details:

    n_reg(μ,T)=n(μ,T) + N_c  N_f Z_2 ∫d^4q/(2π)^4[γ_4  S(q)].


Results for the quark number density are shown in fig:qnd.
Setup I was employed, viz., lattice data was used to fix the parameters.
Clearly, also in the quark number density a transition is seen at μ=700 MeV.





 §.§ Phase diagrams




Solving the present system of equations for higher temperatures, we can also investigate the crossover region and search for a critical point.
For comparison, we also show the same calculations for the gauge groups SU(3) and G_2.
To connect to previous work, we use the setup of Ref. <cit.> which is detailed in tab:phase_diags.
The setup for SU(2) corresponds to setup IIa of tab:parameters.
We recall that the value of the coupling is different for G_2 due to the choice of the gluon propagator; see Ref. <cit.> for details.
First results for SU(2) and G_2 were already shown in Ref. <cit.>.





For the confinement/deconfinement transitions the dual condensate <cit.> is used, which is an order parameter for center symmetry in quenched QCD.
It is computed by introducing generalized U(1) valued boundary conditions for the quarks ψ(x, 1/T) = e^i φψ(x, 0) <cit.> and projecting out the loops with winding number 1:

    Σ=∫_0^2πdφ/2πe^-i φ<  ψψ>_φ dφ.

For the identification of the chiral transitions, the chiral condensate is used.
As definition for the transitions we use the  maxima of the following derivatives:

    χ_ch   = ∂Δ_l,h/∂ T,
    χ_dec   = ∂Σ/∂ T.


The results for the chiral and confinement/deconfinement transitions are shown in fig:phase_diags.
Clearly, all three theories show a similar behavior with a crossover from vanishing chemical potential to a critical point at 0.15 GeV<μ<0.175 GeV.
The lines of the chiral and confinement/deconfinement transitions are very close and merge at the critical points.
Their locations are given in tab:phase_diags.
For all three theories, a first order transition is found beyond the critical points.






§ SUMMARY, CONCLUSIONS AND OUTLOOK




We explored the phase diagram of QC_2D and compared various quantities with results from the lattice.
The employed truncation, which includes direct unquenching effects, uses a model for the quark-gluon interaction and fits for the (quenched) gluon propagator.
For low chemical potential, the results agree well with those from the lattice.
However, for higher chemical potential, we see deviations, which can be traced back to several sources.

One is related to the gluon propagator, for which we see a somewhat stronger dependence on chemical potential than expected from corresponding lattice results.
This, however, can easily be circumvented by using the unquenched gluon propagator from μ=0.
This modification improves the agreement to some extent but is itself insufficient.
One possible source for these deviations is the quark-gluon vertex, which contains only a minimal dependence on temperature and density via the quark dressing functions in the model.
It was observed already earlier for μ=0 that the model parameters need to be changed slightly when describing vacuum physics instead of the crossover at nonvanishing temperature.
Thus, it does not come as a surprise that a dependence on the chemical potential is missing as well.
To test this, we demonstrated successfully that a modified IR strength of the interaction leads to better agreement with lattice results.





Another reason for deviations from lattice results is that we did not include diquarks yet, which are expected to play an important role.
They can be included with the Nambu-Gor'kov formalism  <cit.> which leads to a more complex structure of the quark propagator,

    𝒮_NG(p) = ([ S^+(p) T^-(p); T^+(p) S^-(p) ]).

S^+(p) is the quark propagator from eq:quarkprop and S^-(p) is related to S^+(p) by charge conjugation.
For details of the anomalous propagators T^±(p), we refer to Ref. <cit.>.
The diquark condensate, which is an order parameter for the diquark condensation phase, can be calculated as

    <  q q >  = - Z_2 _q[γ_5 M T^-(q)],

where M = T_2 τ_2 with T_2 (τ_2) carrying the color (flavor) structure <cit.>.

To get a first idea of the influence of diquarks we calculated for a few values of the chemical potential the chiral condensates with and without diquarks.
To compare with previous results, we choose the simplest generalization of the quark-gluon vertex from eq:qglvertModel to the Nambu-Gor'kov setting by neglecting its anomalous components.
As gluon propagator we use the one fixed at μ=0.
The result for the chiral condensate is shown in fig:chiral_cond_diquark.
Not only does the position of the transition move, also the order changes: it is second order now.
The corresponding diquark condensate is shown in fig:diquark.
As expected, it starts to increase at the same point at which the chiral condensate starts to drop.
For now, the calculation with diquarks has a limited resolution and 
more work is required, but they should be included in future work.

In summary, the setup of two-color QCD presents an interesting testbed where we can compare continuum with lattice methods.
Such comparisons are useful to show the path for future extensions of truncations of functional equations.
In particular since state-of-the-art truncations are quite demanding, such hints are valuable.

We also calculated the phase diagrams for QCD and QCD with the gauge group G_2 with the same truncation and found that the truncation behaves very similarly in all three cases.
This hints at a universal behavior which can be exploited by using QCD-like theories as a guide to improve truncations for real QCD.

One of the main results of this work is that the interaction between quarks and gluons as described by the quark-gluon vertex needs to be refined for large densities (and low temperatures) in order to achieve better precision.
However, the vertex is still quite an elusive object and information about its behavior beyond the vacuum is scarce.
First steps toward such a calculation were made only recently <cit.>. 
On the other hand, we can directly infer from the lattice results that the gluon propagator is most likely not what we need to be concerned with: its dependence on chemical potential is rather small.
This is reassuring, because a quantitative calculation of the propagator requires quite an elaborate truncation.
Thus, at the present level of precision, it is a convenient workaround to employ the gluon propagator from lattice calculations.






§ ACKNOWLEDGMENTS

We thank Axel Maas and Ouraman Hajizadeh for useful discussions and for providing lattice data.
We are grateful to Christian Fischer and Philipp Isserstedt for helpful discussions.
HPC Clusters at the University of Graz were used for the numerical computations.
The software programs and packages Mathematica <cit.>, DoFun <cit.>, and CrasyDSE <cit.> were used for deriving and solving numerically the DSEs.
Feynman diagrams were created with Jaxodraw <cit.>. 
Support by the FWF (Austrian science fund) under Contract No. P27380-N27 and through the doctoral program “Hadrons in Vacuum, Nuclei and Stars”, Contract W1203-N16, is gratefully acknowledged.




utphys_mod


