Heavy quark-antiquark free energy and thermodynamics of string-hadron avoided crossings

# Heavy quark-antiquark free energy and thermodynamics of string-hadron avoided crossings

E. Megías Max-Planck-Institut für Physik (Werner-Heisenberg-Institut), Föhringer Ring 6, D-80805 Munich, Germany    E. Ruiz Arriola    L. L. Salcedo Departamento de Física Atómica, Molecular y Nuclear and Instituto Carlos I de Física Teórica y Computacional
July 21, 2019
###### Abstract

The correlation function between two Polyakov loops encodes the free-energy shift due to a pair of separated colour conjugated sources in the hot QCD medium. This is analyzed in terms of a novel Källén-Lehmann spectral representation for the separating distance, implying an increasing and concave free-energy at all temperatures. We express the heavy free-energy shift below the phase transition in QCD in terms of colour neutral purely hadronic states with no explicit reference to quarks and gluons. Good agreement with lattice data is achieved when considering the avoided crossing mechanism underlying string breaking and with standard quenched values of the string tension known from charmonium and bottomonium phenomenology. We also address the role of the corresponding entropy shift and its renormalization group properties.

finite temperature; heavy quarks; free-energy; Polyakov Loop
###### pacs:
11.10.Wx 11.15.-q 11.10.Jj 12.38.Lg

## I Introduction

In a non-abelian gauge theory such as QCD one can meaningfully address specific information on quark and gluon colour charges and their interactions in the hot medium beyond the EoS which could not be posed in the pre-QCD times. Nonetheless, we expect that any physical observable can still be described in terms of manifestly colour singlet hadronic states below the phase transition. We have recently shown this quark-hadron duality at finite temperature for the quark self-energy via the renormalized Polyakov loop Megias:2012kb () using single heavy meson charm and bottom excited spectra (see also Bazavov:2013yv ()). (See  Megias:2013xaa () for higher representations and Arriola:2014bfa () for a pedagogical discussion.) In this letter we verify that the quark-antiquark free-energy can be represented, as naively expected, via purely hadronic states and a single bosonic string potential, with no explicit reference to quarks and gluons. This requires a non-trivial level avoided crossing structure in the heavy-light meson-antimeson spectrum and is done in harmony with a new Laplace-Stietjels spectral representation of the Polyakov correlator, whence concavity properties of the free-energy at any temperature can be trivially deduced. Some motivation for the present work has been presented previously in Arriola:2015gra (); Megias:2015qya ().

The quark-antiquark free energy is ambiguous by a constant which can be fixed at a reference temperature and relative distance between the heavy sources. On the contrary both the entropy shift and the specific heat shift are unambiguous, and one expects them to fulfill renormalization group invariance. In a short note we have profited from this view to display the entropy shift due to a single heavy source based on integrating the specific heat with suitable boundary conditions Megias:2016bhk (). In the present work we extend such an analysis to the corresponding quark-antiquark situation and show that more information than the specific heat is needed in order to reconstruct the free energy.

The paper is structured as follows. In Section II we present a derivation of the Källen-Lehman representation for the free energy of two colour charged conjugated sources. This allows to easily deduce convexity properties for the free energy. In Section III we analyze the renormalization group aspects of free energy, entropy and specific heat and illustrate the constraints for the case of a medium and the modifications following the addition of a colour source and its charge conjugated source at a given separation distance. In Section IV we analyze the avoided crossing structure of the spectrum made of a string and pairs of heavy-light, and singly heavy baryons and the corresponding partition function containing the Hadron Resonance Gas and a string. We highlight the important role played by the string-hadron transition on the light of a comparison with available recent lattice data for the Polyakov loop correlator. Finally in Section V we summarize our points and come to the conclusions.

## Ii Källén-Lehmann spectral representation for the correlation of two Polyakov loops

When two heavy colour charge conjugated sources belonging to the representations and of the colour gauge group are created and placed at a given distance in the hot medium, there arises a free-energy which provides the maximum work the system can exchange with the medium at a fixed temperature. McLerran and Svetitsky suggested 35 years ago to explore these free-energies for the fundamental representation as suitable order parameters for a hadronic-quark gluon plasma phase transition McLerran:1980pk (); McLerran:1981pb () (for an early review see e.g. Svetitsky:1985ye ()).

### ii.1 Standard scheme

The obvious approach to the heavy quark-antiquark free-energy is to consider a compactified Euclidean time to include the finite temperature, so the three coordinates are spatial ones and is time-like (see Fig. 1, left). Denoting the partition function of the system without sources, the partition function with heavy sources a distance apart, and the correlation function between Polyakov loops,

 C(r,T)=⟨trRΩ(r)tr¯RΩ(0)⟩T=ZR⊗¯R(r,T)Z0(T), (1)

where the Polyakov loop, , is a purely gluonic operator, , and is the thermal expectation value. Correspondingly, for the shift in the free-energy

 F(r,T)≡FR⊗¯R(r,T)−F0(T)=−TlogC(r,T). (2)

Cluster decomposition and translational invariance requires factorization (using also colour charge conjugation )

 C(∞,T)=⟨trRΩ⟩2T=e−2F(T)/T (3)

where is the shift in the free-energy produced by introducing a single heavy charge.

To avoid paradoxes it is important to note that, unlike and , the shift is not a true free-energy. This quantity is not extensive and thermodynamic stability does not require the entropy shift

 S(r,T)=−∂F(r,T)∂T, (4)

to be an increasing function of the temperature. Likewise, is not a true partition function (i.e., a sum of states with non negative integer degeneracies) but rather the ratio of two partition functions Luscher:2002qv (). Nevertheless, by the same arguments given in  Megias:2012kb (); Megias:2013xaa () for the Polyakov loop, this ratio should behave, to a good approximation, like a partition function for temperatures below the crossover, that is,

 C(r,T)≈∑ngne−En(r)/T=∫∞−∞dEρ(E,r)e−E/T, (5)

where is non negative. This implies that standard thermodynamic relations apply to in this low temperature regime. Similar remarks apply for and  Megias:2012kb (); Megias:2013xaa ().

### ii.2 Alternative scheme

Additional information can be obtained from the alternative point of view in which are regarded as spatial coordinates and as Euclidean time (see Fig. 1, right). In this picture (the inverse physical temperature) is the size of the system in the spatial direction. The size remains infinite in the and directions. The system is described by a Hamiltonian which produces the evolution in the direction and depends on . Since the size in the direction is unbounded the effective temperature is zero. Therefore, in this picture, the expectation values are taken in the vacuum state of , which we denote .

To carry out a canonical treatment, a convenient choice is the static gauge (where is the gluon gauge field in the direction). This leaves the Polyakov loop in the direction, , together with and , , as degrees of freedom. acts in the Hilbert space as a multiplicative operator. In this view the expectation value of two Polyakov loops, one at and the other at , and a later time , represents the action of the operator (at ) on the vacuum at time followed by evolution to time to connect again with

 C(r,T)=⟨0,T|TrRΩe−rHz(T)TrRΩ†|0,T⟩. (6)

This allows to establish an exact spectral representation valid for any representation and any temperature, by inserting a complete set of eigenstates of , ,

 C(r,T)=∑n|⟨n,T|TrRΩ†|0,T⟩|2e−rwn(T):=∫∞0dwτ1(w,T)e−rw. (7)

The density of states excited by the Polyakov loop, , is a non-negative function with support on ,111Note that the vacuum has zero energy, since replacing by the identity operator would yield for all . so Eq. (7) is a Laplace-Stietjels transformation book:1304498 ().

The sum rule naturally splits into its -independent disconnected component, coming from the vacuum term (i.e., the term in the sum), and the connected component from the excited states:

 C(r,T)=L(T)2+Cc(r,T),L(T)≡⟨trRΩ⟩T. (8)

The disconnected component yields a contribution in .

In addition to the usual scheme and the one considered here, , a third scheme would be being now the time direction. In this case is just the expectation value of an operator, , acting at time on the vacuum (the same vacuum as in the scheme). No excited states are involved in this case.

### ii.3 Convexity properties of the free energy

An interesting consequence of the spectral relation (7) comes from the fact that is positive and has positive support. This automatically implies that is an absolutely monotonic function with respect to , that is,

 (−1)n∂nC(r,T)∂rn≥0,n=0,1,2,… (9)

Further bounds on the free-energy can be derived by rewriting these inequalities as expectation values with as measure:

 (−1)nC∂nC∂rn=1C∫∞0dwτ1e−rwwn≡⟨wn⟩. (10)

In this view because both and the measure are non negative. For positive , the measure is exponentially convergent at large , since that regime of energies is dominated by the Coulomb interaction between the two conjugated charges (see below).

Letting , such inequalities imply,

 f′=⟨w⟩≥0,f′2−f′′=⟨w2⟩≥0,… (11)

The prime denotes derivative with respect to . Tighter bounds can be obtained by using optimized positive polynomials.222The coefficients of the polynomials have been chosen so as to remove higher powers of , i.e., monomials with many ’s. The coefficients in front of and cannot be reduced since in that case the expressions would not be positive definite for and arbitrary positive . Specifically

 ⟨w⟩=f′≥0,⟨(w−⟨w⟩)2⟩=−f′′≥0⟨w(w−⟨w⟩)2⟩=f′′′−f′f′′≥0⟨(w2−2⟨w⟩w+2⟨w⟩2−⟨w2⟩)2⟩=2f′′2−f(iv)≥0. (12)

These inequalities hold for all . The first two relations imply

 ∂F(r,T)∂r≥0,∂2F(r,T)∂r2≤0. (13)

An immediate consequence is that the zero temperature heavy quark-antiquark potential must be an increasing and concave function of ,

 V′(r)≥0,V′′(r)≤0, (14)

a result previously established from reflection positivity bachas1986concavity (); Nussinov:2000kn (). A similar representation to ours on the lattice was proposed long ago Kiskis:1985vf () although the important new implications discussed here were not addressed. A further property from the Laplace-Stietjels representation is that the analytical continuation to the complex plane has no singularities for .

### ii.4 Källén-Lehmann spectral representation

The spectral relation in Eq. (7) can be improved by taking into account full translational and rotational invariances in the three non-compactified directions. Indeed, up to now we have only considered expectation values of in a single state namely, (here is the origin of the plane), however, one can take more general states of the type and still must be positive for any choice of the coefficients . Equivalently, the function

 h(r1,r2;z)≡⟨0|trRΩ(r1)e−zHztrRΩ†(r2)|0⟩,z≥0, (15)

provides the matrix elements of a certain positive operator in ,

 h(r1,r2;z)=⟨r1|O(z)|r1⟩, (16)

and furthermore

 h=h(|r1−r2|2+z2). (17)

In view of the translational invariance, the simplest way to impose the positivity condition on is to work in momentum space where this operator is diagonal. The details of the derivation are given in Appendix A. One finds a tighter exact spectral representation for the connected component of

 Cc(r,T)=∫∞0dμτ(μ,T)e−μr4πr, (18)

for some non-negative density . This is nothing else than the usual Källén-Lehmann spectral representation but in three dimensions: behaves like an ordinary scalar field at zero temperature in a three dimensional space-time with Lorentz invariance. As it turns out, a formula equivalent to the one derived here was already noted in Luscher:2004ib () in the context of Yang-Mills theory with unbroken center symmetry.

The new spectral relation contains the previous one, Eq. (7), through the relation

 τ1(w,T)=L(T)2δ(w)+14π∫w0dμτ(μ,T). (19)

As a consequence we learn that the connected component of is not only positive but also an increasing function of . Also it is important to note that is an energy while is an invariant mass, thus the support of may have a discrete part below the two particle threshold.

When the Polyakov loop is not present, (e.g., Yang-Mills in the confined phase) it is easy to derive tighter conditions on the free energy shift from the new representation in Eq. (18), since then has all the properties previously exploited for . So in this case, again letting , it follows that , , , etc. In particular,

 ∂F∂r≥Tr,∂2F∂r2≤−Tr2(L(T)=0). (20)

The first relation implies that, at finite temperature, must increase at least logarithmically (namely, as ) with the separation, as a direct consequence of the factor in the sum rule, whenever the Polyakov loop vanishes.

The determination of the invariant-mass spectral density seems of interest since this quantity contains much condensed information on the system composed of two heavy conjugated charges at finite temperature. As it turns out, such determination can be carried out for in some simple cases. A first case is

 F(r,T)=−αr (21)

which should describe quenched QED (a free theory) and entails at all temperatures. Let

 Gn(x)≡In(2√x)xn/2=0F1(n+1;x)n!=∞∑k=0xkk!(n+k)!. (22)

Here is the modified Bessel function and is the confluent hypergeometric functions. The functions are positive for positive and grow as for large . In addition, . We rely on the identity

 ∫∞0dxG1(x)e−rx=e1/r−1(r>0). (23)

The correlation function follows from

 τ1(w,T)=δ(w)+αTG1(αTw)θ(w). (24)

Here is the Heaviside step function. In turn,

 τ(μ,T)=4πdτ1,c(μ,T)dμ=4πα2T2G2(αTμ)θ(μ)+4παTδ(μ). (25)

In Eq. (24), gives rise to the Polyakov loop, while in Eq. (25), accounts for the additional falloff in the connected part of for large separations.

Another case where the spectral function takes a closed form is that of a simple model of the type string tension plus Coulomb. Specifically, for the internal energy shift

 U(r,T)=σRr−αRr+U0, (26)

where is a constant. This form of is suited to model the low temperature regime of a Yang-Mills theory. The thermodynamic relations require then the entropy shift to be a function of only. The choice

 S(r,T)=−log(4πrμ0) (27)

allows to fulfill the spectral representation, being an arbitrary scale. The corresponding correlation function takes the form

 C(r,T)=14πrμ0e−(σRr−αR/r+U0)/T. (28)

In this model . Such correlation then follows from

 τ(μ,T)=e−U0/Tμ0[δ(μ−σRT)    +αRTG1(αRT(μ−σRT))θ(μ−σRT)]. (29)

A Polyakov loop can be added in the correlation function without touching but in this case and are no longer given by the previous expressions.

In Eq. (29), the Dirac delta term represents the mass of the lowest-lying state, a string of length , which dominates the long distance tail of the correlation in the string tension plus Coulomb model. On the other hand, the Coulomb-induced term increases as an exponential of and it should saturate the large region of of full QCD, albeit softened by asymptotic freedom. This behaviour ensures the convergence of the measure at large for .

In Yang-Mills theory one should expect a gap in the mass spectrum below the transition temperature, giving rise to a string tension. Indeed, a gap in the spectrum above the vacuum state makes the integration in the spectral sum rule to start at the mass of the lightest state and this term dominates the large behaviour of the correlation function. Assuming that is an isolated non vanishing point in the spectrum (the lightest state is expected to be a single particle state) one finds, for large ,

 C(r,T)≈r→∞C0(T)e−μ1(T)r4πr (30)

and the quantity can be interpreted as a (possibly -dependent) string tension. An increase in should translate into a continuous quenching of the gap which becomes zero at the transition temperature. Above this temperature couples to the vacuum state and a non vanishing Polyakov loop expectation value emerges.

If there were a gap also in full QCD one would have instead

 C(r,T)≈r→∞C0(T)e−μ1(T)r4πr+L(T)2 (31)

but still could be interpreted as a (once again -dependent) string tension. The dominance of the Polyakov loop at large would display the string breaking.

However, a gapless spectrum seems more likely for full QCD. A vanishing gap is more suited to describe the deconfined phase due to the lack of a string tension there. Then by continuity the gap must be zero for all temperatures, since in QCD there is a crossover rather of a true phase transition Aoki:2006we (). Also the fact that is non zero for all temperatures and representations suggests a vanishing gap. Nevertheless, even if strictly speaking, the support of the spectral function fills the positive half-line, it is not excluded that, in the confined phase, this function have a narrow peak around some . In this case one can speculate that a precise definition of a QCD string-tension could still be obtained as a pole in the complex plane.

## Iii Renormalization group

### iii.1 Renormalization of the free energy

In the extraction of physical information from lattice calculations the renormalization of the free-energy is a crucial issue. At distances shorter than the thermal wavelength, , one expects that the medium plays a minor role and, due to asymptotic freedom, perturbative QCD applies  Brambilla:2010xn (); Burnier:2009bk (). This requirement has been emphasized in a series of insightful works Kaczmarek:2002mc (); Zantow:2003uh (). A limitation is that the necessary small values of might not be attainable in current lattice settings. More recently, a full calculation with realistic dynamical quark masses has been carried out for fundamental sources Borsanyi:2015yka (); Bazavov:2016qod (). There, the lattice action is renormalized at zero temperature, i.e., by following lines of constant physics. The only new ultraviolet ambiguity introduced by the Polyakov loop operators is a single - and - independent (although -dependent) additive constant in , i.e. , which in principle can be fixed by setting for conventionally chosen , , and .

We will review below the renormalization group equation (RGE) for the vacuum case Ellis:1998kj (); Shushpanov:1998ce () and the heavy quark situation Chernodub:2010sq (). Our reanalysis is based on the entropy and introduces some minor but key modifications. The upshot is that while one can determine the entropy shift from the RGE in the vacuum and single heavy quark case, one cannot do the same in the case of heavy quark-antiquark sources separated a given distance.

### iii.2 Renormalization group for the entropy

#### iii.2.1 Medium

While the entropy or the free energy are of high theoretical interest, experimentally one can only measure the specific heat. Thus we naturally expect that this quantity be expressed in a manifestly renormalization group invariant way. To fully appreciate this point let us consider the simplest case, of a pure gauge theory which is specified by a renormalization scale and a dimensionless coupling constant , at temperature T and volume . The partition function and its relation to the free energy is defined as

 Z=e−F/T=∫DAe−∫d4xL (32)

where the Lagrangian is given by

 L=14(Gaμν)2. (33)

From the standard thermodynamic relation we have the entropy

 S=−∂F∂T=∂∂T(TlogZ). (34)

As any dimensionless quantity, the entropy must fulfill a functional dependence involving just dimensionless quantities. If we take , , and as the relevant variables we have

 S(g,μ,T,V)=ϕ(g,log(μ/T),log(μV13)). (35)

The renormalization group invariance means that

 μddμS(g,μ,T,V)=0 (36)

whence we obtain

 T∂TS−3V∂VS=β(g)∂gS, (37)

where the beta function is defined as

 β(g)=μ∂g∂μ. (38)

To evaluate the derivative with respect to we rescale the gluon field so that the measure and the action scale as

 DA=D¯A/gN,L=14g2(¯Gaμν)2, (39)

with the number of lattice points with the lattice spacing, and is independent of . Thus,

 ∂gS=∂T(T∂glogZ)=∂T(TVTa4+2g3V⟨(¯Gaμν)2⟩T)=V2g∂T⟨(Gaμν)2⟩T.

Defining the trace of the energy momentum tensor as

 Θ=β(g)2g(Gaμν)2 (40)

we get finally

 ∂T(E−3PV)=∂T[T∫d4x⟨Θ⟩T]. (41)

In the case of infinite volume we just get the energy density and and the volume factors out

 ∂T(ϵ−3P)=∂T⟨Θ⟩T. (42)

Therefore, integrating from to one gets

 ϵ−3P=⟨Θ⟩T−⟨Θ⟩0, (43)

where we have assumed that for . This follows from the low temperature partition function behaviour where the mass of the lightest glueball. The derivation of Ref. Ellis:1998kj (); Shushpanov:1998ce () starts already with the partition function in terms of the scaled fields and thus the final equation does not include the subtracted contribution at zero temperature, which has to be added by hand (see also Chernodub:2010sq ()).

In the full QCD case we have to add quark and anti-quark fields

where the QCD Lagrangian for flavours reads, in terms of the re-scaled gluon field with and ,

 L(x)=14g2(¯Gaμν)2+∑q=u,d,s¯q(i/D+mq)q. (45)

Renormalization group invariance requires the inclusion of the mass terms by assuming an extra dependence on the dimensionless variable , yielding

 μdSdμ = β(g)∂S∂g−∑qmq(1+γq)∂S∂mq−T∂S∂T+3V∂S∂V (46) = 0,

with the beta function and the mass anomalous dimension given by

 β(g)=μdgdμ,γq(g)=−dlogmqdlogμ. (47)

We obtain the same formulas as above with the energy momentum tensor defined as

 Θ≡Θμμ=β(g)2g(Gaμν)2+∑qmq(1+γq)¯qq. (48)

Here, in the low temperature limit and are saturated by the free pion gas which again provides a vanishing contribution of the type .

#### iii.2.2 Heavy Source

For an operator the thermal expectation value is defined as

 ⟨O⟩T=∫DAOe−∫d4xL(x)∫DAe−∫d4xL(x). (49)

A qualification is in order here. Generally, not every thermal expectation value corresponds to a ratio of partition functions. Strictly speaking a partition function requires a spectral decomposition of the form where is a non-negative integer. We will only consider operators which indeed correspond to partition functions (in particular must be a dimensionless projector operator)

 ZOZ=e−ΔFO/T (50)

where we have introduced the free energy shift,

 ΔFO=FO−F. (51)

To this free energy shift there corresponds an entropy shift,

 ΔSO=−∂ΔFO∂T. (52)

The renormalization group equation becomes

 T∂TΔSO−3V∂VΔSO=∂T[T∫Vd4x⟨Θ⟩O,T] (53)

where

 ⟨Θ⟩O,T=⟨ΘO⟩T⟨O⟩T−⟨Θ⟩T. (54)

We can take the operator as the Polyakov loop in any representation and, in the infinite volume case, the volume dependence term drops out. Furthermore, for this choice of the thermal expectation value corresponds to the ratio of two true partition functions (Eq. (1)). We will work in the static gauge, in which the Polyakov loop reads () and we take conventionally . In this case the RGE reads

 T∂TΔSR=∂T[T∫d4x⟨Θ⟩R,T]. (55)

This formula provides the specific heat for placing a colour charge in the representation into the hot medium. Remarkably this equation allows to uniquely determine the entropy after specifying its value at a given temperature.

In the limit of small temperatures we have

 ΔSR(0)=logDR, (56)

where is the ground state degeneracy in the subspace with a colour source in the representation at . In particular, if the source is in the fundamental representation the lowest state is obtained by screening the source by a single light anti-quark, hence for mass-degenerated light quark flavors (the factor coming from the two spin states of the anti-quark)  Megias:2012kb (). The different colour states of the source are combined with those of the light quark to form a colour singlet, so this degree of freedom does not add to the entropy.

In the opposite limit of large temperatures a partition function is just the dimension of the Hilbert space of the system. Hence where is the dimension of the representation , thus

 ΔSR(∞)=logdR. (57)

This expresses the fact that at very high temperatures the colour state of the source is not effectively correlated with the medium, so it just counts additively for the entropy.333Recall that a colour source has no other degrees of freedom than colour, by definition. A heavy quark serves as a source, not only because it does not exchange kinetic energy but also because its spin state is fully decoupled and can be disregarded. [For the free energy defined from the Polyakov loop normalized as , this result takes the form .]

#### iii.2.3 Heavy source correlator

In the case of a Polyakov loop correlator there appears the separation between conjugate sources, , which enters in the RGE by including the dimensionless quantity and yielding the replacement , thus

 T∂TΔS¯R⊗R−r∂rΔS¯R⊗R=∂T[T∫d4x⟨Θ⟩¯R⊗R,T]. (58)

It is noteworthy that the previous equation is free from UV divergences, although the integrand itself can display UV divergences at the heavy source points. This is at variance with a RGE for the free energy Chernodub:2010sq ().

Eq. (58) is a standard first order partial linear differential equation which can be solved by the method of characteristics. Assuming that the r.h.s. is known, the equation provides the variation of along paths in the plane. Specifically

 T∂TS−r∂rS=ϕ(r,T) (59)

is equivalent to

 S(r,T)=S(rT/T0,T0)+∫TT0ϕ(rT/T′,T′)dT′T′. (60)

The determination of from the RGE can then be achieved from the knowledge of the entropy along a line visiting all the values.

The explicit solution Eq. (60) applies immediately if the entropy is known for all at some reference temperature . For instance in the case, at low temperatures we expect

 limT→0ΔS¯QQ(r,T)={0rrc (61)

where is kept constant as , and is the string breaking distance. This region of small is beyond a perturbative calculation. On the other hand, at high temperatures we expect a free theory, with

 limT→∞ΔS¯QQ(r,T)=2log(Nc). (62)

In this limit, lies in the perturbative region since is also small. Nevertheless, likely, any perturbative corrections will be exponentially inflated by the RGE, rendering the determination of for finite unreliable.

### iii.3 Perturbation theory and RGE improvement

#### iii.3.1 Heavy source in the medium

One useful application of the RGE is the derivation of constraints on perturbative results. For instance, the expectation value of the Polyakov loop has been computed to in pQCD in Gava:1981qd (); Brambilla:2010xn (); Burnier:2009bk () and to in Berwein:2015ayt (). To this order, the structure of this quantity is

 1NcL(T)=1+c0g3+(c1+d1logg)g4+(c2+e2log(μ/2πT))g5+O(g6), (63)

hence, for the entropy at NNLO

 ΔSQ=c0g3+(c1+d1logg)g4+(c2−e2+e2log(μ/2πT))g5+O(g6). (64)

The RGE then requires , where

 β(g)=−g3∑n≥0βngn. (65)

A determination of the specific heat

 ΔcQ≡T∂TΔSQ (66)

by direct derivation with respect to with fixed would produce just the LO result (of ). However, such a calculation would disregard the fact that some of the -dependence at higher orders is fixed by the known lower orders through the RGE. The RGE guarantees that to any given order we can differentiate (with fixed ) and then take (e.g. ) or the other way around. The latter method recovers the NNLO result for . Alternatively, one can work with fixed and use the RGE relation

 T∂TO=β(g)∂gO, (67)

which holds for any observable in which is the only physical scale. This gives

 (68)

Note that the function to NNLO is needed.

#### iii.3.2 Heavy source correlator

The perturbative calculation of the Polyakov loop correlator was done by Nadkarni Nadkarni:1986cz () in the regime ( is the Debye mass) yielding a function of for the connected piece. Such dependence vanishes under the RGE for the entropy, up to higher orders in perturbation theory stemming from the energy-momentum tensor contributions. This illustrates explicitly that the RGE of the entropy shift is not sufficient for the full reconstruction of the correlator. In the calculation of Ref. Brambilla:2010xn (), in the regime , a non vanishing contribution to the RGE is explicitly obtained, involving again higher orders. Therefore, in order to test these perturbative results, a different object should be computed, namely, the energy momentum density of two charge conjugated heavy sources at finite temperature. To our knowledge such a calculation has never been carried out in any form in the literature.

## Iv Heavy ¯QQ free energy shift in the confined phase

On physical grounds one expects that below the phase transition the free-energy shifts should be expressed in terms of hadronic colour singlet states, but the precise mathematical formulation of this expectation has never been made clear. We will address in this section the hadronic representation of the heavy free-energy shifts, based on our findings of Sec. II. A hadronic representation for the Polyakov loop was already derived in Megias:2012kb () based on first principle arguments of QCD. The results of this section generalize that study.

### iv.1 ¯QQ spectrum and string breaking

At zero temperature and for fundamental sources, the concept of string tension has played a major role in the formulation and understanding of confinement and colour gauge invariance in the pure Yang-Mills theory. In the zero temperature limit one has for large distances . In QCD however, the string between fundamental charges breaks generating a pair from the vacuum which subsequently decays into hadronic states, and so instead where is the mass of the lowest-lying bound state (a hybrid heavy-light meson). The breaking of the string has been studied on the lattice Bali:2005fu () where the avoided crossing between the and the channels was observed, which is familiar from molecular physics in the Born-Oppenheimer approximation landau1965quantum (). More generally, the state can decay into any of the many excited states of the meson spectrum, or the baryon spectrum, as long as they have the same quantum numbers as the system.444The mechanism of decaying into baryons implies the creation of two pairs of light quarks , leading to the formation of two heavy-light baryons with one heavy quark.

To analyze the free-energy lattice results in the confined phase we will apply the approximate -spectral representation in Eq. (5) and the exact -spectral representation in Eq. (7). Further, we will consider the coupled-channel Hilbert space spanned by (a single state representing the colour sources joined by a string in its ground state) and by pairs of heavy hadrons. Here is a colour source screened with light quarks and gluons to form a singly-heavy hadron (meson or baryon), either low-lying or excited. As in Megias:2012kb (), half of the heavy hadron spin states are spurious, since the colour source has no spin. When this source is simulated with a very heavy quark, its spin decouples due to QCD heavy quark spin symmetry Isgur:1989vq (); Neubert:1993mb (). The energies are modeled as

 V¯QQ(r)=−43αr+σr+c,V(n,m)¯HH(r)=Δ(n)¯H+Δ(m)H, (69)

where (-th heavy-light hadron mass minus heavy quark mass). An additive constant has been included in to account for the ambiguity in the renormalization of the lattice data (see Sec. III.1) In addition, a transition potential can be present.

In principle one should take the infinite quark mass limit () and compute the corresponding spectrum. This is an ambitious program and as a guide we will content ourselves with using existing extensive quark model calculations containing the main essential features describing singly heavy hadrons (mesons and baryons) and relativity, which proves crucial to account for excited states. Obviously this represents an approximation but these models already produce a trace anomaly for , , quarks, below the phase transition, which can hardly be distinguished by the conventional hadron resonance gas using the listed PDG values Arriola:2014bfa (). Besides, one can assess the corresponding uncertainty by comparing the -hadron versus the -hadron spectra and this is in accordance with previous analyses for the Polyakov loop, cf. Megias:2012kb (); Megias:2012hk (). Moreover, we will disregard gluonic string-like excitations as we expect them to have a gap of the order of .

In Fig. 2 we illustrate the situation for and , and using the spectrum of heavy-light hadrons with a charm quark from the Godfrey-Isgur relativized quark model (RQM) Godfrey:1985xj () up to . As mentioned, a source of error is that the heavy quarks in nature have a finite mass. In order to estimate this, we have used as well the spectrum of heavy-light hadrons with a quark from the very same model. Besides reproducing accurately de PDG level density, the RQM satisfactorily describes the Polyakov loop below the crossover and, for the heavy-light system with charm or bottom quarks, yields an exponentially growing spectrum with Polyakov-Hagedorn temperature Arriola:2013jxa (). Before mixing, all levels cross with the quenched potential (note the relatively large gaps between lower levels mixing with the channel). The avoided crossings correspond to states with the same quantum numbers as the , which are also highlighted in Fig. 2 for the RQM using a mixing to be discussed subsequently.

### iv.2 Hadronic representation for the free energy

Relying on the -spectral representation and adding up all coupled-channel states, in the absence of mixing one finds Arriola:2014bfa ()

 e−F(r,T)/T=e−V¯QQ(r)/T+(∑ne−Δ(n)H/T)2. (70)

The last term is just  Megias:2012kb (). The simple form of this term follows from neglecting any interaction between the two heavy hadron states (Eq. (69)).

Using the data of Ref. Borsanyi:2015yka () to extract and , one can determine by inverting Eq. (70). This yields a -independent potential for but with a string tension , almost twice the conventionally accepted value . This is shown in Fig. 3 (up). This disturbing result suggests that determining the string tension or renormalizing the free-energy-shift from short distances (at least in lattice QCD with realistic quark masses and current lattice spacings) may contradict the long distance and well established charmonium phenomenology, based on quenched determinations of the string tension.

In Fig. 3 (down) we display the lattice results of Borsanyi:2015yka () for the heavy free energy along with the calculation from Eq. (70) using the Isgur model with charm quarks and . A better agreement with lattice data is achieved mainly for the lowest temperatures. The disagreement at higher temperatures is in part motivated by the failure of the hadronic representation of the Polyakov loop to describe the lattice data at higher temperatures, cf. Megias:2012kb ().

In order to estimate finite mass effects, the above determination of from Eq. (70) can be repeated using instead quarks in the Isgur spectrum. One finds a maximum relative change in this quantity of for and of for . Such maximum deviations take place at large separations as they reflect the slightly different Polyakov loop predicted by the two versions, already noted in Ref. Megias:2012kb (). The variation is larger at the highest temperature, where the model is less reliable. As we explain subsequently, in the extraction of precise quantities like the mixing potential in Sec. IV.3, we take the Polyakov loop directly from lattice data rather than from a hadronic model. This fixes the large separation regime, hence finite mass effects are restricted to intermediate distances, and they are significantly smaller than the above estimates.

Remarkably, Eq. (70) is fully consistent with the exact -representation, Eq. (7). However, when channel mixing is allowed and the modified eigenvalues are introduced in the -spectral representation, the admissible mixings get constrained by the concavity relations in Eq. (13). For instance, in the low temperature limit the lowest eigenvalue dominates and this implies the conditions and (no such direct constraint applies to the excited states). The two-level setting in Eq. (72) (below) illustrates that cannot be arbitrary: for large enough (the sign is conventional) and so the concavity relations require and .

### iv.3 Avoided crossings

When mixing is switched on, we still use an additive model for , i.e., interactions between the two heavy hadrons are neglected, but non vanishing matrix elements arise between different pairs, as well as transitions between the states and the state. The simple form in Eq. (70) no longer holds since the Hamiltonian is no longer diagonal in the coupled-channel space. If denotes the Hamiltonian (a matrix in coupled-channel space) the formula with mixing becomes

 e−F(r,T)/T=tr(e−V(r)/T)=∑αe−Eα(r)/T, (71)

where denote the corresponding eigen-energies.

To study the phenomenon of avoided crossings in closer detail and paralleling the treatment in Bali:2005fu (), let us introduce a mixing just between the state and the lowest-lying state

 V(r)=⎛⎜ ⎜ ⎜⎝−4α3r+σrW(r)W(r)2Δ⋱⎞⎟ ⎟ ⎟⎠. (72)

We only display the two coupled levels of , this matrix being diagonal for the remaining states. To be more precise, the lowest lying level contains several spin-flavor states. To match the quantum numbers of the state, we only couple the flavorless and spinless combinations. The uncoupled states give a contribution as in Eq. (70) and their role is to saturate the Polyakov loop. Therefore their detailed eigen-energies are not needed. The contribution obtained from diagonalization of the two first states in in Eq. (72) is computed explicitly and the remainder is fixed so that the Polyakov loop value is reproduced.

This simplest mixing model allows to determine point by point from . A crucial consistency check of the approach is that must be -independent. Clearly, our string-hadronic model will not be valid above the critical temperature except for very small separations. Therefore, we expect that extraction of from the data will reflect this short distance medium independence. As we can see from Fig. 4 (up) the short distance behaviour is pretty much independent of the temperature provided we take the standard value for the string tension, . Moreover, in the common range for temperatures an exponential behaviour is obtained (at least above ), in agreement with the functional form proposed in Refs.  Drummond:1998ir (); Drummond:1998he (); Drummond:1998ar (). A fit with to the lattice data at with and (the lowest Isgur state with a charm quark) produces

 m=0.80(38)GeV, (73)

with . Parameters and are highly positively correlated, with a correlation .

For , which corresponds to the lowest Isgur state with a bottom quark, one obtains , , and , with and .

We have checked that including more states in the mixing just reduces the strength of the transition potential but does not change the exponential behaviour, cf. Fig. 4 (down).

Note that the highest temperature used in this analysis, is clearly above the physical value of in QCD which is around . The value of in Ref. Borsanyi:2015yka () is obtained in Aoki:2006br (); Aoki:2009sc () from chiral susceptibilities, namely,