Divergences in the quark number susceptibility : The origin and a cure

Quark number susceptibility on the lattice, obtained by merely adding a $\mu N$ term with $\mu$ as the chemical potential and $N$ as the conserved quark number, has a quadratic divergence in the cut-off $a$. We show that such a divergence already exist for free fermions with a cut-off regulator. While one can eliminate it in the free lattice theory by suitably modifying the action, as is popularly done, it can simply be subtracted off as well. Computations of higher order susceptibilities, needed for estimating the location of the QCD critical point, then need a lot fewer number of quark propagators at any order. We show that this method of divergence removal works in the interacting theory.


I. INTRODUCTION
The phase diagram of the strongly interacting matter described by Quantum Chromodynamics(QCD) has been a subject of intense research in the recent years. Usual weak coupling perturbative approach may work for sufficiently high temperatures. However, the gauge interactions are likely to be strong enough for temperatures close to Λ QCD , the typical scale of QCD, necessitating strong coupling techniques. Lattice QCD is the most successful non-perturbative technique which has provided us with some interesting results pertaining to the phase diagram. It is now fairly well known from independent lattice studies that the transition from the hadron to the quark gluon plasma phase at zero baryon density is a crossover [1][2][3]. At non-zero density, or equivalently nonzero quark chemical potential µ, one has to face a sign problem : quark determinant is complex. This does not allow for an importance sampling based Monte Carlo study. Several ways have been advocated in the recent years to circumvent the sign problem in QCD [4][5][6][7]. For two light quark flavours, a critical end-point is expected from the perturbative studies of model quantum field theories with the same symmetries as QCD [8,9]. If present, the critical-end point would result in the divergence of the baryon number susceptibility. Thus its Taylor expansion [7] at finite baryon density as a series in µ B /T can be used to compute the radius of convergence, and therefore, an estimate of the location of the critical end-point [10,11]. First such estimates of the radius of convergence of the Taylor series have predicted the critical end-point to be at T E /T c = 0.94 and µ B /T E = 1. 8(1) [11]. Recently, a study on a finer lattice has suggested the continuum limit to be around T E /T c = 0.94(1) , µ B /T E = 1.68 (5) [12]. In the heavy-ion experiments, the fluctuations of the net proton number could act as a proxy for the net baryon number. The STAR experiment at Brookhaven National Laboratory has reported the measurements for the fluctuations of the net proton number for a wide range of center of mass energy √ s, of the colliding heavy ions between 7.7 and 200 GeV. At √ s = 19.6 GeV the experimental data is observed [13] to deviate from the predictions of the proton fluctuations for models which do not have a critical end-point, and is similar to the lattice QCD-based predictions [14] for a critical point, signaling its possible presence. It would be thus important to have a thorough understanding of the systematics of the lattice QCD results and make them as much reliable as possible.
In addition to the usual suspects, such as continuum extrapolation or effects due to the finiteness of the lattice spacing, the scale-setting, and the statistical precision of the measurements, a key new important factor is that the radius of convergence estimate requires ratios of as many higher orders of quark number susceptibilities(QNSs) as possible. Currently the state of the art is the eighth order QNS [10,11]. It is very important to verify whether the existing results are stable if ratios of further higher order of QNS are taken into account. In order to calculate the QNS of order m, one has to take the m th -derivative of the free energy with respect to the quark chemical potential. Since the popular method of incorporating the chemical potential on the lattice is through exp(±µa) factors multiplying the forward and the backward temporal gauge links respectively of the fermion operator [15,16], there is an ever increasing proliferation of terms of varying sign as m increases. Their large number and the large cancellations amongst them at a specific order make it difficult to increase m beyond eight at present. Introducing the chemical potential by a µN -term, where N is the corresponding conserved charge, leads to both much fewer terms and lesser cancellations at the same m [17], thereby reducing the computational cost up to 60 % for terms up to the eighth order; more savings ought to accrue by going to even higher orders. More precise Taylor coefficients and more terms in the Taylor expansion can potentially lead to a better control of the QCD equation of state at finite baryon density which will be needed for the analysis of the heavy-ion data from the beam energy scan at RHIC as well as the future experiments at FAIR and NICA.
In this paper, we discuss whether such a linear in µ approach is viable or has unsurmountable problems by comparing with the usual exponential in µ method. In section II, we revisit the number density for non-interacting fermions in the continuum using a cut-off regulator. We point out that divergences appear already in the continuum free theory when the cut-off regulator is taken to infinity contrary to the conventional wisdom. We then discuss the lattice approaches to tackle to this divergence and argue why they should work even after the introduction of the gauge interactions. We verify this in section III by performing continuum extrapolation of the second and fourth order QNS for quenched QCD for the linear approach. This is the most important result of our paper. We discuss its possible consequences and the extensions to higher order QNS.

II. THERMODYNAMICS OF NON-INTERACTING FERMIONS
QCD thermodynamics can be derived from its partition function, written in the path integral formalism [18] as, where ψ,ψ and A µ represent the quark, anti-quark and the gluon fields respectively, whose the color indices are not written explicitly above. µ is the chemical potential for the the net quark number with the corresponding conserved charge being d 3 xψγ 4 ψ. Generalizations to various conserved flavour numbers is straightforward. For simplicity, we will consider only a single flavour with the baryonic chemical potential µ B = 3µ q . Appropriate derivatives of Z lead to various thermodynamical quantities, e.g., the quark number density, or equivalently (1/3) the baryon number density, is defined as, Earlier attempts to discretize the above theory to investigate the finite baryon density physics on a space-time lattice revealed µ-dependent quadratic divergences in the number density and the energy density when the chemical potential is introduced in the quark Dirac operator by multiplying it with the corresponding conserved charge on the lattice. These divergences, which appear as a µ/a 2 term in the expression for the lattice number density with a as the lattice spacing, are present even if the gauge interactions are absent. Through explicit calculation of the number density for non-interacting fermions on the lattice, it was then shown [15,16,19] that suitable modification of the µN term in the action, eliminates these divergent terms on the lattice, and yields a finite a → 0 continuum limit. Numerical studies of the QNS for the interacting theory subsequently confirmed that once the free theory divergences are thus eliminated, no further divergences arise [20,21]. A succinct way to describe all the various actions is to introduce functions f (µa)[g(µa)] as the multiplying factors for the forward (backward) timelike gauge fields on the lattice. While for the naive discretization, f = 1 + aµ and g = 1 − aµ leads to a divergent baryonic susceptibility in the continuum limit, the choice f = exp(aµ) and g = exp(−aµ) does not.
It is now commonly believed that such an exponentiation is necessary on the lattice as the divergence is a lattice artifact. Clearly since all derivatives of f and g are nonzero for the exponential case, whereas only the first derivative is nonzero for the linear case, higher order QNS are a lot simpler for the latter. Furthermore, for fermions with better chiral properties such as the Overlap fermions or the Domain Wall fermions, the exponential form leads to a loss of the chiral symmetry [22] for nonzero µ. Indeed the only chiral symmetry preserving form these fermions have for finite µ and a is the linear form [23] for f and g. This motivates us to revisit the issue of the nature and origin of these divergences when the chemical potential enters linearly instead of the exponential form. As we show below, the divergences are present in the continuum as well, and the lattice regulator simply faithfully reproduces them. While one can employ the freedom of lattice action to eliminate them, it is not necessary. Indeed, one can perhaps employ simpler subtraction methods to eliminate them, as we demonstrate in this paper.

A. Continuum free fermions
Results for the continuum free fermions are easily found in textbooks [18]. We review them below solely with the idea of pointing out explicitly the µ-dependent divergences present in them. For simplicity, we consider only massless fermions though this derivation can be easily extended for finite mass. The expression for the number density for free fermions is easily obtained from the definitions above as where p 2 = p 2 1 + p 2 2 + p 2 3 and ω j = (2j + 1)πT . Here we choose the gamma matrices to be all Hermitian as is common in lattice studies. The continuum convention followed in the standard texts has only γ 4 as Hermitian and the other gamma matrices are anti-Hermitian. The expression in Eq. (3) can be evaluated by the usual trick of converting the sum over energy states to a contour integral. We compute this quantity at T = 0 only. At zero temperature the energy eigenvalues become continuous and lie in the range [−∞, ∞]. The expression for the number density then becomes Under a variable transformation ω + iµ = ω ′ , it can be recast as In order to compute the integral carefully, we impose a cut-off Λ on all the four momenta and use the contour in the complex ω plane in Fig. 1, leading to It may be remarked that for other choice of the Dirac matrices the contour in Fig. 1 is rotated in the complex plane and the poles appear on the imaginary axis at the Matsubara frequencies. While the first term arises from the residue of the pole of the integrand enclosed by the contour, the last three terms in the Eq.(6) arise from closing the contour. The line integral 1 is identically zero because the integrand is an odd function. The sum of the line integrals 2 and 4 is Since Λ >> µ, Λ being the cut-off, expanding the logarithm and retaining the leading term, It is straightforward to do the remaining two momenta integrals in Eq. (6). The first term leads to the usual µ 3 term. However, the sum of the two line integrals 2 and 4 in Eq. (8) yields a µΛ 2 divergence in the expression for number density as below, Note that the leading diverging Λ 3 -contribution is the same for the line integrals 2 and 4 and it does cancel. It is the non-leading µ-dependent term which leads to the Λ 2 divergence above and is usually missed in the textbooks [18]. A detailed comparison of our calculation with that outlined in Ref. [18] is provided in Appendix A. The same exercise, as above, can be repeated for the free fermions on the lattice [15,16,19]. For the sake of brevity, we just emphasize the differences with respect to the continuum case above. For the linear µ case, corresponding the naıve discretization, the number density is, where p 2 = (ma) 2 + sin 2 (ap 1 ) + sin 2 (ap 2 ) + sin 2 (ap 3 ). In the exponential µ case, the term in the brackets above changes to (sin ω n cosh µa + i sinh µa cos ω n ). Following Ref. [24], and employing the earlier introduced functions f (aµ) and g(aµ), both the linear µ and the exponential µ cases can be covered by writing the expression in the bracket as R sin(ω n + iθ) with R 2 = f · g and tanh θ = (f − g)/(f + g). The usual periodicity in the lattice energy is then restored for µ = 0 as well. The same contour as in Fig. 1, with µ → θ and ±Λ → ±π with the appropriate dimensional factors, can then be used to compute the free fermion density. Due to the periodicity of the sine and the cosine functions in the integrand, the line integrals 2 and 4 cancel each other for all f and g. For the linear µ case, it is the line integral 3 that gives rise to the µ/a 2 or µΛ 2 divergence, while R = 1 for the exponential µ form ensures that such a term is absent for it. As pointed out in Ref. [24], any choice of f and g satisfying R = 1 ensures the absence of the divergences.

III. QUENCHED RESULTS ON THE LATTICE
The analytical proof of Refs. [15,16,24] for the lack of divergences in the quark number susceptibility in the exponential µ case , outlined briefly above in the previous section, was for non-interacting fermions. No equivalent proof exists for the interacting fermions even for the exponential case. On the other hand, it is easy to check that for chiral fermions, such as the staggered quarks, the chiral symmetry is maintained for µ = 0. For the linear µ case, one even deals with conserved baryonic currents on the lattice, aµ being the coefficient of the conserved baryon number on the lattice. One therefore expects no further divergences to arise, and no extra renormalization needed, after switching on the gauge interactions in either case. This was explicitly checked by numerical simulations in the exponential (indeed, generically for any f · g = 1) case for QCD in the quenched approximation [21,26]. It was proposed in Ref. [17,25], and again demonstrated for the non-interacting fermions, that the spurious and divergent terms arising in the linear µ case can be evaluated and explicitly subtracted. Clear tests of such a proposal for the interacting case are that the continuum limit of a → 0 of the so-subtracted quark number susceptibilities should i) exhibit no additional divergences and ii) yield the same result as for the exponential µ form. In this section, we report results of our these numerical tests for the linear µ case in quenched QCD and verify that it passes both the tests, as expected. The choice of the quenched approximation was governed by the fact that the corresponding published available results [21,26] for the exponential case make it simpler to compare.
On an N 3 × N T lattice, the temperature is given by T = 1/ (N T a), where a the lattice spacing is governed by the gauge coupling β = 6/g 2 . We employ standard Wilson plaquette action for the gauge fields and use the staggered quarks for our susceptibility determinations with the corresponding Dirac matrix D(aµ) given by Here ma is the quark mass and η's are the usual staggered fermion phases. Our quenched QCD configurations were generated by using the Cabibbo-Marinari pseudo-heatbath algorithm with three SU(2) subgroup update per sweep using the Kennedy-Pendleton updating method. We chose to simulate at two different temperatures and two different quark masses on a variety of lattice sizes, as listed in the Table I, by selecting suitable β values [21,26] such that T was held constant as N T (or equivalently a) was increased (decreased). This enabled us to make a continuum limit extrapolation at both the temperatures. We quote the temperatures in the units of the critical temperature T c corresponding to the first order transition for SU(3), defined by using the order parameter, the Polyakov loop. Although we do not need it explicitly anywhere below, we mention that T c = 276(2) MeV [27] in the continuum limit, using the string tension value to be 425 MeV to set the scale. Noting from Eq. (10) that only the first derivative with respect to aµ, D ′ , is nonzero, the quark number susceptibility for this linear chemical potential action is : This expression is similar to that in Ref. [10] but without the D ′′ term which is identically zero here, as are further higher derivatives with µ. Our notation is same as in Ref. [10] to facilitate comparison of our numerical results below. The traces in the above expression were computed stochastically using Gaussian random vectors. From the Monte Carlo time evolution of the different operators that enter the susceptibility computation, including the two terms above separately, we estimate that the autocorrelation length is much less than 1000 sweeps. In order to ensure that our measurements are statistically independent, they were done on configurations [28] separated by 1000 heatbath sweeps and excluding the first 5000-10000 sweeps for thermalization. Such N conf igs configurations, which varied from 24-100, were employed to obtain the thermodynamic averages. The details of the number of random vectors and number of configurations used at each quark mass and temperature are summarized in the Table I. As discussed in the previous section, the choice of D(aµ) in Eq. (10) leads to a QNS with a term ∝ 1/a 2 for the free case. At each value of lattice cut-off, we computed numerically the coefficient for this 1/a 2 -term for non-interacting fermions on the corresponding N 3 × ∞ lattice and subtracted it from the computed values of the susceptibility in the interacting case. If there are no additional divergences in the interacting theory, one expects the continuum extrapolation in 1/N 2 T ∼ a 2 performed on these subtracted values of the susceptibilities to have a smooth limit. On the other hand, if further divergences do exist in the interacting theory then the 1/a 2 or equivalently for a fixed temperature the N 2 T dependent term would survive and increase rapidly to blow up in the continuum limit. The results for the dimensionless second order susceptibility χ 20 /T 2 , after the subtraction of the free results at 1.25 T c and 2 T c are displayed in Figure 2 for the same physical quark mass m/T c = 0.1, and in Figure 3 for the  smaller m/T c = 0.01. On general grounds, we expect χ 20 /T 2 should behave as Since the values of the second order susceptibilities reduce with increasing N T in all the plots, there is clearly no sign of any divergence in the interacting theory with c 2 (T ) = 0. This is so irrespective of the temperature, and even after lowering the quark mass by a factor of 10. Our results after fitting the data with unknown coefficients c 1 and c 3 are summarized in Table II. As seen in Figure 2, our continuum extrapolated values are in agreement with the corresponding results obtained using the exponential in µ action [26] at these temperatures. We also verified that there is no mass dependent divergent term in the expression for the second order susceptibility by performing a similar continuum extrapolation using a different bare quark mass of m/T c = 0.01, as shown in Figure 3. The second order off-diagonal susceptibility is identical in both the linear and the exponential method. Since it is zero for free fermion, no additional subtraction is necessary for the linear case. Calculations for free fermions show that the fourth order quark number susceptibility has no divergent contributions but an additional finite contribution in the continuum limit for the linear µ-case. Adopting the same procedure for it as well, we subtract the additional obtained free contribution from the corresponding (quenched) lattice QCD determination.
For the fourth order susceptibility, there are one diagonal and two off-diagonal components. Following the convention of [10], these can be written as, , as discussed above. In order to be consistent, a subtraction of such a constant from O 2 in the expressions above should also be done. It can be easily verified that such a substitution in the expressions above does not change them at all since all the free theory divergence terms arising out of it cancel in each expression. This is, of course, consistent with the fact the direct computation of the free diagonal susceptibility χ 40 has no divergence in the continuum limit for the linear µ-case either. It does have a constant a 0 term as an artifact though in the term coming from O 4 . Indeed, it is clear that due to dimensional reasons any difference between the linear and exponential case must be a constant for O 4 . Moreover, all O n for n > 4 will only differ by terms which vanish in the continuum limit, being of order a n−4 . Thus the still higher susceptibilities must all agree in the continuum limit with the exponential µ-case. Inspired by the success for the second order susceptibility, we computed the O 4 for free fermions at the same temperature and subtracted it from the value obtained for the quenched theory at that temperature in order to eliminate the a 0 term artifact. Our results for different N T are displayed in Figure 4. These results for χ 40 also show a convergence with increasing N T . Thus there are indeed no additional divergences in the continuum limit, as anticipated. Moreover, the subtraction of the unphysical artifact appears to have been done correctly on each lattice size. The continuum results for the exponential case are not available in the literature to facilitate a comparison unlike Fig. 2. The results do show a converging trend though.
For the off-diagonal susceptibilities at fourth order, the free theory artifacts are again zero so no subtraction are expected. Indeed, we observe this to be true for χ 22 in Fig. 5. For χ 31 in the right panel of the same figure, it appears somewhat difficult to draw any definite conclusions although a finite continuum limit is suggested.

IV. SUMMARY
Investigations of QCD at finite density, and in particular of the QCD critical point, gain from using the canonical Lagrange multiplier type linear chemical potential term in the fermion actions on the lattice. Preservation of exact chiral invariance on the lattice seems feasible [23] only for such a linear term for the overlap and the domain wall fermions. The higher order quark number susceptibilities needed for locating the critical point using the Taylor expansion approach are easier to compute in the linear case as well. However, it is known [15,16,19,24] since long that the linear term leads to O(1/a 2 ) divergences in the baryon number susceptibility. We have shown that such a divergence exists already in the continuum for a gas of free fermions, and therefore, lattice merely faithfully reproduces it. Using simulations of the quenched QCD with the staggered fermions, we have verified that once the free theory divergence is explicitly subtracted out, the susceptibility has no additional divergence in the continuum limit. This is only to be expected since the conserved charge, or number density, does not get renormalized in an interacting theory. Furthermore, its extrapolated value in the continuum agrees very well at two different temperatures with the similar continuum results for the fermion action using terms exponential in µ, which by construction is free of any such divergences. The higher order susceptibilities were also shown to be free of divergences in quenched QCD, and it was argued why this was to be expected. An important consequence of our result is that this would enable computation of the higher order QNS in a significantly shorter time and with better control of errors, thereby enabling a use of ratios of still higher order QNS in locating the QCD critical point and a more precise equation of state at finite baryon density. These could also be exciting for the heavy ion experiments which have already reported preliminary hints of a possible critical point and are beginning to probe the finite density region in the ongoing and future programs.