Baryon effects on the location of QCD's critical end point
Gernot Eichmann, Christian S. Fischer, Christian A. Welzbacher
BBaryon effects on the location of QCD’s critical end point
Gernot Eichmann, Christian S. Fischer, and Christian A. Welzbacher
Institut f¨ur Theoretische Physik, Justus-Liebig-Universit¨at Gießen,Heinrich-Buff-Ring 16, D-35392 Gießen, Germany. (Dated: October 14, 2018)The location of the critical end point of QCD has been determined in previous studies of N f = 2+1and N f = 2 + 1 + 1 dynamical quark flavors using a (truncated) set of Dyson-Schwinger equationsfor the quark and gluon propagators of Landau-gauge QCD. A source for systematic errors in thesecalculations has been the omission of terms in the quark-gluon interaction that can be parametrizedin terms of baryonic degrees of freedom. These have a potentially large dependence on chemicalpotential and therefore may affect the location of the critical end point. In this exploratory studywe estimate the effects of these contributions, both in the vacuum and at finite temperature andchemical potential. We find only a small influence of baryonic contributions on the location of thecritical end point. We estimate the robustness of this result by parameterizing further dependencieson chemical potential. I. INTRODUCTION
Heavy ion collision experiments at RHIC/BNL and thefuture CBM/FAIR facility probe the phase structure ofQCD at finite chemical potential. One of the major goalsof these experiments is the study of the existence, thelocation and the properties of a critical end point (CEP),where the chiral crossover at small chemical potentialturns into a first-order transition.From a theoretical point of view it remains unclear atthe moment whether this is indeed the case. Lattice cal-culations firmly established the crossover behavior at zerochemical potential, see e.g. [1, 2] and references therein.At finite chemical potential, lattice calculations are ham-pered by the notorious fermion sign problem. Althoughvarious extrapolation methods agree with each other atrather small chemical potential [3–11], for regions in the(
T, µ ) plane with µ B /T > N f = 2, N f = 2 + 1 and N f =2 + 1 + 1 quark flavors has been determined using the(truncated) Dyson-Schwinger equations (DSEs) for thequark and gluon propagators of QCD. A particular fo-cus in these studies has been the inclusion of back-reaction effects of the quarks onto the Yang-Mills sec-tor, which allowed to go beyond simple modeling ofthe gluon part either within the DSE approach [17, 18]or in chiral models like the Nambu-Jona-Lasinio (NJL)model [19], its Polyakov-loop extended versions [20–22]and the Polyakov-loop extended quark-meson (PQM)model [23–25]. At zero chemical potential the inclusionof the quark back-reaction on the gluons produced thecorrect temperature behavior of the quark condensateand led to predictions for the magnitude of unquench-ing effects in the gluon propagator [14], which have beenverified by subsequent lattice calculations [26]. In this approach a CEP has been found at rather largequark chemical potential ( T c , µ cq ) = (115 , µ q = 0may provide for sizable quantitative corrections. In thiswork we focus on a particular class of such corrections,namely vertex corrections that can be parametrized interms of (off-shell) baryons. In general, baryonic back-reaction effects onto the quark propagator provide a di-rect mechanism how the quark condensate may be influ-enced by changes in the baryon’s wave functions such asthe one inflicted e.g. by the nuclear liquid-gas transitionat very small temperatures [12]. These back-reaction ef-fects, however, may very well decrease in size for growingtemperatures and it is an open question whether they arestill important in the region of the QCD phase diagramwhere the putative CEP for the chiral phase transitionis located. In a two-color version of QCD this influencehas been studied in Refs. [27, 28] and found to be crucialto an extent that not only the location but even the veryexistence of a CEP is affected. Whether this is also thecase in the SU (3) theory is an open question that needsto be addressed. We regard the study reported in thiswork as a first step in this direction.The paper is organized as follows. In Sec. II we re-view our truncation scheme of the DSEs for the quarkand gluon propagators. We explain the details of thecorresponding equation for the quark-gluon vertex andidentify the diagrams that can be parametrized in termsof (non-elementary, i.e. composite) hadronic degrees offreedom. We specify an approximation scheme for theseterms that is suitable for an exploratory calculation of itseffects on dynamical chiral symmetry breaking. For sim-plicity we restrict ourselves to the case of QCD ( N c = 3)with two degenerate fermion flavors, N f = 2, and notethat a generalization to the N f = 2 + 1 case is straight-forward but very expensive in terms of CPU time. InSec. III we present our results. We first discuss the effects a r X i v : . [ h e p - ph ] D ec of the resulting baryon loop on the quark propagator inthe vacuum and at finite temperature but zero chemicalpotential. We then present an estimate for the size of theeffects that may be expected for the CEP. We concludein Sec. IV. II. DYSON-SCHWINGER EQUATIONSA. DSEs for the propagators
In order to accommodate the notation already for itsintended purpose later in this work, we specify the quarkand gluon propagators at finite temperature T and quarkchemical potential µ q and indicate the limits T → µ q → S − ( p ) = i p · γ Z + i ˜ ω n γ Z + Z m , with wave function renormalisation Z and bare quarkmass m . The dressed inverse quark propagator S − andthe Landau-gauge gluon propagator D µν are given by S − ( p ) = i p · γ A ( p ) + i ˜ ω n γ C ( p ) + B ( p ) ,D µν ( p ) = P Tµν ( p ) Z T ( p ) p + P Lµν ( p ) Z L ( p ) p (1)with momentum p = ( ω n , p ). The Matsubara frequenciesare ω n = πT (2 n + 1) for fermions and ω n = πT n forbosons, and we use the abbreviation ˜ ω n = ω n + iµ q . Alldressing functions implicitly depend on temperature andchemical potential. The projectors P T,Lµν are transverse( T ) and longitudinal ( L ) with respect to the heat bathand given by P Tµν = (1 − δ µ ) (1 − δ ν ) (cid:18) δ µν − p µ p ν p (cid:19) ,P Lµν = P µν − P Tµν , (2)where P µν = δ µν − p µ p ν /p is the covariant transverseprojector. The limit of zero chemical potential is straight-forward; in the additional zero-temperature limit themomentum p reduces to its usual O (4) − symmetric Eu-clidean form and the wave functions of the quark prop-agator become degenerate, i.e. A ( p ) = C ( p ). Further-more, the transverse (magnetic) and longitudinal (elec-tric) dressing functions of the gluon approach the samelimit in the vacuum, i.e., Z T ( p ) = Z L ( p ) ≡ Z ( p ). Inthe medium there exists also a fourth contribution to theinverse quark propagator, which vanishes in the vacuumand is proportional to ˜ ω n γ p · γ . Due to its negligiblecontribution also at higher temperatures and chemicalpotential we do not consider it throughout this work.The DSE for the quark propagator is shown diagram-matically in Fig. 1. The pieces that need to be de-termined in order to allow for a self-consistent solutionof this equation are the fully dressed gluon propagatorand quark-gluon vertex. Model calculations [17, 18] of-ten use simple ans¨atze for the gluon propagator that do − = − + − − = + FIG. 1. The DSE for the quark propagator (top panel) andthe truncated gluon DSE for N f = 2 QCD (bottom panel).Large blobs denote dressed propagators and vertices, and thehatched circle represents the quenched (lattice) propagator. not take into account the proper temperature and fla-vor dependence of the gluon self-energy. We prefer toinclude these important effects by taking the Yang-Millssector of QCD into account and calculating the back-reaction of the quarks onto the gluon explicitly. Thisframework has been gradually evolved from the quenchedcase, N f = 0 [32, 33], to two-flavor QCD [14, 34] and re-cently to N f = 2 + 1 and N f = 2 + 1 + 1 [14, 16].Such an approach has two distinct advantages oversimple modeling. On the one hand it allows us to tracethe effects of quark masses and flavors as exposed in theColumbia plot. On the other hand, it serves to take intoaccount the effects of chemical potential on the gluonpropagator explicitly, thereby rendering results at finite µ q more reliable. Furthermore, since we have explicit ac-cess to all fundamental degrees of freedom of QCD, i.e.quark, gluon and ghost propagators, we are in a positionto determine the Polyakov loop potential at all valuesof T and µ q and thereby also to study the deconfine-ment transition. This has been exploited in Refs. [14–16]for physical quarks and in Ref. [35] for a study of thesecond-order critical surface of the deconfinement transi-tion of heavy quarks. These studies are complemented bycorresponding ones using the functional renormalizationgroup, see e.g. [36] and references therein.For the gluon DSE, shown in Fig. 1, we use the samesetup as described in detail in Ref. [14]. We use lat-tice input for the quenched propagator at different tem-peratures and determine the temperature and chemicalpotential dependent effects of the quark loop explicitlyusing the quark propagator from its DSE. We work withtwo fermion flavors, N f = 2, in the isospin limit whichallows us to use one and the same quark DSE for bothflavors. The quark loop in the gluon DSE is then sim-ply multiplied by a factor of two to accommodate forboth flavors. The essentials of this setup are collected inApp. B; more details can be found in [14] and shall notbe repeated here for brevity. = + + + + = π, . .. + + + ( . .. ) N, . ..
FIG. 2. The full, untruncated Schwinger-Dyson equation for the quark-gluon vertex [37] is shown diagrammatically in the firstequation. The second equation describes the first terms of an expansion in terms of hadronic and non-hadronic contributionsto the quark-antiquark scattering kernel. In both equations, all internal propagators are fully dressed. Internal dashed lineswith arrows correspond to ghost propagators, curly lines to gluons and full lines to quark propagators. In the second equation,the dotted line describes mesons and the triple line baryons. − = − + NB + B FIG. 3. DSE for the quark propagator, where the quark-gluonvertex is separated into non-baryon and baryon contributions.
B. Baryon effects in the quark DSE andquark-gluon vertex
Let us now focus on the Dyson-Schwinger equation forthe quark-gluon vertex and explicate its structure. Thefull equation is shown in the left part of Fig. 2. It con-tains three one-loop diagrams with a fully dressed quark(solid), ghost (dashed) and gluon line (curly) runningthrough the loop and attached to the external gluon bya corresponding bare vertex. The remaining graph is agluonic two-loop diagram with a bare four-gluon vertex.These diagrams contain four- and five-point Green func-tions that are 1PI with respect to the external legs in the t channel, i.e., they contain no contributions from inter-mediate annihilation of the external quarks into a singlegluon line. The four- and five-point functions can beexpanded in skeleton diagrams with fully dressed inter-nal propagators and primitively divergent vertices [37].For our purposes we concentrate on the first non-trivialdiagram that contains a four-quark amplitude. In itsskeleton expansion, a part of the resulting diagrams canbe re-expressed in terms of Bethe-Salpeter vertices andpropagators of mesons as well as Faddeev-type verticesfor baryons.The result of such an expansion is shown in the sec-ond equation in Fig. 2. The first diagram correspondsto (off-shell) meson exchange between the quark lines. Itis important to note that this meson is not introducedas a new elementary field; it is rather a composite ob-ject of a quark and an antiquark that is described (atleast on-shell) by its Bethe-Salpeter equation (BSE). Thefirst baryon exchange diagram shows up as a two-loopdiagram involving the baryons’s Faddeev amplitude. InRef. [38] the corresponding diagram has been displayedin the quark-diquark approximation already on the levelof the vertex DSE; below, we will introduce this approxi-mation on the level of the quark DSE. Finally, we show a representative non-resonant contribution due to dressedone-gluon exchange. Note that double-counting is triv-ially avoided in this combined expansion in elementaryand effective degrees of freedom due to different quan-tum numbers in the exchange channel.The computation of the hadronic diagrams in Fig. 2 israther involved. The meson-exchange diagram requiresthe solution of a coupled system of the DSE for thequark propagator and a corresponding BSE for the me-son Bethe-Salpeter amplitude. The effect of the pionback-reaction onto the quark in the vacuum has alreadybeen explored to some extent in the context of pion-cloud contributions to light mesons and baryons [38–40]. Even more complex is the baryon-exchange diagram,which involves the computation of a Faddeev-type equa-tion for the baryon bound state. In the present workwe are interested in hadronic effects at finite chemicalpotential, which will primarily show up in the diagramincluding baryons because all elements of this diagram(the quark and baryon propagators as well as the baryonwave function) depend on chemical potential. Since thepion-exchange diagram is the hadronic contribution withminimal explicit dependence upon chemical potential werelegate its explicit study to future work and focus ex-clusively on the baryonic diagram.The resulting quark DSE with an explicit separationinto non-baryon and baryon parts is shown in Fig. 3. Thethree-loop diagram contained therein is hard to evaluatenumerically, especially at finite temperature and chemi-cal potential. To make this diagram tractable, we there-fore introduce an additional approximation and convertthe three-quark Faddeev amplitude into a quark-diquarkBethe-Salpeter amplitude for the baryon. To this end,note that the gluon in the baryon-exchange diagram isattached to the quark on the left, but since it couples toboth quarks symmetrically we can rearrange the diagramas illustrated by the first equation of Fig. 4. The incom-ing and outgoing baryon lines have to be connected by abaryon propagator which is indicated by the open circles.Inserting a separable quark-diquark ansatz for eachthree-body Faddeev amplitude leads to the second equal-ity in Fig. 4. The resulting topologies can be groupedinto two different classes: one where the incoming quarkcouples to a diquark amplitude and one where it cou- +== +++ ++ ++ )( FIG. 4. Baryonic diagram in the quark DSE with the quark-diquark approximation of the baryon’s Faddeev amplitude. + dq + dqN − = − + N B
FIG. 5. Quark DSE with diquark and baryon loop. In these loops the right vertices (circles) are Bethe-Salpeter amplitudeswhereas the left vertices (hatched circles) are effective ones summing all effects from the diagrams in Fig. 4. In the main textwe give arguments why these vertices are well approximated by bare ones. ples to a quark-diquark amplitude. The hatched ampli-tudes are effective and absorb all the remaining objectsin these graphs. As a consequence, the quark DSE takesthe form shown in Fig. 5, which contains a quark-diquarkand a baryon-diquark loop. For brevity we will refer tothem as ‘diquark’ and ‘baryon’ loops in what follows.In both diagrams the vertex appearing on the right isa proper Bethe-Salpeter amplitude, once for a diquarkand once for a baryon in quark-diquark approximation.The hatched vertices on the left carry the same quantumnumbers as their counterparts on the right but representeffective vertices that absorb all effects appearing in themulti-loop diagrams in Fig. 4. Thus, the quark DSE inFig. 5 follows directly from the original equation in Fig. 3if the quark-diquark ansatz for the Faddeev amplitude ismade. We will specify its ingredients in Sec. II C. Fornow, note that we work in the N f = 2 theory in theisospin symmetric limit, which leaves us with isospin sin-glet scalar and isospin-triplet axial-vector diquarks anda degenerate isospin doublet of nucleons.For the non-baryonic part of the quark-gluon vertex(denoted by ’NB’ in Fig. 3) we employ a construction us- In principle there are two further diagrams with a closed quarkloop where the gluon couples to each of the quarks in the loop.However, their contributions cancel each other, i.e., the ‘diquark-gluon vertex’ is zero. ing the first term of the Ball-Chiu vertex that satisfies theAbelian Ward-Takahashi identity [41], multiplied with aninfrared-enhanced function of quark and gluon momentathat accounts for the non-Abelian dressing effects andthe correct ultraviolet running of the vertex. The explicitexpressions are collected in appendix A. In the followingsubsection we complete the discussion of the quark DSEwith the last remaining ingredient, the diquark Bethe-Salpeter amplitudes together with the quark-diquark am-plitude for the baryon.
C. Diquark and baryon amplitudes
In principle, the baryon is a three-quark state anda comprehensive, full treatment of its structure shouldtake this explicitly into account. Indeed, the corre-sponding three-body Faddeev equation has been solvedin Refs. [40, 42–46] and electromagnetic as well as ax-ial form factors have been extracted [46–48]. Owing tothe dynamical formation of diquark correlations insidethe nucleon, a potentially satisfying approximation to thethree-body framework is a description in terms of quarkand diquark degrees of freedom. The BSE for such abaryon with quark and diquark constituents is displayedin Fig. 6. In this approximation, the quark and diquarkinside the nucleon interact via quark exchange and thecorresponding diquark amplitude has to be determinedfrom a separate BSE.Using a simple model for the underlying quark anddiquark propagators and ans¨atze for the diquark am-plitudes, baryon properties in the quark-diquark picturehave been determined in many works, see e.g. [49–52]and references therein. A more fundamental approach isthe use of an underlying quark-gluon interaction, fromwhich all components of such a calculation, the quarkpropagator in the complex momentum plane, the di-quark amplitude from its BSE, the diquark propagatorfrom its scattering equation, and the baryon, are deter-mined consistently without any introduction of furtherparameters. This has been performed in [53–56] usinga well-established rainbow-ladder interaction kernel forthe quark-gluon interaction. In all these calculations itturned out that a satisfactory description of the ground-state properties of the nucleon and ∆ baryon can be ob-tained using scalar and axial-vector diquarks only. Thequark-diquark approximation works well at zero temper-ature and chemical potential; below we assume that thisis still the case at finite T and µ q . Whether that is trueremains to be studied in future work.In fact, if one is only interested in the gross proper-ties of the nucleon even the influence of the axial-vectordiquark may be omitted and both the diquark and thenucleon can be represented by their leading tensor struc-ture. In this approximation, the diquark and nucleonBethe-Salpeter amplitudes are parametrized byΓ dq ( q, P ) = f dq ( q ) γ C ⊗ (cid:15) ABE √ ⊗ s ab , Γ N ( q, P ) = f N ( q ) Λ + ( P ) ⊗ δ AB √ ⊗ t ae . (3)Here, q, P are the relative and total momentum of thebound states, C = γ γ is the charge-conjugation ma-trix, and Λ + the projection operator onto positive-energystates (which we omit in the baryon loop diagram be-cause its purpose is already served by the nucleon prop-agator). We use normalized color wave functions withcapital subscripts and normalized flavor wave functionswith small subscripts; s = √ ( ud † − du † ) = √ iσ withPauli matrices σ i and t = ( uu † + dd † ) = . The solu-tions for the diquark and nucleon amplitudes determinedin the rainbow-ladder framework of Refs. [54–56] are wellparametrized by f dq ( q ) = N dq (cid:18) e − α dq · x + β dq x (cid:19) ,f N ( q ) = N N (cid:18) e − α N · x + β N (1 + x ) (cid:19) (4)with x = q / Λ and the scale Λ = 0 . N dq = 15 . N N = 28 .
4. The parameters are α dq = 0 .
85 and α N = 1 . β dq = 0 .
02 and β N = 0 .
03 for the UV behavior. (cid:1869) (cid:1853)(cid:1854) (cid:1857) (cid:1842) (cid:1985)(cid:1944) (cid:1985) (cid:1960)(cid:1973) - (cid:1985) (cid:1960)(cid:1973) (cid:1985)(cid:1944) FIG. 6. The Bethe-Salpeter equation for the baryon in thequark-diquark approximation.
The remaining question concerns the effective ampli-tudes in Fig. 5 which we did not specify yet. For thosewe resort to a simple approximation: we take them as‘bare’, i.e., we use Eq. (3) with f dq ( q ) = f N ( q ) = 1.This is analogous to the treatment of the pion loop inthe quark DSE in Refs. [38–40]. It can be motivated byestimating the overall strength of the baryon diagram inFigs. (3–4) from its contribution to the quark conden-sate. When connecting the quark lines with the scalar qq vertex (calculated from its inhomogeneous BSE), the re-sulting vacuum bubble C is proportional to the integrated(off-shell) scalar form factor of the nucleon: C = 2 m N (cid:90) d P (2 π ) g S ( P , Q = 0) P + m N . (5)This can be seen by inserting the covariant Faddeev equa-tion for the three-body amplitude in the first line ofFig. 4; the resulting quantity is what appears in baryonform factor diagrams such as in Ref. [46]. The on-shellvalue of the scalar form factor at P = − m N is de-termined by the nucleon sigma term via the Feynman-Hellmann theorem: σ N = m q g S ( − m N ,
0) = m q dm N dm q ≈ m π dm N dm π . (6)The magnitude of C obtained with the experimentalvalue σ N = 45 MeV, together with an integral cutoffat m N = 0 .
94 GeV, is similar to the value obtained fromnumerically tracing the sum of the diquark and baryondiagrams in Fig. 5 with a scalar vertex, however withonly one amplitude dressed in each case. Dressing bothwould overestimate the strength by far due to the nor-malization factors in Eq. (4). Since we are only interestedin the gross effects of baryons on the phase diagram, wetherefore view this as a justified approximation.
D. Quark DSE including hadronic loops
Putting everything together, we will now give explicitexpressions for the diagrams in the quark and gluon DSEs We also calculated the vacuum bubble from the baryon diagramin Fig. 3 directly using a reasonable off-shell ansatz for the three-quark Faddeev amplitude instead of a cutoff; the result is in thesame ballpark. including the hadronic back-reaction diagrams. As willbecome clear below, the hadronic effects on the quarkpropagator enter on the level of ten-percent corrections.From the diagrammatic form of the quark-gluon vertexDSE this is exactly the order of magnitude as expected,since the corresponding diagrams are suppressed by afactor 1 /N c . In the Yang-Mills sector of QCD the to-tal quark effects are on the level of a 1 /N c correction.Therefore, hadronic contributions to the quark-loop dia-gram in the gluon DSE only contribute at 1 /N c and it iswell justified to neglect those in a first exploratory cal-culation. Hence we will use the same truncation for thegluon DSE as in previous works [14, 16, 34].In the quark DSE we take into account the three dia-grams in Fig. 5. If we denote the quark dressing functionsin Eq. (1) collectively by H ( p ) = A ( p ) , B ( p ) , C ( p ) andabbreviate the gluon, diquark and baryon-loop contribu-tions to the quark self-energy by Σ glueH , Σ dqH and Σ baH , theresulting equations read H ( p ) = Z λ H + Σ glueH + Σ dqH + Σ baH , (7)where λ B = m is the bare current-quark mass, λ A = λ C = 1, and Z is the quark wave-function renormal-ization constant. The gluon-dressing loop Σ glueH containsthe unquenched, temperature- and chemical-potential de-pendent gluon propagator together with a model for thequark-gluon vertex [14, 16, 34]; the explicit formulas arerelegated to App. A. The self-energy contributions fromthe diquark and baryon loop are given byΣ dqH ( p ) = 12 (cid:88)(cid:90) q f dq ( q − p ) D dq ( q + p ) q A ( q ) + ˜ ω q C ( q ) + B ( q ) K dqH , Σ baH ( p ) = 13 (cid:88)(cid:90) q f N ( q − p ) D dq ( q − p ) q + ( ω q + 3 iµ q ) + m N K baH (8)with K dqA = p · qp A ( q ) , K baA = p · qp ,K dqC = ˜ ω q ˜ ω p C ( q ) , K baC = ω q + 3 iµ q ˜ ω p , (9) K dqB = − B ( q ) , K baB = − m N . We have already carried out all color and flavor traces.At finite temperature and chemical potential the ar-guments p , q serve as abbreviations for p = ( ω n , p ), q = ( ω m , q ). The (fermionic) Matsubara frequencies aregiven by ω n = πT (2 n + 1) and we write ˜ ω = ω + iµ q with quark chemical potential µ q . The Matsubara sumas well as the integration over the loop three-momentum q is abbreviated by (cid:80)(cid:82) q = T (cid:80) n q (cid:82) d q (2 π ) , and the diquarkpropagator D dq is given by D dq ( q ± p ) = 1( q ± p ) + ( ω q ± ω p + 2 iµ q ) + m dq . (10) In these expressions we take into account the lowest-lying J P = / baryon multiplet for the two-flavorcase, i.e., the nucleon, in the approximation with scalardiquarks only. In principle other baryons may alsocontribute but since they are suppressed by powers of m N /m B with respect to the nucleon their influence is cer-tainly subleading. There is, however, one exception: theparity partner of the nucleon becomes (approximately)mass-degenerate once chiral symmetry is restored, i.e.in the high temperature/density phase. Performing theDirac traces of the corresponding loops, it turns out thatcontributions from mass-degenerate multiplets of paritypartners cancel each other in Σ baB , while they add upin the other two contributions to the quark self-energy.We take this effect qualitatively into account by mul-tiplying the right-hand side of Σ baB with an additionalfactor M ( T, µ q ) /M (0 ,
0) evaluated at zero momentumand lowest Matsubara frequency. Here M = B/A de-fines the renormalization-point independent quark massfunction in the medium. This factor has no effect inthe vacuum but mimics the cancellation of multipletsof parity partners in the chirally restored phase, where M ( T, µ q ) becomes small. A corresponding factor of2 − M ( T, µ q ) /M (0 ,
0) is added to Σ baA and Σ baC . Note thatthe diquark loop, which has been derived from the orig-inal baryon three-body diagram via the quark-diquarkpicture of baryons, contains diquarks only; thus thereare no damping/enhancement factors in Σ dqA,B,C .The remaining unknowns in Eq. (8) are the temper-ature and chemical-potential dependence of the diquarkand baryon masses and Bethe-Salpeter amplitudes. Ide-ally these need to be determined consistently from theirBSEs evaluated at finite T and µ c . This formidable nu-merical task is yet to be performed and relegated to fu-ture work. Here, in this exploratory work, we resort tothe vacuum expressions for the diquark and nucleon am-plitudes as given in Eq. (4), evaluated at four-momentathat include temperature effects in the form of Matsub-ara frequencies and the results for the masses from thecorresponding bound state calculations m N = 0 .
938 GeVand m dq = 0 .
810 GeV. Certainly this can only be a firstapproximation on a qualitative level. In order to gaugethe quantitative effects of including potential changes ofthe baryon and diquark masses and wave functions withchemical potential, we will introduce and discuss addi-tional dependencies on µ q in section III C. III. RESULTSA. Vacuum
Before we consider baryon effects at finite tempera-ture and chemical potential, we first study the impactof the different loops on the strength of dynamical chiralsymmetry breaking in the vacuum. Both diquark andbaryon loops originate from diagrams in the quark-gluonvertex DSE that contain additional quark loops. In gen- -5 -4 -3 -2 -1 p [GeV ] M ( p ) [ G e V ] gluon-dressing loop only+ diquark loop+ baryon loop+ baryon loop + diquark loop 100 150 200 250 T [MeV] ∆ l . h . ( T ) / ∆ l . h . ( ) gluon-dressing loop only+ diquark loop+ baryon loop+ baryon loop + diquark loop+ baryon loop + diquark loop (strength rescaled) gluon dressing only+ diquark + baryon+ baryon + diquark gluon dressing only+ diquark + baryon+ baryon + diquark+ baryon + diquark (rescaled strength)
100 150 200 250
T [MeV] ∆ l . h . ( T ) / ∆ l . h . ( ) gluon-dressing loop only+ diquark loop+ baryon loop+ baryon loop + diquark loop+ baryon loop + diquark loop (strength rescaled) FIG. 7.
Left: quark mass function with and without diquark and baryon loops included.
Right: regularized and normalizedcondensate as a function of temperature with and without diquark and baryon loops. In addition, we display the result withrescaled strength in the quark-gluon interaction, see main text for details. eral these are chirally restoring, as has been discussedin Refs. [38, 57]. On the lattice, the unquenched quarkpropagator has indeed a smaller mass function than thequenched one [58]. While this behavior is expected forthe total sum of all unquenching contributions to thequark propagator, it is not necessarily true for each indi-vidual diagram such as those investigated here. As a firstexercise we therefore study the sign and magnitude ofthe individual effects of each contribution onto the quarkmass function. For the calculation we use two dynamicalquark flavors in the gluon DSE, N f = 2, a renormalizedbare quark mass m ( µ ) = 0 . µ = 80 GeV, and an interaction strength pa-rameter d = 8 .
05 GeV in the quark-gluon vertex (cf.App. A). These values have been taken over from the N f = 2 + 1 theory (set A in Ref. [16], matching cor-responding lattice results) without further adaption forreasons of simplicity.Our results are displayed in the left diagram of Fig. 7.It shows the quark mass function calculated with only theusual gluon-dressing loop, compared to results includingthe diquark and baryon loops individually and in combi-nation. Indeed, all additional contributions are chirallyrestoring, with a larger effect coming from the baryonloop. The total contribution of both loops is ∼ M (0) = 640 MeV to M (0) = 545 MeV, whereas theimpact is immaterial in the large momentum regime. In total, we find that the effect due to baryons is ratherlarge. Lattice QCD finds total unquenching effects in the Note that the considerable size of the generated quark mass isa direct effect of performing a N f = 2-calculation while workingwith scales adapted to the N f = 2 + 1 theory. An additionalback-coupling of the strange quark would reduce the strength ofthe gluon propagator and decrease the quark mass considerably. quark mass function of less then 20 % [58]; thus our bary-onic effects leave almost no room for other unquenchingcorrections like e.g. meson back-coupling effects. Thismay be attributed to the comparably simple approxima-tion of the baryon wave functions used in this work. Onthe other hand, one of the goals herein is to gauge the sys-tematic effects of such contributions onto the QCD phasediagram. In such a study it seems better to over- thanto underestimate the induced systematic corrections. B. Finite temperature
Next we assess the effects of the diquark and baryonloops on the quark condensate evaluated at finite T . Weuse a regularized expression for the condensate,∆ l,h = (cid:104) ¯ ψψ (cid:105) l − m l m h (cid:104) ¯ ψψ (cid:105) h , (11)which eliminates the divergences appearing for non-zerobare quark masses. The definition of (cid:104) ¯ ψψ (cid:105) is givenin Eq. (C1). For the heavy quark mass we choose m h (80 GeV) = 100 MeV. In order to evaluate the cor-responding condensate (cid:104) ¯ ψψ (cid:105) h in the N f = 2 theorywe would need to solve the complete coupled system ofDSEs, Fig. 1, a second time for each temperature andchemical potential. However, the sole purpose of (cid:104) ¯ ψψ (cid:105) h is regularization. Thus it turns out to be sufficient toevaluate this quantity from the quark DSE with modifiedquark mass m l → m h in the bare quark propagator S − ,but keeping the gluon and the quark-gluon vertex (in-cluding baryonic loops) from the light-quark calculation.We have explicitly checked that this procedure is a goodapproximation for some selected values of temperatureand chemical potential and then adopted it throughoutthe phase diagram. The transition temperatures for thechiral crossover are extracted from the maximum of thechiral susceptibility.Our results for µ q = 0 are displayed in the right dia-gram of Fig. 7. Compared to the calculation without di-quark and baryon loops, the additional loops reduce thestrength of dynamical chiral symmetry breaking and thetransition temperature for the chiral crossover reducescorrespondingly. It turns out that this sizable impact ofthe baryonic contributions on the chiral transition can bealmost completely reabsorbed into the vertex truncationby rescaling the strength of the ’NB’-part of the vertex.To this end we modify the parameter d (cf. Eq. (A2))such that the critical temperature does not change at µ q = 0 upon taking baryon loops explicitly into account.This amounts to d = 8 .
05 GeV → d = 8 .
94 GeV .The resulting condensate (the red solid curve in the plot),where baryon and diquark effects are included, recoversto very good accuracy the original shape of the conden-sate.This observation is important for our general strategy.In Ref. [16] the lattice data for the condensate of the N f = 2 + 1 theory have been reproduced point-wise in aformulation using the gluon-dressing loop only, withoutmaking the baryonic degrees of freedom explicit. Here,for N f = 2, we observe that we can reproduce a similarfunctional dependence of ∆ l,h ( T ) using explicit baryonicdegrees of freedom and a rescaled version of the quark-gluon interaction. We regard this as a strong indicationthat the same property holds in the N f = 2 + 1 theory.In the following, we therefore use the N f = 2 theory withrescaled interaction strength parameter d as a templateto study the baryonic effects at finite µ q .While at zero chemical potential all effects can be ab-sorbed into d , this is not a priori clear for the finitechemical potential case because the diquark and baryonloops contain a much stronger explicit dependence on µ q than the gluon dressing loop, as discussed above. In thenext section we will explore the consequences of these ad-ditional contributions for the location of the critical endpoint (CEP). C. Finite temperature and chemical potential
In the Dyson-Schwinger approach to the QCD phasediagram the introduction of (real) quark chemical poten-tial is straightforward; cf. Eq. (1), where the dressingfunctions of the quark propagator become complex. Theimpact of chemical potential is apparent in the quarkpropagator but it also affects the gluon (cf. App. B) aswell as the quark-gluon vertex (cf. App. A) due to theexplicit unquenching procedure. In this section we in-vestigate how this nontrivial influence is modified by thebaryon and diquark loops. We use the maximum of thechiral susceptibility, Eq. (C3), as the definition of the(pseudo-) critical temperature.In Fig. 8 we show the results in the T- µ q plane. Thesolid (black) curve is the result for the unquenched sys- µ q [MeV] T [ M e V ] gluon-dressing loop only+ baryon loop + diquark loop (rescaled strength)+ baryon loop + diquark loop (rescaled strength, prefactor) gluon dressing only+ baryon + diquark (rescaled strength)+ baryon + diquark (rescaled strength, prefactor) FIG. 8. Comparison of the phase diagram for N f =2 includingdifferent types of selfenergy contributions. tem with gluon-dressing loop and no baryonic effects. Forsmall values of the chemical potential µ q = 0 the transi-tion is a cross-over up to the filled circle, which indicatesthe 2nd order critical end point (CEP) at the criticalvalue µ q = µ cq . The two dashed lines emerging from theCEP mark the first order spinodal region for µ q > µ cq . Inthis case we find a critical endpoint at( T c , µ cq ) = (177 , N f = 2. Comparison with thedashed (red) curve, which includes explicit baryonic ef-fects with rescaled vertex strength, leads us to the fol-lowing observations: • a critical endpoint still exists; • the chiral phase transition lines are almost on topof each other; • the critical endpoint is shifted by less than 5 MeVto smaller chemical potential.The first observation is important because in QC D adisappearance of the critical endpoint was observed af-ter introducing the two-color equivalent of baryonic ef-fects [27, 28], see also [29–31]. While a negligible in-fluence of baryonic degrees of freedom on the transitionline is generally expected at low chemical potential, inour case such a behavior also persists for higher µ q andremains true for the critical endpoint. This also impliesthat our original truncation with the unquenched gluon-dressing loop can implicitly absorb baryonic effects atfinite chemical potential, at least those that are capturedby our simple approximation.Potentially important effects beyond our currentscheme are additional dependencies of the baryon’s massand wave function on the chemical potential. Unfortu-nately not much is know in this respect. In Ref. [60] theauthors investigate the thermal properties of baryons at µ q = 0 and find that the amplitudes are almost indepen-dent of T , whereas the masses rapidly increase aroundthe (pseudo-) critical temperature. As a result, baryoniccontributions would decrease in importance in a regionclose to the transition line. Such temperature-dependenteffects would presumably have no impact on our resultssince they can be reabsorbed in the strength of the ’NB’-part of the vertex.In the absence of explicit knowledge, we gauge the im-pact of modifications of the baryon wave function withchemical potential by multiplying the baryon loop witha function f κ ( µ q ) = 1 − µ q / Λ κ a κ ( µ q / Λ κ ) + b κ ( µ q / Λ κ ) , (13)where we make use of recent evaluations of the curvatureof the chiral transition line on the lattice [9–11]. Thiscurvature can be parametrized in terms of a quantity κ , T ( µ B ) T (0) = 1 − κ (cid:18) µ B T ( µ B ) (cid:19) , (14)that characterizes the lowest order in a Taylor expan-sion in the baryonic chemical potential µ B . While recentlattice values for κ range between 0 . . . . . /N f corrections,we adopt the N f = 2 + 1 value κ = 0 . N f = 2-calculation and match the coefficientsΛ κ , a κ , b κ in Eq. (13) such that we reproduce the lat-tice curvature in a region where the lattice can well betrusted. This is possible for Λ κ = 0 .
714 GeV, a κ = − . b κ = 36. The function f κ then has a minimumat µ q ≈
120 MeV. For larger chemical potential we use f κ ( µ q ) = f κ ( µ q = 120 MeV) to make the function mono-tonic. The resulting phase diagram is shown as the dash-dotted (indigo) curve in Fig. 8, with a new location of thecritical end point at( T c , µ cq ) = (197 , . (15)Due to the smaller curvature at low chemical potential,the CEP shifts by ∼
10% towards larger temperatureand ∼
5% towards smaller chemical potential. The ra-tio µ cB /T c changes accordingly from µ cB /T c = 2 . µ cB /T c = 1 .
9. These changes are by no means dramaticbut quantitatively significant. However, this comes at theexpense of a modification of the strength of the baryonloop by more than 50% due to the additional function f κ . Whether such a variation of the baryon wave func-tion and masses with chemical potential is realistic or notneeds to be investigated in the future. It is also by nomeans clear whether baryonic effects are the only possiblesource of a smaller curvature in the present framework.Here we only demonstrated that it is possible in principlethat baryonic corrections can induce such an effect.Qualitatively, it is always the case that the CEP shiftstowards larger temperatures and smaller chemical poten-tial if f κ ( µ q ) is smaller than one for µ q >
0. The opposite effect can be obtained if f κ ( µ q ) is chosen to be larger thanone: the CEP then shifts towards smaller temperaturesand larger chemical potential. Apart from arguments bycomparison with the lattice we see no physical reasona priori why baryon effects should have one effect or theother. Again, this needs to be studied in a more advancedframework. IV. SUMMARY AND CONCLUSIONS
In this exploratory study we extended our existingtruncation of the coupled system for the quark and gluonDyson-Schwinger equations to take explicit baryonic de-grees of freedom into account. This was achieved byconsidering a specific class of diagrams in the Dyson-Schwinger equation for the quark-gluon vertex wheregenuine hadronic contributions can be identified. Uponintroducing baryons through the quark-diquark picture,the baryon diagrams enter as a baryon-diquark loop and aquark-diquark loop in the quark Dyson-Schwinger equa-tion in addition to the gluon-dressing loop. In the N f = 2calculation performed herein we employ the vacuum am-plitudes and masses for the nucleon and the (scalar) di-quark, but we take into account cancellation effects dueto the degeneration of the chiral partner of the nucleonthrough a factor that couples the baryon-diquark loopto the chiral dynamics of the system. With this setupwe performed a calculation of the QCD phase diagramand find that the inclusion of baryon degrees of freedomchanges the location of the critical endpoint only by afew MeV in T and µ q .More drastic effects are possible once a dependence ofthe baryon masses and wave functions on chemical po-tential are taken into account. We estimated these usinga parametrization that reproduces the lattice transitionline at small chemical potential. As a result we find a shiftof the CEP in the 5 . . .
10% range which drives the ratio µ cB /T c slightly below 2. We expect that these results ob-tained in the two-flavor theory still hold qualitatively forthe N f = 2+1 case. This will be explored in future work,where we also strive to determine the chemical potentialdependence of the baryon and diquark masses and wavefunctions explicitly. ACKNOWLEDGEMENTS
We thank Bernd-Jochen Schaefer and Lorenz vonSmekal for fruitful discussions. Furthermore we thankChristian H. Lang for contributions in the early stagesof this project. This work has been supported by theHelmholtz International Center for FAIR within theLOEWE program of the State of Hesse and by the Ger-man Science Foundation DFG under project number TR-16.0
Appendix A: Gluon contribution to quarkself-energy
To make the paper self-contained we collect in this ap-pendix the ingredients of the non-baryonic (NB) part ofthe quark DSE, i.e., the gluon-dressing loop in Fig. 5.Our ansatz for the quark-gluon vertex that appearstherein is given byΓ µ ( p, q ; k ) = γ µ ( δ µi Γ S + δ µ Γ ) Γ( k ) , (A1)whereΓ S = A ( p ) + A ( q )2 , Γ = C ( p ) + C ( q )2 , Γ( k ) = d d + k + xx + 1 (cid:20) β α π ln(1 + x ) (cid:21) δ (A2)with x = k / Λ . Here, p = ( ω p , p ) and q = ( ω q , q )are the fermion momenta and k = ( ω k , k ) is the gluonmomentum. The scales d = 0 . and Λ = 1 . α = 0 . d and Λ controlthe renormalization-group running of the vertex functionfrom the large- into the low-momentum region, d con-trols the strength of the quark-gluon interaction at smallmomenta and therefore the amount of dynamical chiralsymmetry breaking in the hadronic phase. Its value isdiscussed in the main text. The ultraviolet momentumregion is governed by the anomalous dimension of thevertex δ = − Nc N c − N f and β = N c − N f .The explicit expressions for the gluonic parts of the selfenergy are:Σ glueA = Z C F g (cid:88)(cid:90) q Γ( k ) D ( q ) A ( q ) K AA + C ( q ) K AC p , Σ glueB = Z C F g (cid:88)(cid:90) q Γ( k ) D ( q ) B ( q ) K BB , (A3)Σ glueC = Z C F g (cid:88)(cid:90) q Γ( k ) D ( q ) A ( q ) K CA + C ( q ) K CC ˜ ω p , where k = p − q , C F = is the Casimir operator, Z is thewave-function renormalization constant, and Γ has beendefined above. The Matsubara sum as well as the inte-gration over the loop three-momentum q is representedby (cid:80)(cid:82) q = T (cid:80) n q (cid:82) d q (2 π ) . The denominator of the quarkpropagator is given by D ( q ) = q A ( q )+ ˜ ω q C ( q )+ B ( q ) and the kernels K read K AA = Γ S (cid:20) Z L k ω k k (cid:18) p · q − p · k q · kk (cid:19) + 2 Z T k p · k q · kk (cid:21) + Γ Z L k k k p · q ,K AC = (Γ S + Γ ) Z L k p · k k ˜ ω q ω k ,K BB = Γ S (cid:18) Z T k + Z L k ω k k (cid:19) + Γ Z L k k k ,K CA = (Γ S + Γ ) Z L k q · k k ω k ,K CC = Γ S (cid:18) Z T k + Z L k ω k k (cid:19) ˜ ω q − Γ Z L k k k ˜ ω q . (A4) Appendix B: Unquenching the gluon
In the Yang-Mills sector we solve the gluon DSE in-cluding the quark-loop contributions shown in Fig. 1, i.e. D − µν ( k ) = (cid:2) D que.µν ( k ) (cid:3) − − N f (cid:88) f Π fµν ( k ) , (B1)Π fµν ( k ) = g Z f (cid:88)(cid:90) p Tr (cid:2) γ µ S f ( p ) Γ fν ( p, q ; k ) S f ( q ) (cid:3) , with the explicit flavor dependence indicated by the su-perscript f . D que.µν ( k ) denotes the quenched gluon propa-gator which is taken from lattice calculations [59]; the fit-ting procedure has been discussed in Ref. [33]. Γ fν ( p, q ; k )is given in Eq. (A1) but Γ is evaluated for p + q instead of k to ensure multiplicative renormalizability(cf. Eq. (B4) below). The lattice fits for the gluon dress-ing functions Z T,L (cf. Eq. (1)) are given by Z T,L ( k ) = x ( x + 1) (cid:20) (cid:18) ˆ cx + a T,L ( T ) (cid:19) b T,L ( T ) + x (cid:18) β α π ln(1 + x ) (cid:19) γ (cid:21) , (B2)where x = k / Λ and a T,L ( T ), b T,L ( T ) are temperature-dependent fit parameters. In the ultraviolet, the log-arithmic term leads to the perturbative running withanomalous dimension γ = − N c +4 N f N c − N f . The temperature-independent parameters are ˆ c = 5 .
87 and Λ = 1 . T c = 277 MeV for thequenched SU(3) theory, the temperature-dependent pa-1rameters are given by a L ( t ) = (cid:26) . − . · t + 0 . · t if t < . · t − . t > ,a T ( t ) = (cid:26) .
595 + 1 . · t if t < . · t − . t > ,b L ( t ) = (cid:26) . − . · t + 0 . · t if t < . · t + 0 . t > ,b T ( t ) = (cid:26) .
355 + 0 . · t if t < . · t + 0 . t > t := T /T c . Note that since this expression repre-sents the quenched gluon propagator it is independentof chemical potential, which enters the gluon DSE onlythrough the quark loop, and N f = 0 in the anomalousdimension.By contracting Eq. (B1) with the projectors in Eq. (2)one arrives at the equations for the transverse and longi-tudinal parts of the quark loop. Since we are using hardmomentum cutoffs in the numerical integration, theseneed to be carefully regularized to remove quadratic di-vergencies without spoiling the Debye screening masses.This procedure is described in the appendix of Ref. [14].The resulting equations readΠ T,L ( k ,
0) = 2 g Z (cid:88)(cid:90) p Γ( p + q ) D ( p ) D ( q ) K T,L (B4)for the contribution with lowest Matsubara frequency,where p and q = p + k are the quark momenta in theloop and K T = A ( p ) A ( q ) Γ S (cid:18) p · k ) k + 2 p · k − p (cid:19) ,K L = A ( p ) A ( q ) (cid:20) Γ S (cid:18) p · k k · qk − p · q (cid:19) + Γ p · q (cid:21) + B ( p ) B ( q ) (Γ − Γ S ) − C ( p ) C ( q ) (Γ S + Γ ) ˜ ω p . The higher modes are accessed via Π
T,L ( k , ω k ) → Π T,L ( k + ω k , Appendix C: Chiral condensate and definition of T C in a crossover region The quark chiral condensate for a flavor f is given by (cid:104) ¯ ψψ (cid:105) f = Z Z m N c T (cid:88) n (cid:90) d p (2 π ) Tr D [ S f ( p )] , (C1)where Z is the quark wave function renormalizationconstant, Z m the quark mass renormalization constant, N c the number of colors and Tr D indicates the Diractrace. Due to the nonzero quark mass the condensate isquadratically divergent and has to be regularized. Weconsider two ways to define a (pseudo-) critical temper-ature that are both connected to the chiral condensate.The first one is the inflection point of the chiral conden-sate with respect to temperature, T infl.C = max ∀ T (cid:12)(cid:12)(cid:12)(cid:12) ∂ (cid:104) ¯ ψψ (cid:105) f ∂T (cid:12)(cid:12)(cid:12)(cid:12) , (C2)and the second one returns the maximum of the chiralsusceptibility: T χC = max ∀ T (cid:12)(cid:12)(cid:12)(cid:12) ∂ (cid:104) ¯ ψψ (cid:105) f ∂m f (cid:12)(cid:12)(cid:12)(cid:12) . (C3)Both definitions are independent of the regularization ofthe divergent part (finite bare masses). While they givedifferent results for T C in a crossover region, they returnthe same position for the critical end point. [1] S. Borsanyi et al. [Wuppertal-Budapest Collaboration],JHEP (2010) 073 [arXiv:1005.3508 [hep-lat]].[2] A. Bazavov, T. Bhattacharya, M. Cheng, C. DeTar,H. T. Ding, S. Gottlieb, R. Gupta and P. Hegde et al. ,Phys. Rev. D (2012) 054503 [arXiv:1111.1710 [hep-lat]].[3] Z. Fodor and S. D. Katz, JHEP , 014 (2002);[arXiv:hep-lat/0106002];[4] R. V. Gavai and S. Gupta, Phys. Rev. D , 114014(2005); [arXiv:hep-lat/0412035];[5] P. de Forcrand and O. Philipsen, Phys. Rev. Lett. (2010) 152001 [arXiv:1004.3144 [hep-lat]]. P. de For-crand and O. Philipsen, PoS LATTICE2008 (2008)208. [arXiv:0811.3858 [hep-lat]]. C. Bonati, M. D’Elia,P. de Forcrand, O. Philipsen and F. Sanfillippo, PoS LATTICE (2014) 219 [arXiv:1311.0473 [hep-lat]].[6] O. Kaczmarek, F. Karsch, E. Laermann, C. Miao,S. Mukherjee, P. Petreczky, C. Schmidt and W. Soeldner et al. , Phys. Rev. D (2011) 014504 [arXiv:1011.3130[hep-lat]].[7] G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, JHEP (2011) 001 [arXiv:1102.1356 [hep-lat]].[8] P. Cea, L. Cosmai, M. D’Elia, A. Papa and F. Sanfilippo,Phys. Rev. D (2012) 094512 [arXiv:1202.5700 [hep-lat]].[9] C. Bonati, M. DElia, M. Mariti, M. Mesiti, F. Negroand F. Sanfilippo, Phys. Rev. D (2015) 5, 054503doi:10.1103/PhysRevD.92.054503 [arXiv:1507.03571[hep-lat]]. [10] R. Bellwied, S. Borsanyi, Z. Fodor, J. Gnther,S. D. Katz, C. Ratti and K. K. Szabo, Phys. Lett.B (2015) 559 doi:10.1016/j.physletb.2015.11.011[arXiv:1507.07510 [hep-lat]].[11] P. Cea, L. Cosmai and A. Papa, arXiv:1508.07599 [hep-lat].[12] W. Weise, Prog. Part. Nucl. Phys. (2012) 299[arXiv:1201.0950 [nucl-th]].[13] K. Fukushima, J. Phys. G (2012) 013101[arXiv:1108.2939 [hep-ph]].[14] C. S. Fischer and J. Luecker, Phys. Lett. B (2013)1036 [arXiv:1206.5191 [hep-ph]].[15] C. S. Fischer, L. Fister, J. Luecker and J. M. Pawlowski,Phys. Lett. B (2014) 273 [arXiv:1306.6022 [hep-ph]].[16] C. S. Fischer, J. Luecker and C. A. Welzbacher, Phys.Rev. D (2014) 3, 034022 [arXiv:1405.4762 [hep-ph]].[17] S. x. Qin, L. Chang, H. Chen, Y. x. Liu andC. D. Roberts, Phys. Rev. Lett. (2011) 172301[arXiv:1011.2876 [nucl-th]].[18] B. Wang, Z. F. Cui, W. M. Sun and H. S. Zong, FewBody Syst. (2014) 47 [arXiv:1404.3073 [hep-ph]].[19] J. Berges and K. Rajagopal, Nucl. Phys. B (1999) 215 doi:10.1016/S0550-3213(98)00620-8 [hep-ph/9804233]; F. Xu, H. Mao, T. K. Mukherjeeand M. Huang, Phys. Rev. D (2011) 074009[arXiv:1104.0873 [hep-ph]].[20] K. Fukushima, Phys. Lett. B591 (2004) 277-284. [hep-ph/0310121].[21] E. Megias, E. Ruiz Arriola and L. L. Salcedo, Phys. Rev.D (2006) 065005 [hep-ph/0412308].[22] C. Ratti, M. A. Thaler, W. Weise, Phys. Rev. D73 (2006) 014019. [hep-ph/0506234].[23] B. -J. Schaefer, J. M. Pawlowski, J. Wambach, Phys.Rev.
D76 (2007) 074023. [arXiv:0704.3234 [hep-ph]].[24] V. Skokov, B. Stokic, B. Friman and K. Redlich, Phys.Rev. C (2010) 015206 [arXiv:1004.2665 [hep-ph]];V. Skokov, B. Friman and K. Redlich, Phys. Rev. C (2011) 054904 [arXiv:1008.4570 [hep-ph]].[25] T. K. Herbst, J. M. Pawlowski, B. -J. Schaefer, Phys.Lett. B696 (2011) 58-67. [arXiv:1008.0081 [hep-ph]].[26] R. Aouane, F. Burger, E.-M. Ilgenfritz, M. Mller-Preussker and A. Sternbeck, Phys. Rev. D (2013) 11,114502 [arXiv:1212.1102 [hep-lat]].[27] N. Strodthoff, B. J. Schaefer and L. von Smekal, Phys.Rev. D (2012) 074007 [arXiv:1112.5401 [hep-ph]].[28] N. Strodthoff and L. von Smekal, Phys. Lett. B (2014) 350 [arXiv:1306.2897 [hep-ph]].[29] S. Cotter, P. Giudice, S. Hands andJ. I. Skullerud, Phys. Rev. D (2013) 3, 034507doi:10.1103/PhysRevD.87.034507 [arXiv:1210.4496[hep-lat]].[30] T. Boz, S. Cotter, L. Fister, D. Mehta and J. I. Skullerud,Eur. Phys. J. A (2013) 87 doi:10.1140/epja/i2013-13087-6 [arXiv:1303.3223 [hep-lat]].[31] T. Boz, P. Giudice, S. Hands, J. I. Skullerud andA. G. Williams, arXiv:1502.01219 [hep-lat].[32] C. S. Fischer, Phys. Rev. Lett. (2009) 052003;[arXiv:0904.2700 [hep-ph]].[33] C. S. Fischer, A. Maas and J. A. Muller, Eur. Phys. J. C , 165 (2010) [arXiv:1003.1960 [hep-ph]].[34] C. S. Fischer, J. Luecker and J. A. Mueller, Phys. Lett.B (2011) 438 [arXiv:1104.1564 [hep-ph]]; C. S. Fis-cher and J. A. Mueller, Phys. Rev. D (2011) 054013[arXiv:1106.2700 [hep-ph]]. [35] C. S. Fischer, J. Luecker and J. M. Pawlowski, Phys. Rev.D (2015) 1, 014024 [arXiv:1409.8462 [hep-ph]].[36] N. Christiansen, M. Haas, J. M. Pawlowski andN. Strodthoff, arXiv:1411.7986 [hep-ph]; J. M. Pawlowskiand F. Rennecke, Phys. Rev. D (2014) 7,076002 [arXiv:1403.1179 [hep-ph]]; M. Haas, L. Fis-ter and J. M. Pawlowski, Phys. Rev. D (2014)9, 091501 [arXiv:1308.4960 [hep-ph]]; T. K. Herbst,J. M. Pawlowski and B. J. Schaefer, Phys. Rev. D (2013) 1, 014007 [arXiv:1302.1426 [hep-ph]].[37] W. J. Marciano and H. Pagels, Phys. Rept. (1978)137.[38] C. S. Fischer, D. Nickel and J. Wambach, Phys. Rev. D (2007) 094009 [arXiv:0705.4407 [hep-ph]].[39] C. S. Fischer and R. Williams, Phys. Rev. Lett. (2009) 122001 [arXiv:0905.2291 [hep-ph]].[40] H. Sanchis-Alepuz, C. S. Fischer and S. Kubrak, Phys.Lett. B (2014) 151 [arXiv:1401.3183 [hep-ph]].[41] J. S. Ball and T. W. Chiu, Phys. Rev. D , 2542 (1980).[42] G. Eichmann, R. Alkofer, A. Krassnigg and D. Nicmorus,Phys. Rev. Lett. (2010) 201601 [arXiv:0912.2246[hep-ph]].[43] H. Sanchis-Alepuz, G. Eichmann, S. Villalba-Chavezand R. Alkofer, Phys. Rev. D (2011) 096003[arXiv:1109.0199 [hep-ph]].[44] H. Sanchis-Alepuz and C. S. Fischer, Phys. Rev. D (2014) 9, 096001 [arXiv:1408.5577 [hep-ph]].[45] H. Sanchis-Alepuz and R. Williams, arXiv:1504.07776[hep-ph].[46] G. Eichmann, Phys. Rev. D (2011) 014014[arXiv:1104.4505 [hep-ph]].[47] G. Eichmann and C. S. Fischer, Eur. Phys. J. A (2012)9 [arXiv:1111.2614 [hep-ph]].[48] R. Alkofer, G. Eichmann, H. Sanchis-Alepuz andR. Williams, arXiv:1412.8413 [hep-ph].[49] M. Oettel, G. Hellstern, R. Alkofer and H. Reinhardt,Phys. Rev. C (1998) 2459 [nucl-th/9805054].[50] M. Oettel, R. Alkofer and L. von Smekal, Eur. Phys. J.A (2000) 553 [nucl-th/0006082].[51] I. C. Cloet, G. Eichmann, B. El-Bennich, T. Klahnand C. D. Roberts, Few Body Syst. , 1 (2009)[arXiv:0812.0416 [nucl-th]].[52] J. Segovia, I. C. Cloet, C. D. Roberts and S. M. Schmidt,Few Body Syst. (2014) 1185 [arXiv:1408.2919 [nucl-th]].[53] G. Eichmann, A. Krassnigg, M. Schwinzerl andR. Alkofer, Annals Phys. , 2505 (2008)[arXiv:0712.2666 [hep-ph]].[54] G. Eichmann, I. C. Cloet, R. Alkofer, A. Krassniggand C. D. Roberts, Phys. Rev. C (2009) 012202[arXiv:0810.1222 [nucl-th]].[55] D. Nicmorus, G. Eichmann and R. Alkofer, Phys. Rev.D (2010) 114017 [arXiv:1008.3184 [hep-ph]].[56] G. Eichmann and D. Nicmorus, Phys. Rev. D (2012)093004 [arXiv:1112.2232 [hep-ph]].[57] C. S. Fischer, P. Watson and W. Cassing, Phys. Rev. D (2005) 094025 [hep-ph/0509213].[58] W. Kamleh, P. O. Bowman, D. B. Leinweber,A. G. Williams and J. Zhang, Phys. Rev. D (2007)094501 [arXiv:0705.4129 [hep-lat]].[59] A. Maas, J. M. Pawlowski, L. von Smekal and D. Spiel-mann, Phys. Rev. D (2012) 034037 [arXiv:1110.6340[hep-lat]]. [60] K. l. Wang, Y. x. Liu, L. Chang, C. D. Roberts andS. M. Schmidt, Phys. Rev. D (2013) 7, 074038[arXiv:1301.6762 [nucl-th]]. [61] A. Maas, J. Wambach and R. Alkofer, Eur. Phys. J. C42