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

TL;DR

This study investigates how the static quark free energy varies with baryon chemical potential in $N_f=2+1$ QCD across a broad temperature range, revealing sensitivity to the chiral crossover and aligning with perturbative predictions at high temperatures.

Contribution

It provides the first detailed analysis of the dependence of static quark free energy on $$ in $N_f=2+1$ QCD using both Taylor expansion and analytic continuation methods.

Findings

The $$-susceptibility peaks around 150 MeV, indicating sensitivity to the chiral crossover.

Results are consistent between Taylor expansion and analytic continuation approaches.

High-temperature behavior matches perturbative theoretical expectations.

Abstract

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.

Tables

Table 1. Table 1: List of parameters used in the Monte-Carlo simulations for the study of the susceptibility χ Q , μ B 2 subscript 𝜒 𝑄 superscript subscript 𝜇 𝐵 2 \chi_{Q,\mu_{B}^{2}} , chosen so as to stay on a line of constant physics at the physical point, using a spline interpolation of the data in Refs.  Borsanyi:2010cj ; Borsanyi:2013bia .

$N^3\times N_t$	$\beta$	$a\left[\mathrm{fm}\right]$	$T\left[\mathrm{MeV}\right]$	${\mu }_I/\left(\pi T\right)$
${24}^3\times 6$	3.4500	0.2835	116	$0,0.04,\dots ,0.32$
”	3.4789	0.2631	125	$0,0.04,\dots ,0.32$
”	3.5085	0.2436	135	$0,0.04,\dots ,0.32$
”	3.5246	0.2332	141	$0,0.04,\dots ,0.64$
”	3.5421	0.2222	148	$0$
”	3.5585	0.2121	155	$0$
”	3.5695	0.2055	160	$0$
”	3.5800	0.1993	165	$0,0.04,\dots ,0.32$
”	3.5923	0.1923	171	$0$
”	3.6172	0.1787	184	$0,0.04,\dots ,0.32$
”	3.6746	0.1515	217	$0,0.04,\dots ,0.32$
”	3.7305	0.1310	251	$0,0.04,\dots ,0.32$
”	3.7829	0.1153	285	$0,0.04,\dots ,0.32$
”	3.8300	0.1034	318	$0,0.04,\dots ,0.32$
”	3.8749	0.0936	351	$0,0.04,\dots ,0.32$
”	3.9184	0.0856	384	$0,0.04,\dots ,0.32$
”	3.9608	0.0788	417	$0,0.04,\dots ,0.32$
”	4.0019	0.0729	451	$0,0.04,\dots ,0.32$
”	4.0798	0.0635	518	$0$
”	4.1506	0.5622	585	$0$
”	4.2200	0.0504	652	$0$
”	4.2797	0.4574	719	$0$
”	4.3297	0.0418	786	$0$
”	4.3778	0.0386	853	$0$
”	4.4284	0.0357	920	$0$
”	4.4808	0.0333	987	$0$
”	4.5317	0.0312	1054	$0$
”	4.5764	0.0293	1121	$0$
${32}^3\times 8$	3.5835	0.1973	125	$0,0.04,\dots ,0.32$
”	3.6100	0.1827	135	$0,0.04,\dots ,0.64$
”	3.6245	0.1749	141	$0,0.04,\dots ,0.64$
”	3.6417	0.1666	148	$0,0.04,\dots ,0.64$
”	3.6570	0.1591	155	$0$
”	3.6700	0.1541	160	$0,0.02,\dots ,0.16$
”	3.6800	0.1494	165	$0,0.02,\dots ,0.24$
”	3.6925	0.1442	171	$0,0.04,\dots ,0.32$
”	3.7250	0.1333	185	$0$
”	3.8525	0.0982	251	$0$
”	4.1678	0.0546	451	$0,0.04,\dots ,0.32$
”	4.2560	0.0476	518	$0$
”	4.3255	0.0422	585	$0$
”	4.3899	0.0378	652	$0$
”	4.4586	0.0343	719	$0$
”	4.5273	0.0314	786	$0$
”	4.5861	0.0289	853	$0$

Equations

$$L(\mathbf{r})=\frac{1}{N_c}\mathcal{P}\exp \left(ig{\int }_0^{1/T}d\tau A_0(\mathbf{r},\tau )\right)\text{leavevmode},$$

$$2F_{Q}=-{T}\log|\left\langle\mathrm{Tr}L\right\rangle|^{2}\,.\leavevmode\nobreak\$$

$$\frac{\Delta F_Q(T,{\mu }_B)}{T}=-\log \left(\frac{|⟨\mathrm{Tr}L⟩(T,{\mu }_B)|}{|⟨\mathrm{Tr}L⟩(T,\text{leavevmode}0\text{leavevmode})|}\right),$$

$$\chi_{Q,\mu_{B}^{2}}\equiv-\frac{\partial^{2}(F_{Q}/T)}{\partial(\mu_{B}/T)^{2}}\bigg{|}_{\mu_{B}=0}\,.$$

$$\mathcal{Z}=\int \mathcal{D}U{e}^{-{\mathcal{S}}_{\text{YM}}}\prod _{f=u,d,s}\det {\left[M_{\text{st}}^f(U,{\mu }_f)\right]}^{1/4}\text{leavevmode},$$

$${\mathcal{S}}_{\text{YM}}=-\frac{\beta }{3}\sum _{i,\mu \ne \nu }\left(\frac{5}{6}{W}_{i;\mu \nu }^{1\times 1}-\frac{1}{12}{W}_{i;\mu \nu }^{1\times 2}\right)$$

$$M_{\text{st}}^f(U,{\mu }_f)$$

$$\frac{\Delta F_Q(T,{\mu }_B,\beta )}{T}=-\log \left(\frac{|⟨\mathrm{Tr}L⟩(T,{\mu }_B,\beta )|}{|⟨\mathrm{Tr}L⟩(T,\text{leavevmode}0\text{leavevmode},\beta )|}\right)\text{leavevmode},$$

$$\frac{|⟨\mathrm{Tr}L⟩(T,{\mu }_B)|}{|⟨\mathrm{Tr}L⟩(T,\text{leavevmode}0\text{leavevmode})|}=1+\frac{1}{2}{\chi }_{Q,{\mu }_B^2}{\left(\frac{{\mu }_B}{T}\right)}^2+\mathcal{O}\left({\left(\frac{{\mu }_B}{T}\right)}^4\right)$$

$$\frac{\partial^{2}}{\partial(\mu_{B}/T)^{2}}\frac{|\langle\textnormal{Tr}L\rangle(T,\mu_{B})|}{|\langle\textnormal{Tr}L\rangle(T,\leavevmode\nobreak\ 0\leavevmode\nobreak\ )|}\bigg{|}_{\mu_{B}=0}=-\frac{\partial^{2}(F_{Q}/T)}{\partial(\mu_{B}/T)^{2}}\bigg{|}_{\mu_{B}=0}$$

$${\chi }_{Q,{\mu }_B^2}$$

$$\frac{|⟨L⟩({\mu }_{B,I}){|}^2}{|⟨L⟩(0){|}^2}=1-2{\chi }_{Q,{\mu }_B^2}\left[1-\cos \left(\frac{{\mu }_{B,I}}{T}\right)\right].$$

$$\frac{|⟨L⟩({\mu }_{B,I}){|}^2}{|⟨L⟩(0){|}^2}=1-{\chi }_{Q,{\mu }_B^2}{\left(\frac{{\mu }_{B,I}}{T}\right)}^2+l_4{\left(\frac{{\mu }_{B,I}}{T}\right)}^4.$$

$${\chi }_{Q,{\mu }_B^2}=\frac{p_0}{1+{\left[(T-T_L)/p_1\right]}^2},$$

$$F_Q(T)=-\frac{C_F}{2}\frac{g^2}{4\pi }{m}_D(T)\text{leavevmode},$$

$$m_D^2(T)=\frac{1}{3}\left(N_c+\frac{N_f}{2}\right)g^2{T}^2\text{leavevmode}.$$

$$m_D^2(T,{\mu }_B)=m_D^2(T)\left[1+\frac{3N_f}{2N_c+N_f}{\left(\frac{{\mu }_B}{3\pi T}\right)}^2\right]\text{leavevmode},$$

$$\displaystyle\chi_{Q,\mu_{B}^{2}}\big{|}_{T\to\infty}$$

$$g^{-2}(T)=2{\beta }_0\log \frac{2\pi T}{\Lambda }{\beta }_0=\frac{11N_f-2N_c}{48{\pi }^2}\text{leavevmode},$$

$$\chi_{Q,\mu_{B}^{2}}\big{|}_{T\to\infty}=p_{0}\left[\log\frac{2\pi T}{\Lambda}\right]^{-3/2}\leavevmode\nobreak\ ,$$

$$\displaystyle\partial_{\mu}^{2}\big{|}\langle\mathrm{Tr}L\rangle\big{|}^{2}$$

$$⟨\mathrm{ReTr}L⟩=\frac{1}{\mathcal{Z}}\int \mathcal{D}U{e}^{-{\mathcal{S}}_{\text{YM}}}\mathrm{ReTr}L\prod _f\det {\left[M_{\text{st}}^f\right]}^{\frac{1}{4}},\text{vspace}-0.3em$$

$$\partial_{\mu}\prod_{f}\det\left[M_{\text{st}}^{f}\right]^{\frac{1}{4}}=\bigg{(}\sum_{f}n_{f}\bigg{)}\prod_{f}\det\left[M_{\text{st}}^{f}\right]^{\frac{1}{4}},$$

$$n_f=\frac{1}{4}\mathrm{Tr}\left[{M_{\text{st}}^f}^{-1}{\partial }_{\mu }{M}_{\text{st}}^f\right].$$

$${\partial }_{\mu }⟨\mathrm{ReTr}L⟩=⟨n\mathrm{ReTr}L⟩-⟨\mathrm{ReTr}L⟩⟨n⟩$$

$${\partial }_{\mu }⟨n\mathrm{ReTr}L⟩$$

$${\partial }_{\mu }⟨n⟩$$

$${\partial }_{\mu }{n}_f=\frac{1}{4}\mathrm{Tr}\left[{\left({M_{\text{st}}^f}^{-1}{\partial }_{\mu }{M}_{\text{st}}^f\right)}^2-{M_{\text{st}}^f}^{-1}{\partial }_{\mu }^2{M}_{\text{st}}^f\right]\text{leavevmode}.$$

$${\partial }_{\mu }^2|⟨\text{Tr}L⟩{|}^2$$

Full text

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 Kaczmarek:2002mc ; Petreczky:2004pz ; Kaczmarek:2007pb ; Kaczmarek:2005gi ; Borsanyi:2015yka , from the expectation value of the Polyakov loop and of its correlators. The Polyakov loop is defined in the continuum as

[TABLE]

where $T$ is the temperature, $\mathcal{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

[TABLE]

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 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 aefks ; afks ; betal ; tchot ; tchot2 ; tchot3 . 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 Bazavov:2016uvm , while other related observables show a signal closer to the chiral transition temperature: this is the case for the Polyakov loop entropy $S_{Q}=-\partial F_{Q}/\partial T$ Bazavov:2016uvm or the so-called transverse susceptibility related to fluctuations in the imaginary part of the Polyakov loop Lo:2013etb ; Lo:2013hla . 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 $\mu_{B}$. The modification of the heavy quark free energy due to $\mu_{B}$, $\Delta F_{Q}(T,\mu_{B})\equiv F_{Q}(T,\mu_{B})-F_{Q}(T,0)$, is given by the following expression

[TABLE]

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. Doring:2005ih and more recently in Ref. Andreoli:2017zie for QCD with physical quark masses.

One expects the dependence of $\Delta F_{Q}(T,\mu_{B})$ on $\mu_{B}$ to be sensitive to the location of the transition. Indeed, if the Polyakov loop were an exact order parameter then its dependence on $\mu_{B}$ should become singular at $T_{c}$, because $\mu_{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 $\mu_{B}$, the first non-trivial derivative to investigate the associate pseudocritical behavior is the mixed susceptibility

[TABLE]

Early simulations of $N_{f}=2$ QCD have shown that this quantity has a broad peak in a region close to $T_{c}$ Doring:2005ih . More recent simulations, performed for $N_{f}=2+1$ QCD discretized via stout-staggered fermions with physical quark masses Andreoli:2017zie , were limited to a temperature range $T\leavevmode\nobreak\ \gtrsim\leavevmode\nobreak\ 180$ MeV, showing nevertheless a peculiar behavior pointing to a seeming divergence for $T\sim 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. Andreoli:2017zie 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 $\chi_{Q,\mu_{B}^{2}}$ around $T_{c}$, and to compare results obtained at high $T$ with perturbative predictions. Since lattice simulations at non-zero $\mu_{B}$ are not feasible, because of the sign problem, we employ both Taylor expansion and analytic continuation from simulations at imaginary $\mu_{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 $\mu_{f}$ ($f=u,d,s$) coupled to quark number operators, ${\mathcal{Z}}(T,\mu_{u},\mu_{d},\mu_{s})$, in a setup for which $\mu_{u}=\mu_{d}=\mu_{s}=\mu_{B}/3$, corresponding to a purely baryonic chemical potential. The path integral formulation of ${\mathcal{Z}}(T,\mu_{B})$, discretized via improved rooted staggered fermions and adopting the standard exponentiated implementation of the chemical potentials Hasenfratz:1983ba ; Gavai:1985ie , reads

[TABLE]

where

[TABLE]

is the tree-level Symanzik improved action Weisz:1982zw ; Curci:1983an ($W_{i;\mu\nu}^{n\times m}$ stands for the trace of the $n\times m$ rectangular parallel transport in the $\mu$-$\nu$ plane and starting from site $i$), and the staggered fermion matrix is defined as

[TABLE]

where $U_{i;\nu}^{(2)}$ are two-times stout-smeared links, with isotropic smearing parameter $\rho=0.15$ Morningstar:2003gk . Bare parameters have been set so as to stay on a line of constant physics Aoki:2009sc ; Borsanyi:2010cj ; Borsanyi:2013bia , 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_{\pi}\simeq 135\leavevmode\nobreak\ \text{MeV}$.

The main observable we are interested in is the Polyakov loop and its dependence on $\mu_{B}$. In particular, as already described above, the ratio of Polyakov loops at different baryon chemical potentials gives access to the $\mu_{B}$-dependent part of the free energy density, $\Delta F_{Q}(T,\mu_{B})\equiv F_{Q}(T,\mu_{B})-F_{Q}(T,0)$,

[TABLE]

and if the ratio is taken for Polyakov loops measured at the same value of the inverse bare coupling $\beta$ 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. Hasenfratz:1983ba and adopted in the present investigation. That means that the dependence of $\Delta F_{Q}(T,\mu_{B},\beta)$ on $\beta$ is expected to be limited to finite UV corrections to continuum scaling.

It would be interesting to study the dependence of $F_{Q}$ on $\mu_{B}$ in the whole range of physically relevant values of $\mu_{B}$, however our investigation will be limited to the region of small $\mu_{B}/T$ and, in particular, to the susceptibility $\chi_{Q,\mu_{B}^{2}}$ defined in Eq. (4), which can be directly related to the Polyakov loop ratio of Eq. (8) by the formula

[TABLE]

since, from Eq. (8), one has

[TABLE]

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 $\mu_{B}\neq 0$. Present strategies to partially circumvent the sign problem are reliable only in a limited range of small $\mu_{B}/T$, where they lead to controllable systematic errors; Taylor expansion tay1 ; tay2 ; tay3 ; tay4 and analytic continuation from simulations at imaginary chemical potential alford ; lomb99 ; fp1 ; dl1 ; azcoiti ; chen ; Wu:2006su ; NN2011 ; giudice ; ddl07 ; cea2009 ; alexandru ; cea2012 ; Karbstein:2006er ; cea_other ; sanfo1 ; Takaishi:2010kc ; cea_hisq1 ; corvo ; nf2BFEPS ; bellwied ; gunther ; gagliardi ; Bornyakov:2017upg 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 Bonati:2018nut .

In the analytic continuation approach, the baryon chemical potential is taken to be purely imaginary, $\mu_{B}=i\mu_{B,I}$, the path-integral measure staying real and positive for $\mu_{B,I}\neq 0$. Within our numerical setup, adding a non-zero $\mu_{B,I}$ can be rephrased in terms of a rotation of temporal boundary conditions of the quark fields by a factor $\exp(i\mu_{I}/T)$, where $\mu_{I}=\mu_{B,I}/3$ is the imaginary part of the quark chemical potential. The value of the Polyakov loop is measured for several values of $\mu_{I}$ at fixed temperature, then numerical data are fitted to the analytic continuation of some suitable ansatz for the dependence on $\mu_{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-\mu_{B,I}$ phase diagram, which is sketched in Fig. 1. In general one can prove, combining $\mu_{B,I}$ translations with gauge field center transformations, that the theory is $2\pi$-periodic in $\mu_{B,I}/T$ rwpaper . This periodicity is smoothly realized for $T<T_{c}$: there a Fourier expansion is the most natural choice dl1 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}\simeq 210$ MeV in the continuum limit for $N_{f}=2+1$ QCD with physical quark masses Bonati:2016pwz ), the periodicity is realized in a non-analytic way, with first order phase transition lines (RW-lines) crossed for $\mu_{B,I}/T=(2k+1)\pi$ 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 $\mu_{B,I}/T<\pi$, however the dependence of the Polyakov loop modulus is well approximated by an even power law expansion in $\mu_{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 $\mu_{B,I}/T$ from 0 to $\pi$ 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 $\mu_{B,I}$ non-trivial, thus in fact restricting the region of $\mu_{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 $\chi_{Q,\mu_{B}^{2}}$ directly at $\mu_{B}=0$, following its definition in Eq. (4). In particular, after some computations (which are reported in the appendix), one writes $\chi_{Q,\mu_{B}^{2}}$ as a combination of correlators involving the Polyakov loop and fermionic terms. The expression is

[TABLE]

where $n=n_{u}+n_{d}+n_{s}$ is the total quark number and $n^{\prime}$ is its derivative with respect to $\mu_{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.

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}\times 6}$ and on a ${32^{3}\times 8}$ lattices using a Rational Hybrid Monte-Carlo algorithm rhmc1 ; rhmc2 ; rhmc3 . A summary of the parameters adopted in our simulations, together with details on the strategy chosen in each case, is reported in Tab. 1.

In the cases in which the susceptibility $\chi_{Q,\mu_{B}^{2}}$ has been measured 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^{\prime}$ have been computed through stochastic noisy estimators Dong:1993pk , 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\times 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 $\chi_{Q,\mu_{B}^{2}}$ 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}\times 6$ lattice as a function of $\mu_{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 $\mu_{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 $\chi^{2}/{\rm d.o.f.}$ regression parameter close to one:

[TABLE]

This allows to determine $\chi_{Q,\mu_{B}^{2}}$. We have considered in the final error also the variability which is obtained by adding a further term in the Fourier expansion111 Notice that the parametrization in Eq. (12) changes if other Fourier terms are added, since in this case $\chi_{Q,\mu_{B}^{2}}$ takes contributions from all Fourier coefficients., i.e. a term proportional to $\cos({2\,\mu_{B,I}}/{T})$.

In the high-temperature regime, instead, we have adopted a polynomial expansion truncated to the quartic term in $\mu_{B,I}$, i.e.

[TABLE]

In all cases the fit range has been limited by the location of the pseudocritical value of $\mu_{B,I}$ for the given temperature, as extracted from data reported in Refs. corvo , 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 $\chi^{2}/{\rm d.o.f.}\sim 1$.

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 $\chi_{Q,\mu_{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 $\chi_{Q,\mu_{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

[TABLE]

where $T_{L}$ indicates the pseudo-critical temperature related to the observable $\chi_{Q,\mu_{B}^{2}}$. 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\pm 1.2\leavevmode\nobreak\ \textnormal{MeV}$ and $T_{L}=147.7\pm 1.4\leavevmode\nobreak\ \textnormal{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\%$ Borsanyi:2010cj ; Borsanyi:2013bia . 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}\simeq 155$ MeV.

IV Comparison with perturbation theory

Finally, it is interesting to discuss the fate of $\chi_{Q,\mu_{B}^{2}}$ in the large $T$ limit. At zero baryon chemical potential, $F_{Q}(T)$ is expected to decrease unboundedly as $T$ increases, a well-known behavior predicted by weak-coupling calculations Gava:1981qd ; Berwein:2015ayt and observed also on the lattice in many studies Doring:2005ih ; Bazavov:2013yv ; Borsanyi:2015yka ; Bazavov:2016uvm . At leading order, its expression in the high temperature regime is given by

[TABLE]

where $C_{F}=(N_{c}^{2}-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

[TABLE]

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,\mu_{B})$ is obtained performing an expansion of the Debye mass for small values of the chemical potential Doring:2005ih ; Kaczmarek:2007pb . The result is the appearance of a quadratic dependence on $\mu_{B}$, $F_{Q}(T,\mu_{B})=F_{Q}(T)\,m_{D}(T,\mu_{B})/m_{D}(T)$ where

[TABLE]

Inserting this expression in Eq. (4) one finds

[TABLE]

Consequently, since the coupling runs to zero at large $T$, the susceptibility $\chi_{Q,\mu_{B}^{2}}$ 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}=-\partial F_{Q}/\partial T$ which, asymptotically, is expected to behave as $S_{Q}\sim-F_{Q}/T$ Berwein:2015ayt ; Bazavov:2016uvm , in agreement with our calculation.

In order to check the consistency of these predictions with lattice results, we have extended the computation of $\chi_{Q,\mu_{B}^{2}}$ 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 Kaczmarek:2004gv

[TABLE]

where $\beta_{0}$ is the first coefficient of the QCD $\beta$-function, which is independent of the renormalization scheme. Inserting this expression in Eq. (18) one obtains

[TABLE]

where $p_{0}$ is a pre-factor which is independent of the renormalization scheme and whose value is $p_{0}\sim 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}\times 6$ ($\chi^{2}/{\rm d.o.f.}=0.73$) and the $32^{3}\times 8$ ($\chi^{2}/{\rm 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 $\Lambda$ 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 $\mu_{B}$ and that we have parameterized in terms of the susceptibility $\chi_{Q,\mu_{B}^{2}}$. 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 $\mu_{B}$, obtaining consistent results.

Results for $\chi_{Q,\mu_{B}^{2}}$ 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 $\mu_{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 $\chi_{Q,\mu_{B}^{2}}$ 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 $\chi_{Q,\mu_{B}^{2}}$ 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 $\chi_{Q,\mu_{B}^{2}}$. 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.

Acknowledgements.

Numerical simulations have been performed on the MARCONI machine at CINECA, based on the agreement between INFN and CINECA (under project INF18_npqcd), at the Scientific Computing Center at INFN-PISA, and on the COKA cluster at the University of Ferrara and INFN-Ferrara based on the GPU code developed in Refs. incardona ; ferrarapisa ; ferrarapisa2 . FN acknowledges financial support from the INFN HPC_HTC project. *

Appendix A Computation of $\chi_{Q,\mu_{B}^{2}}$

The expression of the curvature $\chi_{Q,\mu_{B}^{2}}$ 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 $\partial_{\mu}\equiv\partial/\partial(\mu/T)$ to the numerator, which is the only part depending on the chemical potential, one has

[TABLE]

where $\mu=\mu_{B}/3$ is the common chemical potential for all flavors, and the last line in brackets indicates terms where real and imaginary parts of the Polyakov loop are exchanged. The expectation values entering this expression can be written as

[TABLE]

where a similar expression holds for $\langle\mathrm{ImTr}L\rangle$ and $\mathcal{Z}$ is the partition function defined in Eq. (5). Since the Polyakov loop does not depend explicitly on the chemical potential, all dependence on $\mu$ 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

[TABLE]

where $n_{f}$ is the quark number operators related to each different flavor,

[TABLE]

Setting $n=\sum_{f}n_{f}$ one can rewrite the derivative of the expression in Eq. (22) as

[TABLE]

and the same is true also for $\langle\mathrm{ImTr}L\rangle$. Further application of the derivative $\partial_{\mu}$ leads to new correlators involving the quark number $n$ or its derivative $n^{\prime}=\partial_{\mu}n$. Indeed, one finds that

[TABLE]

where $n^{\prime}=\sum_{f}n_{f}^{\prime}$ and $n_{f}^{\prime}=\partial_{\mu}n_{f}$ with

[TABLE]

Finally, joining and re-arranging all the pieces appearing in Eq. (21), the following expression is found

[TABLE]

The curvature $\chi_{Q,\mu_{B}}$ is obtained by normalizing this formula with the square module of $\langle\mathrm{Tr}L(0)\rangle$ and evaluating the ratio at zero chemical potential, see Eq. (10). As a result, the expression above simplifies since, for $\mu=0$, both the quark number $\langle n\rangle$ and $\langle\mathrm{ImTr}L\rangle$ vanish because of charge conjugation symmetry. Then, re-arranging the remaining terms the definition in Eq. (11) is found.