title: Sensitivity of the Polyakov loop and related observables to chiral symmetry restoration
authors: David Anthony Clarke Olaf Kaczmarek Frithjof Karsch Anirban Lahiri Mugdha Sarkar
text: 
  
  
  
  
  
   
  
  
  
  
  
  
  Fakultät für Physik, Universität Bielefeld, D-33615
  Bielefeld, Germany
  
   Key Laboratory of Quark & Lepton Physics (MOE) and Institute
  of Particle Physics, Central China Normal University, Wuhan 430079, China
  Fakultät für Physik, Universität Bielefeld, D-33615
  Bielefeld, Germany
  
  Fakultät für Physik, Universität Bielefeld, D-33615
  Bielefeld, Germany
  
  Fakultät für Physik, Universität Bielefeld, D-33615
  Bielefeld, Germany
  
  Fakultät für Physik, Universität Bielefeld, D-33615
  Bielefeld, Germany
  
  
  
  
  
  
  
  
  
  While the Polyakov loop is an order parameter of the deconfinement transition
  in the heavy quark mass regime of QCD, its sensitivity to the deconfinement of
  light, dynamical quarks in QCD is not apparent. On the other hand, the 
  quark mass dependence of the Polyakov loop is sensitive to the appearance
  of a chiral phase transition. 
  Using lattice QCD calculations in the staggered fermion discretization 
  scheme at finite values of the lattice spacing, aT=1/8, we
  show here, for the first time, that the Polyakov loop expectation
  value, and the heavy quark free energy extracted from it, behave like 
  energy-like observables in the vicinity of the chiral phase
  transition temperature T_c. Consistent with scaling behavior of 
  energy-like observables in the 3-d, Ø(2) universality class,
  the quark mass derivatives  
  diverge in the chiral limit at T_c while the temperature derivatives
  stay finite. The latter will develop a characteristic 
  spike at T_c. This, however, may be resolved only in calculations with quark masses being two orders of magnitude
  smaller than those currently accessible in lattice QCD calculations. 
  
  
  
  11.10.Wx, 11.15.Ha, 12.38.Aw, 12.38.Gc, 12.38.Mh, 24.60.Ky, 25.75.Gz, 25.75.Nq
  
  Sensitivity of the Polyakov loop and related observables to chiral symmetry restoration
      Mugdha Sarkar
      December 30, 2023
  =======================================================================================
  
  
  
  
  
  Introduction.—
  The Lagrangian of quantum chromodynamics (QCD), the theory describing
  interactions controlled by the strong force, possesses 
  exact global symmetries only in the massless (chiral) and infinite quark mass 
  (pure gauge) limits. 
  The latter case has been extensively exploited 
  in lattice QCD calculations 
  to discuss the deconfinement
  phase transition in pure gauge theories 
  <cit.> and its imprint 
  in the heavy quark
  sector of QCD. Also in the light quark mass region the rapid change 
  of the Polyakov loop expectation value, P, as function of the 
  temperature, characterized by an inflection point in its T-dependence,
  is often taken as an indication for the occurrence of deconfinement.
  However, studies with (almost) physical light up and down quark masses and
  
  improved discretization schemes for the QCD Lagrangian,
  performed closer to the continuum limit,
  in general show that the QCD transition is a smooth crossover, and no evidence
  for an inflection point in the vicinity of the chiral transition
  temperature is found
  <cit.>.
  
  
  In the limit of  vanishing values of the two light quark masses 
  the chiral flavor symmetry, (2)_L×(2)_R, gets restored above a 
  temperature T_c, giving rise to a chiral phase transition 
  <cit.>.
  In QCD the non-zero light quark 
  masses, m_l, are small on the scale of relevant temperatures, 
  e.g. T_c. 
  Thus the light quark chiral condensate, ψ̅ψ, is a good 
  indicator for the occurrence of a phase transition in the chiral limit of QCD. 
  The maxima in either the quark mass or temperature derivatives of
  ψ̅ψ diverge in the chiral limit, and for 
  m_l>0 the positions of these maxima define pseudo-critical 
  temperatures that converge to T_c in the chiral limit. 
  
  The Polyakov loop is a purely gluonic observable that is 
  trivially invariant under chiral transformations in the fermion sector of the 
  QCD Lagrangian. As far as critical behavior close to a second 
  order phase transition point is concerned, it thus may be expected that 
  P as well as the heavy quark free energy, F_q=-TlnP,
  behave like any other “energy-like" operator that may appear in an
  effective Hamiltonian describing e.g. QCD thermodynamics in the
  vicinity of the chiral phase transition. One may thus expect that 
  P as well as F_q/T are sensitive to critical 
  behavior arising from this transition. 
  
  Here we present results on the temperature and quark mass dependence of 
  P and F_q/T close to the chiral limit. We show that both of them
  reflect properties of “energy-like" observables in the vicinity of T_c
  and discuss the resulting chiral limit behavior. 
  As we will perform calculations at a non-zero value of the lattice 
  spacing, using the Highly Improved Staggered Quark (HISQ) action, the relevant 
  symmetry group for the discussion of universal scaling properties is the 
  Ø(2) rather than Ø(4) group as it will be the case in the 
  continuum limit of 
  lattice QCD. However, as will become clear in the following, none of the
  qualitative features that will arise from the energy-like behavior of the
  Polyakov loop in theories with global Ø(N) symmetry will depend on this
  difference.
  
  
  
  
  
  
  Polyakov loop and the heavy quark free energy.—
  For lattice QCD in a finite Euclidean space-time volume 
  N_σ^3× N_τ, the Polyakov loop, P_x⃗, and its
  spatial average, P, 
  
      _x⃗≡1/3∏_τ U_4(x⃗,τ) 
        , 
        P≡1/N_σ^3 e^N_τ c(g^2)∑_x⃗_x⃗,
  
  are given in terms of (3)-valued field variables, 
  U_4(x⃗,τ), 
  defined on the temporal link (μ= 4) originating at a Euclidean 
  space-time point (x⃗,τ). The volume V=(N_σ a)^3 and inverse 
  temperature  T^-1=N_τ a are given in terms of the lattice spacing a.
  The bare Polyakov loop, , has been renormalized
  using renormalization constants, c(g^2),
  determined in Ref. <cit.> (Table V) for the regularization
  scheme used also in this work.
  
  The heavy quark free energy, F_q/T, 
  characterizes the behavior of correlation functions between static quark 
  and anti-quark sources at infinite distances <cit.>,
  
      F_q(T,H)     =    -TlnP 
      = -T/2lim_|x⃗-y⃗|→∞ln⟨ P^†_x⃗ P^†_y⃗⟩  .
  
  Here H=m_l/m_s parametrizes the quark mass dependence in terms
  of the ratio of degenerate light quark masses m_l≡ m_u=m_d and 
  the strange quark mass m_s.
  For the analysis of the quark mass dependence of the heavy quark free
  energy, we 
  calculate the re­normalization-scheme-independent mixed susceptibility
  
      ∂ F_q(T,H)/T/∂ H   =    -
      1/P∂P/∂ H≡
      - χ_mP/P  ,
  
  with the quark mass derivative of P given by
  
      χ_mP   ≡   ∂P/∂ H
      =   P  ·Ψ - PΨ  .
  
  Here Ψ≡1/2m̂_s  M_l^-1 denotes the extensive 
  observable
  defining a dimensionless combination of the 2-flavor light quark chiral 
  condensate in terms of the light quark, staggered fermion matrix M_l;  
  m̂_s is the bare strange quark mass in lattice units, 
  m_s/T= m̂_s N_τ.
  
  Using fit results for F_q/T and its derivative with respect to the 
  quark mass we will also be able to determine the derivatives of 
  P and F_q/T with respect to T. 
  Similar to the H-derivatives, these derivatives are closely related to each other,
  
      T_c ∂ F_q(T,H)/T/∂ T   =    -T_c/P∂P/∂ T  .
  
  
  
  
  
  Polyakov loop and chiral symmetry restoration.—
  Within Wilson's renormalization group approach 
  <cit.>, thermodynamics in the vicinity of a critical
  point can be described by an effective Hamiltonian, which is defined in a
  multi-dimensional space of operators (observables). These operators may be
  invariant under the global symmetry that gets broken at the critical point or 
  may break this symmetry explicitly. In the former case the operator is said to 
  be energy-like, while in the latter case it is magnetization-like. 
  In QCD the 2-flavor, light quark chiral
  condensate is a typical magnetization-like operator.
  
  
  The Polyakov loop is invariant 
  under chiral transformations of the quark fields. Its expectation value,
  P, as well as the heavy quark free energy, F_q/T,
  thus are energy-like observables.
  We expect that they are sensitive to the chiral phase transition to
  the extent that they receive non-analytic (singular) contributions in 
  addition to analytic (regular) 
  terms. 
  In the vicinity of the critical point, i.e. close to (T,H)=(T_c,0),
  non-analytic contributions are universal scaling functions of a scaling variable
  z=z_0 t  H^-1/βδ with t=(T-T_c)/T_c and z_0, T_c being 
  non-universal constants. Energy-like observables receive contributions from the
  scaling function, 
  f'_f(z)= df_f(z)/ dz, of the 3-d, Ø(N) universality class 
  <cit.>, which is the derivative of the scaling function f_f(z) that characterizes the singular part of the logarithm of the partition function. 
  For the heavy quark free energy we use the scaling ansatz,
  
      F_q(T,H)/T = A H^(1-α)/βδ f'_f(z) 
      		+f_ reg(T,H)   ,
  
  where A is another non-universal constant and
  the critical exponents are
  β, δ and α=2-β (1+δ).
  The regular contribution is an analytic function, which close 
  to (T_c,0) can be given as a Taylor series with even powers 
  in H,
  
      (T,H)= ∑_i,j  a^r_i,2j t^i H^2j≡∑_j p^r_2j(T) H^2j  .
  
  
  Using eqs. (<ref>) and (<ref>) the Polyakov loop expectation 
  value may be written as
  
      P_T,H = exp(-
      		A H^(1-α)/βδ f'_f(z) 
      		-f_ reg(T,H) )   .
  
  Note that for  (T,H) close to (T_c,0)
  eq. (<ref>) reduces to the usual 
  non-exponential scaling ansatz.
  
  
  
  Using the known asymptotic behavior of f_f(z) for 
  z→±∞ <cit.> 
  we obtain from eq. (<ref>) for the quark mass 
  dependence of  F_q/T  at fixed T and for 
  small H,
  
      F_q(T,H)/T∼
      a^-(T)+A p_s^-(T) H    , T< T_c 
      
      a^r_0,0+A a_1 H^(1-α)/βδ   , T=T_c 
      
      a^+(T)+p^+(T) H^2    , T> T_c
        ,
  
  with a^±(T)=A a_s^±(T)+(T,0) as well as p^+(T)=A p_s^+(T)+p^r_2(T) 
  receiving contributions from both the singular and regular terms. 
  For T≤ T_c the 
  dominant quark mass dependence arises from the singular 
  term only. In particular, we have
  
      a_s^±(T)    =    (2-α)  z_0^1-α c_0^± t |t|^-α  , 
      
      p_s^-(T)    =    (2-α-βδ) (-z_0 t)^1-α- βδ  ,  
      
      p_s^+(T)    =    (2-α- 2 βδ) c_1^+ (z_0 t)^1-α-2 βδ ,
  
  where we followed the notation of Ref. <cit.> with c_0^±,
  c_1^+
  and a_1 denoting coefficients appearing in the parametrization of the 
  scaling function f_f(z).
  
  For the  3-d, Ø(4) universality class these coefficients are given in 
  <cit.>. As we will present here only results from calculations for one non-zero value of the lattice spacing, a=1/8T, and do not perform a continuum extrapolation, we will use 
  
  the critical exponents and scaling functions of the 3-d, Ø(2) 
  universality class. We have extracted these from <cit.>, 
  where the scaling functions are parametrized in the Widom-Griffiths form. 
  We use a_1=0.4734, c_0^-=2.447,  c_0^+ = 2.728,  c_1^+=-0.678,
  along with β=0.349, δ=4.780 
  <cit.>.
  
  Making use of the relation between scaling functions, f_f(z) 
  and that of the order parameter, f_G(z) (see for instance 
  <cit.>), 
  we obtain in the vicinity of T_c,
  
      ∂ F_q(T,H)/T/∂ H = -A H^(β -1)/βδf'_G(z) 
                   +∂ f_ reg(T,H)/∂ H  .
  
  The derivative of F_q/T with respect to T is given by
  
      T_c∂ F_q(T,H)/T/∂ T =  
      A z_0 H^-α / βδ f”_f(z) + T_c ∂ f_ reg(T,H)/∂ T  .
  
  As α / βδ < 0 for the 3-d, Ø(2) as well
  as Ø(4) universality classes, the derivatives of
  P and F_q/T with respect to T do not diverge at T_c
  in the limit H→ 0. 
  On the other hand, the corresponding quark mass derivatives (eq. (<ref>))
  are expected to diverge when approaching the chiral limit at T_c 
  as well as when approaching T_c from below at H=0.
  A possible influence of the chiral phase transition on the
  quark mass dependence of P and F_q/T is thus much 
  easier to establish than on their temperature dependence. 
  
  
  
  
  
  
  Computational setup and data analysis.—
  In this study we analyze properties of (2+1)-flavor QCD where 
  the strange quark mass has been kept fixed to its physical value
  and the two degenerate light quark masses are varied in the range 
  H=1/20 to 1/160, corresponding to a pion mass
  160 MeV≳ m_π≳ 58 MeV <cit.>.
  The analysis is performed on sets of gauge field configurations that 
  had been generated
  previously by the HotQCD collaboration  
  <cit.> on 
  lattices of size 32^3× 8 (H=1/20, 1/27), 40^3× 8 
  (H=1/40) and 56^3× 8 (H=1/80, 1/160).
  
  using highly improved staggered quarks (HISQ) <cit.>
  and the tree-level improved Symanzik gauge action.
  These data have previously been
  used in calculations of pseudo-critical temperatures <cit.> and the chiral phase transition temperature of (2+1)-flavor QCD <cit.>.
  We also have generated additional configurations for H=1/40, 1/80. 
  For the scale setting we use the kaon decay constant 
  obtained in calculations with the HISQ action, i.e.
  f_K=156.1/√(2) MeV <cit.>.
  
  Observables have been calculated on lattices
  with temporal extent N_τ=8 and N_σ/N_τ = 4-7.
  
  For  H =1/80 we have results for the entire range of aspect ratios, N_σ/N_τ. 
  An analysis of finite volume effects on P shows no
  significant volume dependence of  P in the entire
  temperature range. Results for different N_σ/N_τ
  agree within errors
  of about 1%. We thus can safely neglect any finite volume
  corrections to our results for P. 
  The statistics used in this study can be found in <cit.>,
  along with more details for the finite volume dependence of P.
  
  
  
  
  
  Results.—
  In Fig. <ref> (left) we show results for the
  derivative of F_q(T,H)/T with respect to H as function of T.
  It is apparent  that  ∂ (F_q/T)/∂ H
  decreases with quark mass at high temperature and increases at low
  temperature. This is consistent with an approach towards zero at high 
  T and a non-vanishing, strongly temperature dependent constant
  at low T. Such a pattern is in accordance with the expected quadratic
  dependence on H for T>T_c and the linear dependence for T<T_c
  given in eq. (<ref>).
  Although errors are large for the results obtained with
  the smallest light quark mass, H=1/160, it is evident that 
  ∂ (F_q/T)/∂ H has maxima at
  T∼ 145-150 MeV which are close to the chiral phase transition
  temperature on lattices with temporal extent N_τ=8 determined
  in <cit.>, T_c^N_τ=8=144(2) MeV. With decreasing H they
  approach T_c^N_τ=8 and the peak height increases.
  
  In Fig. <ref> (middle) we show ∂ (F_q/T)/∂ H
  re­scaled with the appropriate power of H expected
  from the Ø(2) scaling ansatz and plotted versus z. 
  The good scaling behavior
  suggests that ∂ (F_q/T)/∂ H will indeed diverge in the chiral limit and that H-dependent
  contributions to F_q/T, arising from regular terms, are small compared to
  those coming from the singular part. This motivated a fit ansatz for the
  regular term that is independent of H, i.e. we use (T,H=0) in all our fits. 
  
  A fit to ∂ (F_q/T)/∂ H thus only involves
  the singular term of eq. (<ref>). We performed this
  3-parameter fit for all data sets with H≤ by either including or leaving out the data for
  H=1/27. These fits are shown in
  Fig. <ref> (left), and the corresponding fit parameters,
  A, T_c, z_0, are given in Table <ref>.
  In particular, we find that T_c, obtained from these fits, as well as 
  the non-universal scale parameter z_0, agree
  well with earlier fit results for chiral susceptibilities
  <cit.> in (2+1)-flavor QCD.
  The rescaled fits are 
  also shown in Fig. <ref> (middle).
  
  
  
  
  
  
  Data for F_q(T,H)/T are shown in Fig. <ref> (right).
  They have been fitted to the scaling
  ansatz eq. (<ref>) using only the 
  constant and linear H-independent
  terms in the regular part as fit parameters and keeping fixed
  
  the three non-universal constants, determined in the 
  previous step, in the singular part.
  
  
  
  
  The T-range and data 
  included in the fit are shown in the inset. The  resulting fit 
  parameters, a^r_0,0, a^r_1,0,  are 
  given in Table <ref>. In this figure we also show the 
  static-light meson contribution to F_q/T calculated in the hadron-gas approximation <cit.>. In chiral perturbation theory this also 
  gives a linear dependence on H at low temperature <cit.>.
  
  
  
  Once we have determined all five fit parameters for F_q(T,H)/T, we can plug them
  into eq. (<ref>) to arrive at a parameter-free description of the T
  and H dependence of P. The thus determined curves are shown in
  Fig. <ref>. As seen in the inset, they agree well with P
  data near T_c^N_τ=8, which suggests the behavior of P is explained
  well by chiral scaling behavior in this region and serves
  as a consistency check of our approach.
  
  In Figs. <ref> (right) and <ref> (top) we also 
  show the chiral limit (H=0) results for F_q/T and P. For the
  former this is a sum of regular and singular contributions,
  
      F_q(T,0)/T = a^r_0,0 + t ( a^r_1,0 + A^± |t|^-α)   ,
  
  with A^±=(2-α) z_0^1-α c_0^± A (see
  eq. (<ref>)).
  Although at T_c the contribution to the slope is entirely given by
  the regular term a_1,0^r, close to T_c this
  contribution gets to a large extent canceled by
  the singular contributions, A^± |t|^-α. 
  This is the origin of the well known spike in
  specific-heat like observables (2^ nd
  derivatives with respect to T) in the  Ø(N)
  universality classes. In the chiral limit our fit
  results suggest the appearance of such a spike in the
  temperature derivatives of  F_q(T,0)/T as well as
  P. For the former we obtain
  from eq. (<ref>), 
  
      T_c ∂ (F_q(T,0)/T)/∂ T =  a^r_1,0( 1 + R^±  
             |t|^-α)    ,
  
  with R^±=(1-α)A^±/a^r_1,0. R^+ is given in Table <ref> and 
  R^+/R^-= c^+_0/c^-_0=1.12(5) is a universal ratio <cit.>. 
  This makes it evident that already for |t|=0.01 the slope of -F_q/T is about a factor 5 smaller than at T_c . 
  
  The basic features found in our analysis of 
  (2+1)-flavor QCD are quite similar to those found in 
  the analysis of 3-d, Ø(2) symmetric spin models 
  <cit.>. Also in that case a large 
  cancellation of contributions arising from regular
  and singular terms is found; the spike in the specific 
  heat, C_V, is concentrated in a temperature interval 
  of about 1% around T_c, and C_V changes by almost 
  a factor 10 in this temperature interval. In 
  Fig. <ref> we show a comparison of the
  scaling behavior of -T_c ∂ (F_q(T,0)/T)/∂ T, 
  obtained for the case of (2+1)-flavor QCD, and C_V 
  in Ø(2) spin models at H=0 <cit.>.
  The inset shows the development of this sharp peak as the quark mass 
  decreases. It makes clear that this feature becomes visible only for
  H being substantially smaller than the region H ≃ 10^-2 that is
  accessible in current lattice QCD calculations. 
  
  
  
  
  Conclusions.–
  We have examined the quark mass dependence of the Polyakov loop 
  expectation value and the heavy quark free energy extracted from it.
  We provided evidence for the influence of 
  chiral symmetry restoration that manifests in the singular behavior
  of the quark mass derivatives of P and F_q/T and arises from the 
  energy-like 
  behavior of these observables with respect to chiral transformations. 
  These derivatives  diverge at T_c in the chiral limit, consistent with
  the expected behavior for energy-like observables in the 3-d, Ø(2) 
  universality class. 
  
  We showed that at finite values of the lattice spacing
  the relative distribution between singular and regular contributions to 
  the energy-like variables in QCD seems to be
  similar to that in the 3-d, Ø(2) spin model. In particular, 
  in the chiral limit the very narrow spike showing up in specific-heat like 
  observables, which results
  from a partial cancellation of singular and regular contributions combined
  with the existence of a quite small, negative critical exponent α,
  is expected to show up also in the T-derivatives of F_q/T and P.
  
  As the critical exponent α also is negative and small in the 
  Ø(4) universality class, similar behavior of F_q/T and P
  is expected to persist in the continuum limit.
   However, as |α| is an order of magnitude larger 
  compared to the Ø(2) case, the spike may become broader
  and may also be more prominent already for larger quark masses.
  This may be tested also in effective model calculations.
  
  All data from our calculations, presented in the figures of this paper, can
  be found in <cit.>.
  
  
  
  
  
  
  § ACKNOWLEDGMENTS
  
  This work was supported by the Deutsche Forschungsgemeinschaft
  (DFG, German Research Foundation) Proj. No. 315477589-TRR 211; and by 
  the German Bundesministerium für Bildung und Forschung through
  Grant No. 05P18PBCA1.
  Numerical calculations have been made possible through PRACE grants at CSCS, 
  Switzerland, and grants at the Gauss Centre for Supercomputing and NIC-Jülich, Germany. These grants provided access to resources on Piz Daint at CSCS as well as on JUQUEEN and JUWELS at NIC. This research also used awards
  of computer time at NERSC, Berkeley Laboratory, provided by the DOE ASCR 
  Leadership Computing Challenge (ALCC) program. 
  We thank the HotQCD Collaboration for providing access to their latest data
  sets and for many fruitful discussions.
  
  
  
  
  
  
  