Dependence of the static quark free energy on $\mu_B$ and the crossover temperature of $N_f = 2+1$ QCD

We study the dependence of the static quark free energy on the baryon chemical potential for $N_f = 2+1$ QCD with physical quark masses, in a range of temperature spanning from 120~MeV up to 1~GeV and adopting a stout staggered discretization with two different values of the Euclidean temporal extension, $N_t = 6$ and $N_t = 8$. In order to deal with the sign problem, we exploit both Taylor expansion and analytic continuation, obtaining consistent results. We show that the dependence of the free energy on $\mu_B$ is sensitive to the location of the chiral crossover, in particular the $\mu_B$-susceptibility, i.e. the linear term in $\mu_B^2$ in the Taylor expansion of the free energy, has a peak around 150 MeV. We also discuss the behavior expected in the high temperature regime based on perturbation theory, and obtain a good quantitative agreement with numerical results.


I. INTRODUCTION
Heavy quark free energies have been used as a probe of the confining properties of strong interactions since the early days of lattice QCD simulations. They can be extracted, after proper renormalization [1][2][3][4][5], from the expectation value of the Polyakov loop and of its correlators. The Polyakov loop is defined in the continuum as where T is the temperature, P is the path-ordering operator, and N c is the number of colors. On the lattice, this object is constructed by taking the product of gauge links winding along the compactified Euclidean temporal direction. The square module of its trace is the asymptotic value of the unsubtracted correlator between Polyakov loops: it is related to the static quark free energy F Q by the formula In the pure gauge theory the Polyakov loop is an exact order parameter for color confinement/deconfinement, which becomes non-zero only in the deconfined phase and signals the spontaneous breaking of center symmetry. This is usually associated with the possibility of separating two static color charges at arbitrarily large distances without paying an infinite amount of free energy.
In full QCD the situation is different: the creation of dynamical quark-antiquark pairs makes the free energy * Electronic address: massimo.delia@unipi.it † Electronic address: francesco.negro@davigonicoloso.edu.it ‡ Electronic address: andrea.rucci@pi.infn.it § Electronic address: francesco.sanfilippo@roma3.infn.it of static quark pairs finite at any distance even in the confined phase. In fact, dynamical quarks break center symmetry explicitly, so that the Polyakov loop is not an exact order parameter any more and its expectation value is different from zero even in the confined phase.
In the presence of physical quark masses chiral symmetry is surely a relevant symmetry, even if not exact, and the chiral condensate and its susceptibility are usually adopted as probes to locate the pseudo-critical temperature of QCD, which is found to be around 155 MeV [6][7][8][9][10][11]. Still, the Polyakov loop shows a rapid rise at a similar temperature scale, signalling the passage to a deconfined regime with screened color interactions.
Whether deconfinement and chiral symmetry restoration take place at exactly the same temperature is yet not clear and maybe not even a well founded question. The Polyakov loop susceptibility shows a peak around 200 MeV [12], while other related observables show a signal closer to the chiral transition temperature: this is the case for the Polyakov loop entropy S Q = −∂F Q /∂T [12] or the so-called transverse susceptibility related to fluctuations in the imaginary part of the Polyakov loop [13,14]. Since the QCD transition is actually a crossover, it is quite natural to expect that different observables yield different locations of the pseudo-critical temperature. Yet, the information coming from different probes can be useful to better understand the connection between different phenomena taking place around the crossover region.
The purpose of the present study is to give a closer look at static quark free energies, in particular by exploring their dependence on the baryon chemical potential µ B . The modification of the heavy quark free energy due to µ B , ∆F Q (T, µ B ) ≡ F Q (T, µ B ) − F Q (T, 0), is given by the following expression which does not need renormalization if the two Polyakov loops in the ratio are computed at the same ultraviolet (UV) scale. This quantity has been studied in Ref. [15] and more recently in Ref. [16] for QCD with physical quark masses. One expects the dependence of ∆F Q (T, µ B ) on µ B to be sensitive to the location of the transition. Indeed, if the Polyakov loop were an exact order parameter then its dependence on µ B should become singular at T c , because µ B is a relevant parameter which modifies the location of T c . A remnant of this behavior must be present even when the Polyakov loop is not an exact order parameter and, since the free energy is an even function of µ B , the first non-trivial derivative to investigate the associate pseudocritical behavior is the mixed susceptibility Early simulations of N f = 2 QCD have shown that this quantity has a broad peak in a region close to T c [15]. More recent simulations, performed for N f = 2 + 1 QCD discretized via stout-staggered fermions with physical quark masses [16], were limited to a temperature range T 180 MeV, showing nevertheless a peculiar behavior pointing to a seeming divergence for T ∼ 150 MeV.
The purpose of the present study is to extend the investigation for N f = 2 + 1 QCD with physical quark masses to a wider temperature range, going from 120 MeV up to 1 GeV. We consider the same stout-staggered discretization adopted in Ref. [16] and two different sets of lattice spacings, corresponding to Euclidean temporal extensions N t = 6 and N t = 8, in order to estimate the impact of systematic errors related to the UV cutoff. The extended range of temperatures will permit us both to investigate the pseudocritical behavior of χ Q,µ 2 B around T c , and to compare results obtained at high T with perturbative predictions. Since lattice simulations at non-zero µ B are not feasible, because of the sign problem, we employ both Taylor expansion and analytic continuation from simulations at imaginary µ B in order to properly cover the whole temperature range: for temperatures where both methods are used we obtain consistent results.
The paper is organized as follows: in Section II we review our numerical methods and the observables explored in this study; results are presented in Section III and, finally, in Section V, we draw our conclusions.

II. NUMERICAL SETUP AND OBSERVABLES
We have considered the finite temperature partition function for N f = 2 + 1 QCD with chemical potentials µ f (f = u, d, s) coupled to quark number operators, Z(T, µ u , µ d , µ s ), in a setup for which µ u = µ d = µ s = µ B /3, corresponding to a purely baryonic chemical potential. The path integral formulation of Z(T, µ B ), discretized via improved rooted staggered fermions and adopting the standard exponentiated implementation of the chemical potentials [17,18], reads where is the tree-level Symanzik improved action [19,20] (W n×m i;µν stands for the trace of the n × m rectangular parallel transport in the µ-ν plane and starting from site i), and the staggered fermion matrix is defined as where U (2) i;ν are two-times stout-smeared links, with isotropic smearing parameter ρ = 0.15 [21]. Bare parameters have been set so as to stay on a line of constant physics [22][23][24], with equal light quark masses, m u = m d = m l , a physical strange-to-light mass ratio, m s /m l = 28.15, and a physical pseudo-Goldstone pion mass, m π ≃ 135 MeV.
The main observable we are interested in is the Polyakov loop and its dependence on µ B . In particular, as already described above, the ratio of Polyakov loops at different baryon chemical potentials gives access to the µ B -dependent part of the free energy density, and if the ratio is taken for Polyakov loops measured at the same value of the inverse bare coupling β and of the bare quark masses, then no further renormalization is expected, at least when the chemical potential is inserted on the lattice with the prescription introduced in Ref. [17] and adopted in the present investigation. That means that the dependence of ∆F Q (T, µ B , β) on β is expected to be limited to finite UV corrections to continuum scaling. It would be interesting to study the dependence of F Q on µ B in the whole range of physically relevant values of µ B , however our investigation will be limited to the region of small µ B /T and, in particular, to the susceptibility χ Q,µ 2 B defined in Eq. (4), which can be directly related to the Polyakov loop ratio of Eq. (8) by the formula since, from Eq. (8), one has The reason of the limitation to small chemical potentials is the well known sign problem of QCD at finite density, which makes standard Monte-Carlo simulations unfeasible when µ B = 0. Present strategies to partially circumvent the sign problem are reliable only in a limited range of small µ B /T , where they lead to controllable systematic errors; Taylor expansion [25][26][27][28] and analytic continuation from simulations at imaginary chemical potential  are the most widely used techniques. In this investigation we employ both of them, since in part of our wide temperature range the statistical or systematic errors of one technique are less under control, so that a direct comparison with the other technique improves the overall reliability of the results; this combined strategy has revealed successful in other cases, like for the determinations of the curvature of the pseudo-critical line [53].
In the analytic continuation approach, the baryon chemical potential is taken to be purely imaginary, µ B = iµ B,I , the path-integral measure staying real and positive for µ B,I = 0. Within our numerical setup, adding a non-zero µ B,I can be rephrased in terms of a rotation of temporal boundary conditions of the quark fields by a factor exp(iµ I /T ), where µ I = µ B,I /3 is the imaginary part of the quark chemical potential. The value of the Polyakov loop is measured for several values of µ I at fixed temperature, then numerical data are fitted to the analytic continuation of some suitable ansatz for the dependence on µ B , thus fixing the corresponding parameters. Despite its simplicity, this method has some limitations and drawbacks, its systematic errors being related essentially to the arbitrary ansatz for the fitting function.
The choice of the fitting function and the related systematics can be different depending on the value of the temperature, as dictated by the non-trivial symmetries and phase structure of the T −µ B,I phase diagram, which is sketched in Fig. 1. In general one can prove, combining µ B,I translations with gauge field center transformations, that the theory is 2π-periodic in µ B,I /T [54]. This periodicity is smoothly realized for T < T c : there a Fourier expansion is the most natural choice [32] and, moreover, a picture based on the Hadron Resonance Gas (HRG) model suggests an ansatz where the first few terms of the expansion are dominant, unless one is close enough to T c .
On the contrary, at high T , in particular for T > T RW (where T RW ≃ 210 MeV in the continuum limit for N f = 2 + 1 QCD with physical quark masses [55]), the periodicity is realized in a non-analytic way, with first order phase transition lines (RW-lines) crossed for µ B,I /T = (2k + 1)π and k integer: the phase of the Polyakov loop is an order parameter for such transitions, at which the systems switches from one center sector to the other. That limits the range of chemical potentials available for analytic continuation to µ B,I /T < π, however the dependence of the Polyakov loop modulus is well approximated by an even power law expansion in µ B,I , with the lowest order terms becoming more and more dominant as the temperature is increased. The intermediate region, T c < T < T RW , is the one where systematic errors can be more severe. In this region, moving in µ B,I /T from 0 to π one crosses the analytic continuation of the pseudocritical line: even if this is not a true transition but just a crossover, it can make the dependence on µ B,I non-trivial, thus in fact restricting the region of µ B,I where different ansatzs give consistent results; moreover, such a region is smaller and smaller as T c is approached from above.
A second possibility, which can be put in the general framework of the Taylor expansion approach, is to measure χ Q,µ 2 B directly at µ B = 0, following its definition in Eq. (4). In particular, after some computations (which are reported in the appendix), one writes χ Q,µ 2 B as a combination of correlators involving the Polyakov loop and fermionic terms. The expression is where n = n u + n d + n s is the total quark number and n ′ is its derivative with respect to µ B . Even though the measure of this quantity is well defined and seemingly straightforward for all temperatures, in practice its computation involves many noisy estimators and therefore turns out to be numerically expensive, especially in the region around and below T c . In view of the above considerations, the strategy chosen in this work has been to adopt analytic continuation for all temperatures below T c and for most temperatures above T RW , while in the region T c < T < T RW we have adopted both Taylor expansion and analytic continuation, in order to have better control over systematics.  , chosen so as to stay on a line of constant physics at the physical point, using a spline interpolation of the data in Refs. [23,24].
Monte-Carlo simulations have been performed for two different values of N t in order to estimate the impact of UV corrections, in particular on a 24 3 × 6 and on a 32 3 × 8 lattices using a Rational Hybrid Monte-Carlo algorithm [56][57][58]. A summary of the parameters adopted in our simulations, together with details on the strategy chosen in each case, is reported in Tab. II. In the cases in which the susceptibility χ Q,µ 2 B has been mea-sured through Taylor expansion, sets of about 10 4 configurations separated by 10 molecular dynamics trajectories have been analyzed for each run, and fermionic observables such as the quark number n and its derivative n ′ have been computed through stochastic noisy estimators [59], in particular using up to 256 Z 2 random noise vectors per measurement. In the cases in which analytic continuation has been adopted, we have performed around 5 × 10 3 molecular dynamics trajectories for each value of the imaginary chemical potentials. The data analysis has been performed by means of a blocked jackknife resampling in all cases.

III. RESULTS
Let us start by discussing the determination of χ Q,µ 2 B by analytic continuation. As an illustrative example, in Fig. 2 we report the average values of the squared modulus of the Polyakov loop on the 24 3 ×6 lattice as a function of µ B,I and for some of the explored temperatures. For the sake of readability, we have reported separately determinations at high and low T , normalizing data by the value at µ B,I = 0 only in the latter case. At low temperatures, as a matter of fact, we have found that a single cosine term is sufficient to correctly describe our data for all explored temperatures, i.e. with values of the χ 2 /d.o.f. regression parameter close to one: This allows to determine χ Q,µ 2 B . We have considered in the final error also the variability which is obtained by adding a further term in the Fourier expansion 1 , i.e. a term proportional to cos(2 µ B,I /T ).
In the high-temperature regime, instead, we have adopted a polynomial expansion truncated to the quartic term in µ B,I , i.e.
In all cases the fit range has been limited by the location of the pseudocritical value of µ B,I for the given temperature, as extracted from data reported in Refs. [48], and appropriate systematic uncertainties have been added to the fit parameters, which take into account the variability under changes of the fitted range. We have found that the quartic coefficient l 4 is not needed to obtain reasonable fits (and turns out to be compatible with zero when included) for temperatures T > T RW , while for lower temperatures it is definitely needed in order to get In the region above T c , where the pseudo-critical behavior is more pronounced, and in some cases also for the same temperatures at which analytic continuation has been used, we adopted the Taylor expansion method, measuring directly the value χ Q,µB through the formula in Eq. (11). The computation, especially close to T c , turned out to be numerically expensive and, in general, the uncertainties associated to the measures obtained by this method are larger than those extracted by analytic continuation. Nevertheless, in this way no source of systematics is present and, at least at our level of precision, the estimations make the picture clear enough. Moreover, for the temperatures where both methods are available, a reasonable agreement is observed.
The whole collection of results, including all temperatures and both sets of lattice spacings, N t = 6 and N t = 8, is reported in Fig. 3. The dependence on N t appears to be small, confirming that, even if no continuum extrapolation is performed in this study, finite UV cutoff corrections are not large. The susceptibitliy χ Q,µB grows rapidly in the crossover region near T c , where it exhibits a well-defined peak. The location of the peak can be determined quantitatively by modelling the observed behavior near the maximum. In particular, we have adopted a Lorentzian function, defined as where T L indicates the pseudo-critical temperature related to the observable χ Q,µ 2 B . This ansatz well describes the peak structure for both values of N t : best-fit curves are shown in Fig. 4 and yield T L = 143.4 ± 1.2 MeV and T L = 147.7 ± 1.4 MeV respectively for N t = 6 and N t = 8. The uncertainties include systematics related to the choice of the fit range, but not those associated with the determination of the lattice spacing, which are of the order of 2 − 3% [23,24]. Similar results are obtained using a different fitting ansatz, like a purely quadratic function of T . The small N t -dependence observed for T L points to a continuum limit around 150 MeV, which is very close to T c ≃ 155 MeV.

IV. COMPARISON WITH PERTURBATION THEORY
Finally, it is interesting to discuss the fate of χ Q,µ 2 B in the large T limit. At zero baryon chemical po- tential, F Q (T ) is expected to decrease unboundedly as T increases, a well-known behavior predicted by weakcoupling calculations [60,61] and observed also on the lattice in many studies [5,12,15,62]. At leading order, its expression in the high temperature regime is given by where C F = (N 2 c − 1)/2N c is the Casimir operator in the foundamental representation and m D (T ) is the Debye screening mass which, at the leading order is In the dense medium, screening effects are amplified and the value of the single quark free energy grows indefinitely (in module). In the very large temperature limit, at leading order, the expression of F Q (T, µ B ) is obtained performing an expansion of the Debye mass for small values of the chemical potential [3,15]. The result is the appearance of a quadratic dependence on µ B , Inserting this expression in Eq. (4) one finds Consequently, since the coupling runs to zero at large T , the susceptibility χ Q,µ 2 B vanishes asymptotically as g 3 . This means that, in this regime, a finite baryon density does not affect the in-medium static quark free energy, its contribution being overrided by the thermal fluctuations. Notice that the same proportionality to g 3 at high T is shown also by static quark entropy S Q = −∂F Q /∂T which, asymptotically, is expected to behave as S Q ∼ −F Q /T [12,61], in agreement with our calculation.
In order to check the consistency of these predictions with lattice results, we have extended the computation of χ Q,µ 2 B to higher temperatures, adopting the Taylor expansion method which in this regime is not particularly expensive. Results are shown in Fig. 5. In order to obtain a quantitative prediction from Eq. (18), we need to insert the dependence of the coupling constant g(T ) on the temperature, which at the leading order in perturbation theory is given by [63] g −2 (T ) = 2β 0 log 2πT where β 0 is the first coefficient of the QCD β-function, which is independent of the renormalization scheme. Inserting this expression in Eq. (18) one obtains where p 0 is a pre-factor which is independent of the renormalization scheme and whose value is p 0 ∼ 0.019 in our case, where N c = 3 and N f = 3 (we assume that the three quark flavors can be considered as practically degenerate in this temperature regime). The slow decrease shown by the lattice data is well described, both for N t = 6 and N t = 8, by Eq. (20), the fitted value of p 0 being 0.021(2) and 0.014(3), respectively, for the 24 3 × 6 (χ 2 /d.o.f. = 0.73) and the 32 3 × 8 (χ 2 /d.o.f. = 0.46) lattice: we consider such an agreement more than satisfactory, given that only the leading order has been considered; it is interesting to notice that also the values obtained for the Λ parameter are reasonable and of the order of 100 MeV.

V. CONCLUSIONS
In this study we have investigated the dependence of the static quark free energy on the baryon chemical potential in a wide temperature range, considering in particular the leading order dependence, which is quadratic in µ B and that we have parameterized in terms of the susceptibility χ Q,µ 2 B . The investigation has been carried out by lattice simulations of N f = 2 + 1 QCD discretized via stout-staggered fermions with physical quark masses. Both analytic continuation and Taylor expansion have been adopted to avoid the sign problem at non-zero µ B , obtaining consistent results.
Results for χ Q,µ 2 B have been found to be compatible, in the high temperature regime, with predictions obtained in perturbation theory. The dependence of the static quark free energy on µ B which vanishes as a power law in the gauge coupling g(T ), precisely as g 3 , i.e. logarithmically with the temperature T . Numerical results are consistent both with the power law behavior in g and with the predicted prefactor.
At low temperatures χ Q,µ 2 B presents instead a well defined peak located around 150 MeV, i.e. roughly compatible with the crossover temperature T c corresponding to the restoration of chiral symmetry. If the Polyakov loop were an exact order parameter for the deconfinement transition, one would expect a singular behavior for χ Q,µ 2 B at the critical temperature. Therefore, the rough coincidence of the two temperatures points once again to a strong connection between chiral symmetry and deconfinement dynamics, even within a crossover scenario.
Our results have been obtained for just two sets of lattice spacings, corresponding to N t = 6 and N t = 8. Future studies should extend the investigation to larger values of N t so as to achieve a continuum extrapolation for χ Q,µ 2 B . However, present results show only modest changes as N t is changed from 6 to 8, so that no significant modifications of our conclusions are expected in the continuum limit. The expression of the curvature χ Q,µ 2 B is obtained by computing the second derivative of the ratio between square modules of the Polyakov loop, as in Eq. (10). Applying the derivative operator ∂ µ ≡ ∂/∂(µ/T ) to the numerator, which is the only part depending on the chemical potential, one has ∂ 2 µ TrL 2 = 2 (∂ µ ReTrL ) 2 + 2 ReTrL ∂ 2 µ ReTrL + ReTrL ↔ ImTrL , where a similar expression holds for ImTrL and Z is the partition function defined in Eq. (5). Since the Polyakov loop does not depend explicitly on the chemical potential, all dependence on µ is carried by the Dirac matrix. That means that the derivative operator will act only on the fermionic part of the functional integral, which appears also in the denominator. One has and the same is true also for ImTrL . Further application of the derivative ∂ µ leads to new correlators involving the quark number n or its derivative n ′ = ∂ µ n. Indeed, one finds that ∂ µ n ReTrL = n 2 ReTrL − n ReTrL n + n ′ ReTrL ∂ µ n = n 2 − n 2 + n ′ , (A. 6) where n ′ = f n ′ f and n ′ f = ∂ µ n f with Finally, joining and re-arranging all the pieces appearing in Eq. (A.1), the following expression is found ∂ 2 µ | TrL | 2 = 2 n ReTrL 2 + 6 ReTrL 2 n 2 − 8 ReTrL n ReTrL n + 2 ReTrL n 2 ReTrL − 2 ReTrL 2 n 2 + 2 ReTrL n ′ ReTrL − 2 ReTrL 2 n ′ + {ReTrL ↔ ImTrL} . (A.8) The curvature χ Q,µB is obtained by normalizing this formula with the square module of TrL(0) and evaluating the ratio at zero chemical potential, see Eq. (10). As a result, the expression above simplifies since, for µ = 0, both the quark number n and ImTrL vanish because of charge conjugation symmetry. Then, re-arranging the remaining terms the definition in Eq. (11) is found.