Polyakov loop fluctuations in the presence of external fields

We study the implications of the spontaneous and explicit Z(3) center symmetry breaking for the Polyakov loop susceptibilities. To this end, ratios of the susceptibilities of the real and imaginary parts, as well as of the modulus of the Polyakov loop are computed within an effective model using a color group integration scheme. We show that the essential features of the lattice QCD results of these ratios can be successfully captured by the effective approach. Furthermore we discuss a novel scaling relation in one of these ratios involving the explicit breaking field, volume, and temperature.

We study the implications of the spontaneous and explicit Z(3) center symmetry breaking for the Polyakov loop susceptibilities. To this end, ratios of the susceptibilities of the real and imaginary parts, as well as of the modulus of the Polyakov loop are computed within an effective model using a color group integration scheme. We show that the essential features of the lattice QCD results of these ratios can be successfully captured by the effective approach. Furthermore we discuss a novel scaling relation in one of these ratios involving the explicit breaking field, volume, and temperature.

I. INTRODUCTION
Understanding deconfinement and chiral symmetry restoration, as well as exploring their far-reaching consequences [1][2][3] remain challenging in the study of heavy ion collisions. A robust description of the phenomenon is necessary for reliably analyzing many observables of experimental interests, such as fluctuation observables, transport coefficients or the production rate of photon and dilepton [1,4,5].
In the limit of infinitely heavy quarks, deconfinement can be identified with the spontaneous breaking of the Z(3) center symmetry [27,28]. The Polyakov loop [29][30][31][32] serves as an order parameter for the phase transition. In some effective models [16][17][18][19][20] a potential is constructed to describe its behavior. It is possible to constrain the parameters of the potential using the LQCD results on the thermodynamic pressure [33,34] and the Polyakov loop [35] in a pure gauge theory. In particular, the latter dictates the locations of the minimum of the potential at different temperatures.
There is another class of independent observables which are sensitive to the Z(3) center symmetry -the susceptibilities of the real and imaginary parts, as well as of the modulus of the Polyakov loop [36]. These quantities measure the fluctuations of the order parameter field. To describe them in an effective model, not only the location, but also the curvatures around the minimum of the Polyakov loop potential have to be adjusted [37].
In a pure gauge theory, ratios of these susceptibilities have been demonstrated [36] to be excellent probes of deconfinement. They display a θ-function like behavior across the transition temperature T d , with well defined low temperature limits deducible from general theoretical constraints and the Z(3) symmetry.
In a recent study [38] these ratios have been computed in numerical simulations of LQCD with 2+1 light flavors. Unfortunately, the task of extracting useful information from these quantities is more involved than originally thought. Most importantly, many pertinent fea-tures of the ratios are smoothed out in the presence of dynamical fermions, as well as after prescribing a renormalization. In addition, the results are still marred by issues of renormalization scheme dependence and it is far from clear how to connect them to calculations made in an effective model.
Despite these difficulties, we stress that there are strong theoretical motivations for studying and understanding these ratios. For one thing, the widely used order parameter, i.e. the renormalized Polyakov loop computed by LQCD, is a renormalization scheme dependent quantity [35,39]. This calls into question the physical relevance of the deconfinement features deduced from it, for example, the transition temperature T d [10,11] extracted from its inflection point. It is therefore crucial to study the deconfinement phenomenon from the perspective of these additional observables, and investigate whether a coherent picture can be obtained. They can also be used to signal the strength of the explicit symmetry breaking field.
In this paper we compute the susceptibility ratios within an effective model. This allows a transparent study of how aspects of center symmetry breaking, explicit and spontaneous, manifests in the ratio observables. The approach also provides some simple explanations to many features of the LQCD results.
The article is organized as follows: In Sec. II we review the derivation of the Polyakov loop susceptibilities using the color group integration approach. The method is illustrated by computing one of the ratios, R A , in the presence of explicit symmetry breaking field for a Gaussian model. In Sec. III we present the effective Polyakov loop potential for this study and analyze the explicit breaking field and volume dependence of the model susceptibility ratios. A novel scaling relation for R A is also presented. In Sec. IV, we compare the model results with LQCD calculations. In Sec. V we present the conclusion.

A. formalism
In this work we compute the various susceptibility observables using a color group integration scheme [37, 40, arXiv:1801.08040v1 [hep-ph] 24 Jan 2018 41]. The partition function in this approach is expressed as where (x, y) stands for the real and imaginary 1 part of the Polyakov loop, promoted to a (homogeneous) classical field degree of freedom. The actual Polyakov loop potential (U [x, y]) of choice will be presented in Sec. III. Expectation value of an arbitrary operator within this approach is computed via Thus, e.g., the expectation value of the Polyakov loop can be readily obtained from Staying within the real sector and considering explicit symmetry breaking along the real axis imply y = 0. However, fluctuations of the order parameter can be explored along the longitudinal (real) and transverse (imaginary) directions, as well as that of its absolute value: From these, two independent ratios are derived: Note that χ A = χ L + χ T and R A = 1 + R T . In the large volume limit, it can be shown that the two susceptibilities in Eqs. (4)-(5) approach the meanfield results where C is the correlation matrix [37,42], defined as 1 In this study we shall stay exclusively in the real sector and there is no ambiguity in identifying the longitudinal and transverse directions with the real and imaginary axes, respectively.
This gives a transparent interpretation of the susceptibilities as the inverse of curvatures of the effective Polyakov loop potential. Note that all of these quantities are to be evaluated at (x, y) → (x 0 , y 0 ), determined from the gap equations On the other hand, χ A in Eq. (6) does not have a valid mean-field limit. Nevertheless, it can be readily computed in the current color group integration scheme.

B. Gaussian model with an explicit symmetry breaking field
The Gaussian model has proved useful for understanding the low temperature behavior of the susceptibility ratios in a pure gauge system. Inserting a potential of the form in Eq. (1), the following non-trivial manifestations of the Gaussian limit can be derived [36]: These low temperature relations have been verified by lattice calculations in a pure gauge theory.
To investigate how these ratios behave in QCD with dynamical quarks, we extend the discussion to include finite explicit symmetry breaking. To this end, we perform a substitution It is straightforward to derive an exact expression for R A . 2 . The result reads, with where I n (x) is the modified Bessel function of the first kind of the n-th order, and   . We adopt the following scheme in presenting our results: different colors (and line types) correspond to different explicit breaking strengths; while different symbols denote different volumes. In this numerical study, we fix α = 1, h0 = 1, and V0 = (6.9 fm) 3 .
Especially we extract the following important limits: This generalizes Eq. (13) to the case of a finite explicit symmetry breaking.
To gain some familiarity with the ratio R A , we examine the quantity as a function of temperature: (1) at fixed volume but for different h, and (2) at fixed h but for different values of the volume. In this numerical study, we fix α = 1, h 0 = 1, and V 0 = (6.9 fm) 3 . The results are shown in Fig. 1. As expected, the R A ratio interpolates between the two known limits [36]: from the Z(3)-symmetric phase (R A = 2−π/2) to the Z(3)-broken phase (R A = 1). At a fixed volume, increasing the breaking strength h makes the ratio approach unity at lower temperature. The quantity also exhibits a strong volume dependence, as seen in Fig. 1 right.
All the results presented in Fig. 1 originate from the single expression Eq. (15). In the Gaussian model, the breaking field h always enters via a combination of volume V and the parameter α dictated by Eq. (17). This leads to a difficult situation that as V → ∞, R A → 1, regardless of the value of h and temperatures. To obtain useful information from this quantity, it is necessary to work in a finite volume setting. Alternatively, one can study R A as function of the scaling variable ξ. We shall revisit some of these issues for the full effective model in Sec. III.

A. effective Polyakov loop potential
The Gaussian model discussed above, though generalized to include an explicit symmetry breaking field, does not describe the spontaneous Z(3) symmetry breaking. 3 To examine the susceptibility ratios in a setting that is relevant to QCD, a Polyakov loop potential [37], capable of handling the latter aspect will be employed: Here M H (l,l) is the SU (3) Haar measure The temperature dependent model parameters (A, B, C, D) are given in Ref. [37] and will not be repeated here. Note that to implement the color group integration scheme in Eq. (1), we use  . We adopt the same presentation scheme as in Fig. 1 for the results. The volume independence of RT is evident: results with different symbols (but same color) all fall on the same line. In this study, h0 is given in Eq. (22), and V0 = (6.9 fm) 3 . The "PNJL" line denotes the result of a PNJL model for 2+1 light flavors described in Sec. IV. Also shown are the LQCD results at flow time f = f0 extracted from Ref. [38].
The potential U G in Eq. (19) is particularly suited for the current study. Most importantly, the known susceptibilities at zero explicit breaking are reproduced by construction. This is not the case for other commonly used Polyakov loop potentials [16][17][18]. For example, the polynomial potential introduced in Ref. [16] leads to the result R T > 1 for T > T c , which is another manifestation of the "negative susceptibility" problem discussed in Ref. [14,22]. Imposing the Haar measure to the potential [17,18] effectively restricts the Polyakov loop to the target region and thus improves the theoretical description. In fact, the present model builds on this observation and further constrains the curvatures of the potential using the available LQCD results [36] on the susceptibilities in a pure gauge theory.
Furthermore, we consider a linear explicit breaking term h × x. The functional form of the breaking field h is known within the approximation of a one-loop expansion of the fermionic determinant [42,43]. As we aim at understanding the ratios on the qualitative level, we shall employ the following basic form for h = c × h 0 , with where For quark masses we choose m l = 5 MeV for up and down quarks and m s = 100 MeV for the strange quarks. 4 Eq. (23) improves on the previous study [42] by taking the quantum statistics of quarks into account. In practice this leads to an increase of ≈ 10% in the total strength h 0 over the one imposing a further Boltzmann approximation. We shall also allow for an arbitrary prefactor c to manipulate the strength of h. In fact, an attempt will be made to infer its magnitude from the LQCD results on the ratio observables.
In a previous study [42] we have calculated the critical strength of the breaking field, h crit. , for the phase transition to turn from the first order to the second order, that is, the critical end point (CEP) for the Z(3) transition: The breaking field h 0 in Eq. (22) exceeds this limit for all temperatures of interest, meaning that a crossover transition is expected. This is evident in the ratio observables computed in the full model, as shown in Fig. 2. Starting with the ratio R A , shown in Fig. 2 left, we first notice the similarity between the full model results and those from the Gaussian model. Indeed, the R A ratio interpolates between the two known theoretical limits: (2 − π/2 ≈ 0.43) and 1. The expected behaviors from varying the breaking strength and the volume are also verified.
Turning now to the ratio R T , shown in Fig. 2 right, the immediate observation is the volume independence  of the quantity. This is evident from the fact that results with different symbols (but same color) all fall on the same line. This suggests that the finite volume V 0 = (6.9 fm) 3 we selected is sufficiently large. Indeed, we have checked that the ratio R T approaches the mean field value, dependent only on the intensive variables T and h. Increasing the breaking strength h makes the ratio deviates from the known Z(3)-symmetric limit (R T = 1) at lower temperature. On the other hand, the value of R T at large temperatures is not dictated by the Z(3) symmetry. Instead, it can be related to the color screening properties of the QCD medium.
Before proceeding to compare the effective model calculations with LQCD (see Sec. IV), we first discuss an interesting observation of R A ratio, namely, a scaling relation inspired by the Gaussian model.

B. scaling relation of RA
The Polyakov loop potential U G employed in Eq. (19) is clearly non-Gaussian. Nevertheless, we can consider a generalized double-Gaussian approximation to the potential: where the model parameters α 1 , α 2 , andh are constructed to match with These coefficients can be readily obtained by a mean-field calculation of the original potential at vanishing breaking field. It is clear that the approximation scheme operates by constructing local double-Gaussian potential along the line of minima of U G . The advantage of performing such an expansion is that it allows a direct computation of the ratio R A with an equation analogous to the single Gaussian limit studied previously in Eq. (15). The generalized equation reads where F can be obtained with an integral involving the modified Bessel function. (Details in the appendix A) According to this equation, the functional dependence in (T, h, V ) of R A can be uniquely determined by the scaling variable ξ(T, h, V ) and R T (T, h), via This translates to the following: Provided that the generalized Gaussian approximation is valid, all the data point of R A (T, h, V ) will collapse on a single universal line when plotted against ξ, with the choice a "physical" R T = R T (T, h). A direct numerical computation confirms that it is indeed the case, and the result is shown in Fig. 3 left. While the observation is theoretically interesting, it also indicates a rather limited information contained in this observable. For example, the key information about the magnitude of the explicit breaking field can as well be extracted from R T . Nevertheless, Eq. (28) may serve as a useful diagnostic for analyzing R A .

IV. COMPARISON WITH LQCD RESULTS AND A PNJL MODEL
As discussed in the introduction, the renormalized Polyakov loop computed by LQCD is a renormalization scheme dependent quantity [38,39]. It obscures the physical relevance of the derived deconfinement features, e.g. the T d extracted from the inflection point, and complicates the comparison of LQCD results with those obtained in an effective approach.
One of the original motivations for introducing the susceptibility ratios as probe of deconfinement is the removal of both the cutoff and the scheme dependence.
The assumption is that if the Polyakov loop susceptibilities are renormalized the same way as the Polyakov loop, the multiplicative renormalization factor will be canceled against each other.
Contrary to this expectation, recent studies [38,39] report a substantial cutoff dependence in these ratios in QCD with 2+1 light flavors. This is evident from the N τdependence observed in the "bare data" of R A and R T in Ref. [38]. This seems to suggest that renormalizing the Polyakov loop alone does not guarantee the renormalization of the susceptibilities.
In the effective model, the behavior of R T is largely determined by the explicit breaking field h. It is possible, therefore, that the cutoff dependence observed in R T can be associated with the cutoff dependence of h. In fact, it is non-trivial to obtain a continuum extrapolation of the explicit breaking strength h from LQCD that is suitable for comparison with effective model [42,[44][45][46][47].
Furthermore, if we mimic the N τ -dependence of the ratio observables as a change of prefactor in h in the effective model, we can reproduce the same trend in the ordering of curves of R T and R A , namely, from top to bottom for R A in increasing N τ and the reverse order for R T . (See Fig. 17 and 18 of Ref. [38].) This suggests that the two sets of "bare" data are connected, and the connection may be due to h. Using the gradient flow method [38,39], it is possible to renormalize the susceptibilities and the ratio observables. We have selected LQCD results with the "f = f 0 " flow time to compare with our effective model calculations. They are presented in Fig. 2. We note that a reasonable agreement can be attained if we choose an explicit breaking field of strength ≈ (1 − 2) × h 0 and a physical volume of V 0 ≈ (6.9 fm) 3 .
It is straightforward to extend this study to incorporate effects from the spontaneous chiral symmetry breaking. For this purpose, we study a (2+1)-PNJL model, combining the NJL model in Refs. [48,49] and the Polyakov loop potential U G in Eq. (19). The computation is similar to the one presented in Sec. II, and we simply show the major results in Fig. 2 and 3. These results can be readily understood by studying an effective Z(3) breaking strength for the PNJL model. It can be computed via Eq. (22), except for using the constituent quark masses in lieu of the current ones. The behavior of such an effective breaking strength compared to h 0 , i.e. the prefactor c PNJL (T ) for the PNJL model, is shown in Fig. 4. It is evident that the PNJL model leads to a smaller explicit Z(3) breaking, and the observables R T and R A behave accordingly. However, the LQCD results of ratio observables seem to indicate a stronger breaking strength. This suggests that a naive implementation of the coupling between quarks and the Polyakov loop may be inadequate and a more sophisticated treatment including the back-reaction of dynamical quarks on the gauge sector, which could modify the Polyakov loop po-tential U G , may be necessary.
Furthermore, there is still substantial "flow time" dependence in these observables, reflecting a further renormalization prescription dependence. 6 In particular, the large flow time (f = 3f 0 ) result of R T shows a relatively low value (≈ 0.7) at low temperatures, instead of the expected Z(3)-symmetric limit of unity. It is still possible to describe such R T within our effective model, though it requires a rather large h. This also naturally explains the observation that R A → 1 at large flow time. We however find this situation unsatisfactory, since it is more natural to expect the physical breaking field h to be free of the renormalization scheme dependence. And it is this quantity that we hope to extract from the LQCD.
For reference we also compute the Polyakov loop in the effective model and compare with the LQCD result. This is shown in Fig. 3 right. Unlike the case of the ratio observables, we see that the effective model essentially fails to describe the LQCD result. Similar discrepancy has been reported by the matrix model [13]. This may again be due to the scheme dependence and we shall explore this topic in more detail in a future publication.

V. CONCLUSION
This study has demonstrated how features of deconfinement emerge in the ratios of Polyakov loop susceptibilities within an effective model. For the ratio R A , we find a characteristic volume dependence along with temperature and the explicit Z(3) symmetry breaking strength h. In a Gaussian approximation scheme, all these can be subsumed into a single scaling equation. For the ratio R T , we find a minimal volume dependence, which makes it a robust probe of the strength of the explicit breaking term.
On a qualitative level, the effective model is capable of describing many features of the LQCD results. These include the low and high temperature limits of the ratios, and the connection between R A and R T , possibly via h. Nevertheless, it is important to bear in mind that the LQCD results, which we compare our model results to and estimate the effective strength of h from, still suffer from a renormalization scheme dependence. This is an urgent issue to be tackled to achieve a meaningful comparison of effective approaches with LQCD, and will be pursued in future works.
It would also be interesting to investigate how the susceptibility observables behave at large baryon densities [50] and react to other external fields, e.g. a magnetic field [51,52]. Since the curvatures dictate how reluctant the system is to deviate from the equilibrium position in the presence of external disturbances, we expect these 2 ⇡/2 with where K n is the modified Bessel function of the second kind of the n-th order.
It is instructive to study the limits of ξ 1 and ξ 1 for a general R T : where Here K, E are the complete elliptic integral of the 1st and 2nd kind respectively, defined as For R T = 0, we obtain the SU(2) [36,53] limit: A schematic plot of R A (ξ, R T ) in Eq. (A6), as a function of ξ for a given R T (T ), is illustrated in Fig. 5 left, together with the functional dependence on R T of various expansion coefficients in Fig. 5 right.