Fluctuations and Correlations of net baryon number, electric charge, and strangeness: A comparison of lattice QCD results with the hadron resonance gas model
Abstract
We calculate the quadratic fluctuations of net baryon number, electric charge and strangeness as well as correlations among these conserved charges in (2+1)flavor lattice QCD at zero chemical potential. Results are obtained using calculations with tree level improved gauge and the highly improved staggered quark (HISQ) actions with almost physical light and strange quark masses at three different values of the lattice cutoff. Our choice of parameters corresponds to a value of MeV for the lightest pseudo scalar Goldstone mass and a physical value of the kaon mass. The three diagonal charge susceptibilities and the correlations among conserved charges have been extrapolated to the continuum limit in the temperature interval . We compare our results with the hadron resonance gas (HRG) model calculations and find agreement with HRG model results only for temperatures MeV. We observe significant deviations in the temperature range and qualitative differences in the behavior of the three conserved charge sectors. At MeV quadratic net baryon number fluctuations in QCD agree with HRG model calculations while, the net electric charge fluctuations in QCD are about 10% smaller and net strangeness fluctuations are about 20% larger. These findings are relevant to the discussion of freezeout conditions in relativistic heavy ion collisions.
I Introduction
The low energy runs currently being performed at the Relativistic Heavy Ion Collider (RHIC) (1) aim at an exploration of the QCD phase diagram at nonzero temperature () and baryon chemical potential () through the measurement of fluctuations of conserved charges, , net baryon number, electric charge and strangeness. For the former, first results have been published by the STAR collaboration (2) and preliminary results on the latter two have been presented at conferences. A central goal of these experiments is to search for the existence of the QCD critical point, a second order phase transition point, that has been postulated to exist at nonvanishing baryon chemical potential in the  phase diagram of QCD (3); (4). Fluctuations in conserved charges can probe this critical point, the endpoint of a line of first order phase transitions that extends to large baryon chemical potential at nonzero quark masses. More generally, the study of fluctuations at any value of the baryon chemical potential probe thermal conditions in a medium and provide information on the critical behavior of QCD (5).
At zero baryon chemical potential, already, the analysis of fluctuations of conserved charges and their higher order cumulants provides important information about the relation between the QCD chiral phase transition at vanishing light quark masses, the crossover temperature at physical quark masses and the freezeout conditions observed in heavy ion experiments (5); (6); (7). Thus, calculations of conserved charge fluctuations at will provide unique information on freezeout conditions at the LHC where the baryon chemical potential is small, . The quadratic fluctuations of the net baryon number characterize the width of the probability distribution which has been measured at RHIC (2) and has recently been analyzed in the framework of the hadron resonance gas (HRG) model (8). Although the HRG model provides a rather satisfactory description of global hadron yields at chemical freezeout (9), its ability to describe detailed properties of strongly interacting matter such as fluctuations of conserved charges and, in particular, their higher order cumulants is not obvious. The foundation for the HRG model is given by the Dashen, Ma, Bernstein theorem (12), which shows that the partition function of strongly interacting matter can be described by a gas of free resonances, if the system is sufficiently dilute and the resonance production is the dominant part of the interaction among hadrons (13); (14). At very low temperature and high baryon number density, where nonresonant nucleonnucleon interactions become important, as well as at high temperature where strongly interacting matter undergoes a transition to the quarkgluon plasma regime and partonic degrees of freedom become dominant, the HRG model is expected to be a poor approximation to the thermodynamics of strongly interacting matter. To what extent the HRG model provides a good description of strongly interacting matter needs to be explored in detail by comparing model calculations with first principal (lattice) QCD calculations. The latter provides the complete description of QCD thermodynamics at all values of the temperature and ultimately should set the standard for the interpretation of experimental results on strong interaction thermodynamics. Unfortunately, this is, at present, not fully possible. For instance, the only satisfactory way to specify the thermal conditions in a heavy ion experiment at the time of hadronization is through the comparison of experimental data with HRG model calculations (15). It has recently been pointed out that in the future this may be overcome by comparing experimental data on conserved charge fluctuations directly with lattice QCD calculations (5). Also in order to establish such an approach more firmly it is important to understand and quantify to what extent lattice QCD calculations and HRG model calculations agree and in which temperature regime the latter provides a reasonable approximation to strong interaction thermodynamics.
Indeed, some deviations in ratios of higher order cumulants of baryon number fluctuations, calculated within the HRG model, from experimental results for cumulants of net proton fluctuations have been observed (8). Whether these deviations can be accounted for within QCD or are of more technical origin related to the restricted phase space in which experimental observations have been performed is an open question. In any case, a more detailed analysis of the thermal conditions achieved in heavy ion experiments is important. Lattice QCD calculations of fluctuations of conserved charges in equilibrium thermodynamics provide a base line for such discussions. Studies of fluctuations and higher order cumulants (16) may reveal differences between HRG model calculations and QCD thermodynamics that will appear close to criticality in the QCD phase diagram.
In this paper we will use the HRG model in its simplest version, i.e. as a sum of noninteracting, pointlike particles. This is known to provide an accurate description of a dilute, strongly interacting hadron gas (13). The need for taking into account residual interactions, for instance through the introduction of an intrinsic size of the hadrons (10) has been discussed. It has also been noted that the inclusion of higher mass resonances and an improvement in the strangeness sector of the HRG model may be needed to adequately describe pronounced features of hadron production such as the enhancement in the ratio (11). The advantage of the simplest version, however, is that the HRG model in this form is parameter free, while any improvement on this model will introduce further parameters without bringing us closer to the actual underlying theory, QCD, in a controlled way.
Quadratic fluctuations of conserved charges are closely related to quark number susceptibilities (17). Fluctuations of net baryon number, electric charge and strangeness, as well as correlations among them, have been analyzed in previous lattice QCD calculations (18); (19); (20); (21) and have also been used to characterize properties of the relevant thermodynamic degrees of freedom at low as well as high temperature (16); (22); (23). The generic forms of their temperature dependence and their scaling properties are understood in terms of universal properties of the QCD partition function and its derivatives in the vicinity of the QCD chiral phase transition (7); (20). To make use of this knowledge in a quantitative comparison with experimental results, lattice QCD calculations close to the continuum are needed.
In this paper we present an analysis of fluctuations in,
and correlations among, conserved charges using numerical calculations
in (2+1)flavor QCD at three values of the lattice
cutoff
Ii Fluctuations of conserved charges from lattice QCD; the hadron resonance gas and the ideal gas limit
To calculate fluctuations of baryon number (), electric charge () and strangeness () from (lattice) QCD we start from the QCD partition function with nonzero light (, ) and strange quark () chemical potentials. The quark chemical potentials can be expressed in terms of chemical potentials for baryon number (), strangeness () and electric charge (),
(1) 
The starting point of the analysis is the pressure given by the logarithm of the QCD partition function,
(2) 
Fluctuations of the conserved charges and their correlations in a thermalized medium are then obtained from its derivatives evaluated at ,
(3)  
(4) 
with and . Explicit expressions for the calculation of these susceptibilities in terms of generalized light and strange quark number susceptibilities are given in (20).
As all these derivatives are evaluated at , the expectation values of all net charge numbers , with (), denoting the number of particles (antiparticles), vanish, i.e., . The susceptibilities, i.e., the quadratic fluctuations of the charges, are then given by
(5) 
ii.1 The hadron resonance gas
We will compare results for fluctuations and correlations defined by Eqs. (3) and (4) with hadron resonance gas model calculations. The partition function of the HRG model can be split into mesonic and baryonic contributions,
(6)  
where the partition function for mesonic () or baryonic () particle species with mass is given by,
(7)  
Here upper signs correspond to mesons and lower signs to baryons; denotes the energy of particle , is its degeneracy factor and its fugacity is given by
(8) 
With these relations it is straightforward to calculate susceptibilities and charge correlations in the HRG model using Eqs. (3) and (4).
We note that a HRG model is defined by specifying the resonance
spectrum used to construct the partition function in
Eq. (6).
We use all hadron resonances with masses GeV
listed by the particle data group (PDG) in their 2010 summary
tables
ii.2 The ideal gas limit
In the infinite temperature limit, the grand canonical QCD partition function reduces to that of an ideal gas of quarks and gluons. In this limit, quark mass effects, including those in the strange quark sector, are negligible and we may compare our (2+1)flavor QCD calculations with a free quarkgluon gas of 3flavor QCD (StefanBoltzmann (SB) gas). This is given by (29)
(9) 
where the first two terms give the contributions of the gluon and the quark sectors for vanishing chemical potentials. After expressing the flavor chemical potentials in terms of , and as given in Eq. (1), it is straightforward to read off the ideal gas values for diagonal and offdiagonal susceptibilities. These are listed in Table 1.
1/3  0  1/3  
0  2/3  1/3  
1/3  1/3  1 
Iii Lattice calculations
In order to analyze fluctuations of conserved charges, we perform calculations on gauge field configurations generated in our study of the finitetemperature transition in (2+1)flavor QCD (27). These calculations were performed with the HISQ/tree action for three values of the lattice cutoff corresponding to lattices with temporal extent and and a spatial lattice extent . The calculations cover a temperature range, . In this temperature range the line of constant physics is defined by tuning the strange quark mass to its physical value and setting the light quark masses to , which correspond to a pion mass MeV.
It is important to note that the definition of physical quark and/or pion masses is not straightforward in calculations with staggered fermions at nonzero values of the lattice spacing due to taste symmetry breaking. Taste symmetry breaking, which is a consequence of the doubling problem in the staggered formulation, gives rise to sixteen pseudoscalar mesons corresponding to the sixteen elements of the Clifford algebra, of which only one, with taste matrix , behaves as a Goldstone particle at finite lattice spacing. Even in the chiral limit, the other 15 modes receive masses of which vanish only in the continuum limit. Thus, as a result of taste symmetry breaking, these sixteen modes contribute to observables with different masses. The same is true of all other states, but the problem is most severe for Goldstone modes. The size of these effects has been quantified for the staggered fermion discretization scheme exploited here (HISQ/tree) by measuring the masses of the sixteen taste pions and by defining a rootmeansquared mass (27),
(10) 
Taste symmetry breaking also affects the generation of the ensemble of gauge configurations on which measurements are made. To simulate the desired number of flavors, the fourthroot of the staggered fermion determinant is taken for each flavor. While this “rooting” trick corrects for the number of flavors in the continuum limit, at finite lattice spacings, the masses of states contributing to the partition function are not degenerate.
With the improved staggered fermion action (HISQ/tree) used in our calculation, taste violations, while strongly suppressed, are still large (27). In fact, our analysis showed that the HISQ/tree action has the smallest taste violations compared to the other improved staggered actions (p4, asqtad and stout) used in finite temperature calculations. Nevertheless, at values of the cutoff corresponding to the transition region ( MeV) on the three different lattices analyzed by us, the RMS masses vary from MeV on our finest lattices () to MeV on the coarsest lattice (), in contrast to the Goldstone pion mass MeV (27). We also find that, to a good approximation, the difference of and the physical pion mass is proportional to the square of the lattice cutoff at a fixed value of the temperature, . In this paper we show that the cutoff effects resulting from a much heavier RMS mass influence most strongly the electric charge fluctuations as these are most sensitive to the contributions from pions. It is worth pointing out, for comparison, that to obtain an RMS mass of about MeV achieved with the HISQ/tree action on lattices will require with the asqtad and stout actions (27); (21).
We typically analyzed 10,000–20,000 gauge field configurations per parameter set for , 10,000–14,000 configurations for and up to 6000 configurations for lattices. Measurements of all operators needed to calculate quadratic fluctuations have been performed every 10 hybrid Monte Carlo time units for and every 10 or 20 time units for lattices. Most of our quark number susceptibilities on the and lattices were calculated on GPUclusters and we used 500–1500 random source vectors for the analysis. The data were analyzed using 400 random source vectors at low temperatures and 100 at high temperatures.
When performing lattice QCD calculations at nonzero temperature, we have to control (at least) two different sources of cutoff errors. On the one hand there is the intrinsic cutoff dependence of the observables calculated at nonzero temperature at a certain value of the cutoff . We reduce these by working with treelevel improved actions in the gauge as well as the fermion sector. This improvement also insures that at high temperature the cutoff dependence of the basic operators entering our calculations is small. An additional cutoff dependence arises due to the choice of the zero temperature observable used to set the scale for all finite temperature measurements. We investigate two different scalesetting observables: the length scale extracted from the slope of the static quark potential,
(11) 
and the kaon decay constant . The temperature in these units is and . The efficacy of these observables in setting the scale is complementary in many respects (27). In particular, we note that has a mild dependence on the quark masses and is well defined even in the infinite quark mass limit. The kaon decay constant, on the other hand, has a sizeable quark mass dependence, and at lattice spacings used in this study, it is sensitive to taste symmetry violations in the hadron spectrum. To convert to physical units, we use fm (30) and the latest PDG value for the kaon decay constant, MeV. In Fig. 1, we show the difference in the two estimates of temperature over the interval relevant to the calculations performed here. The leading order correction contributing to this difference is .
As noted in (27), we also find that in the study of charge fluctuations it is advantageous to use the scale as it absorbs a significant fraction of the cutoff effects, , the cutoff effects are similar and cancel to a large extent in the ratio of hadron masses and temperature, . However, it should be stressed that any one observable ( or ) cannot eliminate cutoff effects in all observables equally well. We elaborate on this point in more detail in the appendix.
Iv Fluctuations
iv.1 Fluctuations in baryon number, strangeness, and electric charge
We start our discussion of the fluctuations of baryon number, strangeness and electric charge by summarizing the data obtained on lattices with temporal extent and in Tables 2, 3 and 4 and discussing their scaling behavior. The continuum extrapolation and comparison to the HRG model as well as the asymptotic high temperature ideal gas results are discussed in detail in the next subsection.
(MeV)  

5.900  124.17  0.0223(15)  0.1173( 7)  0.0601( 6)  0.0090( 6)  0.0067( 6)  0.0256( 2) 
6.000  138.22  0.0436(28)  0.1829(12)  0.1080(10)  0.0207(11)  0.0114( 9)  0.0437( 3) 
6.025  141.98  0.0545(22)  0.2037(16)  0.1241(14)  0.0258(12)  0.0144( 6)  0.0491( 6) 
6.050  145.83  0.0590(23)  0.2283( 7)  0.1426(11)  0.0283(14)  0.0154( 6)  0.0572( 4) 
6.075  149.79  0.0745(26)  0.2551(14)  0.1672(14)  0.0382(13)  0.0182( 6)  0.0645( 5) 
6.100  153.85  0.0881(24)  0.2885(14)  0.1983(14)  0.0472(13)  0.0205( 6)  0.0756( 4) 
6.125  158.01  0.1044(19)  0.3216(18)  0.2317(19)  0.0580(11)  0.0225( 4)  0.0868( 5) 
6.150  162.28  0.1256(22)  0.3614(31)  0.2741(26)  0.0729(15)  0.0264( 4)  0.1006( 9) 
6.175  166.66  0.1466(52)  0.3956(32)  0.3200(54)  0.0887(29)  0.0290(25)  0.1157(20) 
6.195  170.25  0.1548(32)  0.4176(20)  0.3513(30)  0.0987(23)  0.0280( 5)  0.1263( 6) 
6.215  173.90  0.1730(34)  0.4427(25)  0.3937(38)  0.1155(25)  0.0288( 5)  0.1391( 7) 
6.245  179.52  0.1881(29)  0.4666(23)  0.4422(36)  0.1324(20)  0.0278( 5)  0.1549(10) 
6.285  187.27  0.2048(68)  0.4978(30)  0.5116(64)  0.1504(50)  0.0272(32)  0.1806(31) 
6.341  198.61  0.2272(20)  0.5198(24)  0.5893(40)  0.1845(17)  0.0214( 2)  0.2024(13) 
6.354  201.34  0.2326(53)  0.5261(25)  0.6098(72)  0.1936(46)  0.0195(24)  0.2083(30) 
6.423  216.33  0.2517(14)  0.5466(12)  0.6881(23)  0.2212(13)  0.0152( 1)  0.2334( 7) 
6.488  231.33  0.2639(28)  0.5586(15)  0.7420(43)  0.2404(26)  0.0118(13)  0.2508(18) 
6.515  237.81  0.2676( 9)  0.5629( 9)  0.7576(17)  0.2472( 9)  0.0102( 1)  0.2552( 5) 
6.550  246.45  0.2672(41)  0.5678(19)  0.7787(54)  0.2528(36)  0.0072(19)  0.2629(24) 
6.664  276.43  0.2777(11)  0.5714(11)  0.8101(23)  0.2667(11)  0.0055( 1)  0.2717( 5) 
6.800  316.10  0.2806(18)  0.5718(10)  0.8248(36)  0.2722(18)  0.0042( 8)  0.2761(14) 
6.950  365.18  0.2862(27)  0.5714(15)  0.8374(57)  0.2794(26)  0.0034(14)  0.2790(21) 
7.150  440.31  0.2803(23)  0.5666(12)  0.8395(43)  0.2783(23)  0.0010(11)  0.2804(17) 
(MeV)  

6.195  127.69  0.0083(58)  0.1377(19)  0.0690(13)  0.0050(21)  0.0018(20)  0.0321( 7) 
6.245  134.64  0.0332(52)  0.1711(18)  0.0933(17)  0.0151(23)  0.0090(16)  0.0391( 7) 
6.260  136.79  0.0390(28)  0.1868(20)  0.1054(11)  0.0183(10)  0.0104(12)  0.0436( 4) 
6.285  140.45  0.0417(40)  0.2092(19)  0.1212(18)  0.0198(20)  0.0110(11)  0.0507( 3) 
6.315  144.95  0.0668(38)  0.2455(23)  0.1490(31)  0.0315(29)  0.0177( 7)  0.0588( 6) 
6.341  148.96  0.0673(49)  0.2692(38)  0.1718(39)  0.0345(29)  0.0164(12)  0.0686( 7) 
6.354  151.00  0.0705(63)  0.2908(23)  0.1878(31)  0.0359(35)  0.0173(15)  0.0760( 7) 
6.390  156.78  0.0988(54)  0.3362(37)  0.2385(40)  0.0557(31)  0.0215(11)  0.0914(10) 
6.423  162.25  0.1235(28)  0.3811(15)  0.2894(23)  0.0731(18)  0.0252( 7)  0.1081( 6) 
6.445  165.98  0.1444(27)  0.4161(30)  0.3398(19)  0.0919(15)  0.0263( 8)  0.1239(11) 
6.460  168.57  0.1610(24)  0.4287(19)  0.3623(30)  0.1034(17)  0.0288( 5)  0.1294( 9) 
6.488  173.49  0.1731(34)  0.4534(24)  0.4058(44)  0.1170(28)  0.0281( 4)  0.1444(10) 
6.515  178.36  0.1948(24)  0.4794(28)  0.4606(35)  0.1390(19)  0.0279( 5)  0.1608(11) 
6.550  184.84  0.2085(17)  0.5034(19)  0.5177(30)  0.1572(13)  0.0257( 4)  0.1803(10) 
6.575  189.58  0.2190(17)  0.5171(10)  0.5586(19)  0.1719(13)  0.0236( 2)  0.1933( 4) 
6.608  196.01  0.2291(23)  0.5321(18)  0.6052(39)  0.1874(22)  0.0209( 1)  0.2089(10) 
6.664  207.32  0.2465(22)  0.5511(13)  0.6711(35)  0.2123(20)  0.0171( 2)  0.2294( 8) 
6.800  237.07  0.2748(14)  0.5786(14)  0.7796(27)  0.2540(13)  0.0104( 1)  0.2628( 7) 
6.950  273.88  0.2901( 7)  0.5963(10)  0.8443(18)  0.2782( 7)  0.0060( 1)  0.2831( 5) 
7.150  330.23  0.2984( 5)  0.6073( 9)  0.8829(13)  0.2919( 5)  0.0032( 1)  0.2955( 4) 
(MeV)  

6.70  143.25  0.0404(135)  0.2741( 64)  0.1586( 53)  0.0227(59)  0.0088(65)  0.0679(33) 
6.74  149.03  0.0764( 67)  0.3178( 52)  0.1939( 46)  0.0348(40)  0.0208(32)  0.0796(25) 
6.77  153.48  0.1031( 68)  0.3492( 86)  0.2365( 59)  0.0581(42)  0.0225(34)  0.0892(33) 
6.80  158.05  0.1056( 83)  0.3721( 46)  0.2663( 53)  0.0632(47)  0.0212(40)  0.1015(27) 
6.84  164.31  0.1520( 42)  0.4358( 45)  0.3531( 59)  0.0975(30)  0.0272(20)  0.1278(22) 
6.88  170.77  0.1522( 64)  0.4471( 40)  0.3808( 66)  0.1009(46)  0.0257(30)  0.1399(30) 
6.91  175.76  0.1920( 46)  0.4943( 36)  0.4650( 63)  0.1361(42)  0.0280(21)  0.1645(28) 
6.95  182.59  0.2034( 55)  0.5030( 59)  0.5067( 90)  0.1506(55)  0.0264(28)  0.1780(46) 
6.99  189.64  0.2273( 73)  0.5273( 69)  0.5779( 87)  0.1820(66)  0.0227(35)  0.1979(48) 
7.03  196.91  0.2470( 49)  0.5483( 36)  0.6319( 87)  0.2006(45)  0.0232(24)  0.2157(34) 
7.10  210.21  0.2547( 33)  0.5671( 27)  0.7004(106)  0.2213(38)  0.0167(16)  0.2395(30) 
7.15  220.15  0.2650( 25)  0.5744( 26)  0.7389( 44)  0.2366(32)  0.0142(12)  0.2512(24) 
7.28  247.91  0.2822( 19)  0.5902( 15)  0.8063( 37)  0.2646(22)  0.0088( 8)  0.2708(15) 
Figure 2 shows results for the baryon number susceptibility with the temperature scale set using (left hand panel) and (right hand panel).
In both cases, we show the results from the HRG model including all resonances with mass GeV. The noticeable differences between the left and right hand panels of Fig. 2 are due to setting of the temperature scale and we find that the cutoff effects are smaller when a scale based on is used. As pointed out above, this feature has also been noted in the analysis of chiral observables (27); (31). In principle, both scales should lead to identical results in the continuum limit, however, the continuum extrapolation is much better controlled when the cutoff effects are small. We therefore use the temperature scale in the rest of the paper unless stated otherwise.
In Fig. 3, we show results for fluctuations of net strangeness (left) and net electric charge (right) using to set the scale for the temperature. The figure also shows the corresponding HRG results. Both, the strangeness and the baryon number (see Fig. 2) fluctuations agree with HRG results for temperatures below MeV and show small cutoff effects. The electric charge fluctuations deviate significantly from HRG results at all temperatures and show large cutoff dependence. These large cutoff effects in can mostly be explained as due to discretization errors in the lattice hadron spectrum. The dominant contribution to at low temperatures comes from pions, while receives the leading contribution from kaons and from nucleons. As discussed in Section III, the pion spectrum is strongly affected by taste symmetry violations in the staggered formulation, and it has been shown in earlier work that the agreement in between lattice QCD results and HRG model calculations can be improved using a distorted spectrum in the HRG analysis (32); (34).
To confirm this observation, we show in Fig. 3 (right panel) the HRG results obtained after replacing the physical pion mass by the RMS pion mass relevant to our calculations on the , 8 and 12 lattices, respectively. These cutoff effects leading to the spectral distortion vary not only with but also with temperatures at fixed . We therefore parametrize the cutoff dependence of , at each , using a cubic polynomial fit to the data given in (27). This allowed us to estimate the modified HRG result as a function of temperature, as shown in Fig. 3 by short solid lines (colored lines). The data confirm that below MeV, the numerical results for are well described by this minimally modified HRG model and the major part of the difference is indeed due to the taste symmetry breaking effects in the pion sector. Since all other heavier states contributing to the HRG model are assumed to take on their physical values, the fluctuations of strangeness and baryon number are not influenced by this modification of the HRG model. We discuss these features, together with continuum extrapolations, for the data shown in Fig. 2 and 3 in more detail in the following subsection.
It is worthwhile to clarify our discussion of the comparison of lattice QCD data with the HRG model. When we say that the HRG is a good approximation to QCD, we refer to the value of the susceptibility, as is traditional. At the temperature where this agreement fails, we observe that not only the value but the slope also deviates significantly. Our data also indicate that the curvature starts to deviate 20–30 MeV earlier depending on the observable. These derivatives of the susceptibility are related to higher moments, which are increasingly less well captured in the HRG analysis and have not been calculated in our lattice simulations.
iv.2 Continuum extrapolation and approach to the hadron resonance gas estimates
In this section, we analyze cutoff effects at fixed values of the temperature for three values of the lattice spacing, and . We perform this analysis for temperature scales defined in terms of both and in order to quantify systematic effects at a given lattice spacing and to demonstrate consistency between the estimates in the continuum limit. Our simulations on the lattices with the three different values have not been done at the same values of the temperature. As a result, in order to perform continuum extrapolations at fixed temperature, we have used cubic spline interpolations of our data throughout this paper. We propagate errors on the spline parameters. In addition to estimating statistical errors based on our entire data sample, we also use the difference of spline interpolations performed on two independent subsets constructed by choosing every second Tpoint (even and odd points), as an estimate for the systematic errors. The final error on our data is obtained by adding this error estimate and the purely statistical error obtained from the full data set in quadrature.
We first discuss the continuum extrapolation for the three susceptibilities shown in Figs. 2 and 3 at high temperature, MeV. In this regime, cutoff effects are generally small, which, to some extent, is due to the fact that our numerical calculations have been performed with an improved action with small cutoff dependence of thermodynamic observables in the infinite temperature ideal gas limit (33). In Fig. 4, we show continuum extrapolations for all three susceptibilities, , at three representative values of the temperature, , and MeV. In each case, and for both temperature scales, and , the fits show that the cutoff effects are consistent with corrections and, over the limited range of , all three susceptibilities can be extrapolated to the continuum with an Ansatz that includes corrections linear in . The continuum extrapolated results obtained with the two temperature scales agree within errors, and the results obtained on the lattices are a good approximation to these. This also is true for lower temperatures. However, in this case extrapolations linear in are no longer sufficient for the electric charge and strangeness fluctuations. Systematic effects at start to become important. This is evident from the data sets at MeV, which are shown in Fig. 4 as well. We note that also at this temperature, which is the lowest temperature for which we perform continuum extrapolations, extrapolations based on the and temperature scales are in good agreement. Having demonstrated consistency of the continuum estimate obtained using and , we, as stated previously, use the scale from in the rest of the paper because the slope in the fits is smaller.
The data for the net charge fluctuations in the temperature interval 120–250 MeV, results of the linear extrapolation for , and quadratic extrapolations for and are shown in Figs. 5 and 6. In Fig. 6(right) we also show the ratio of net baryon number and electric charge fluctuations. The continuum extrapolation shown for this quantity has been obtained from the corresponding extrapolations for and .
Continuum extrapolations in the crossover and the low temperature regions require additional considerations because the three different conserved charge susceptibilities show different sensitivities to cutoff effects. In order to quantify differences from the HRG model results in this temperature regime, and in order to clarify the extent to which the HRG model provides a good description of QCD results, we analyze the ratios , and , in Fig. 7. We find that, while baryon number fluctuations start to agree with HRG model results for MeV, the net strangeness fluctuations become larger than the HRG values for temperatures below MeV and then approach the HRG values from above at MeV. At MeV, the differences are still (1020)%.
The electric charge fluctuations show much larger deviations from the HRG model as is evident from Fig. 7. In particular, below MeV, the cutoff dependence in is large and extrapolations including just leading order corrections fail. As discussed in Section IV.1, this, to a large extent, is due to the severe cutoff dependence of the pion spectrum, , the anomalously large RMS pion mass suppresses fluctuations in the electric charge and has a much smaller effect on the baryon and strangeness charges. In short, a continuum extrapolation without including the effects of taste symmetry breaking is insufficient.
The distorted HRG model, which modifies the log of the partition function by replacing by in the pion contribution, , however, does describe the data well. In general, the HRG model defined in Eq. 7 suggests that, in this temperature regime, cutoff effects in any quantity may be accounted for by an exponential Ansatz of the form
(12) 
which, at high temperatures where cutoff effects become small, reduces to the linear fit in , i.e., and also incorporates the next to leading order quadratic corrections . We, therefore, analyze our data for in the transition region, , from the low to high temperature phase of (2+1)flavor QCD using fits linear and quadratic in as well as the exponential Ansatz given in Eq. (12). With our current statistical accuracy we are, however, not sensitive to cutoff effects beyond . In fact, all our fits performed with the exponential ansatz are consistent within errors with fits based on the quadratic ansatz. We, therefore, do not discuss the exponential fits any further in this paper. In the case of we find that linear and quadratic fits agree within errors and lead to a less than unity in the entire range of temperatures MeV. We therefore use the linear fits to perform continuum extrapolations for . For strangeness and electric charge susceptibilities we use quadratic extrapolations in the entire temperature range, although, as discussed above, we do not observe systematic differences between linear and quadratic extrapolations for MeV. Using the latter for our continuum extrapolations, however, leads to more conservative error estimates.
A summary of our continuum extrapolations for and the data in the low temperature region is also shown in Fig. 7. Due to taste symmetry breaking, the data show significant dependence on for MeV. To understand this cutoff effect we compare, in the bottom panel in Fig. 7, results for with a modified HRG model in which the physical pion mass is replaced by (i) the dependent RMS mass, HRG()/HRG(), and (ii) a pion mass of MeV, HRG( MeV)/HRG(), corresponding to the light quark mass actually used in our calculations. Since the line showing HRG( MeV)/HRG() is much closer to unity compared to the other three, we confirm that the errors due to simulating at this slightly heavier pion mass are significantly smaller than the cutoff effects leading to a much heavier RMS mass.
In the interesting temperature range relevant to the discussion of freezeout conditions in heavy ion collisions, , we find that the continuum extrapolated electric charge fluctuations are (1020)% smaller than even the modified HRG model calculation with MeV. For temperatures MeV, the data start to agree with the modified HRG results with and continuum extrapolations using the quadratic Ansatz start to agree with the HRG result.
Strangeness fluctuations on the other hand, both for the data and the continuum extrapolated values, stay systematically above the hadron resonance gas result in the temperature range 150–190 MeV. We, therefore, expect this feature to survive the continuum extrapolation. Below MeV, fluctuations in the strangeness charge show an dependence, which is most likely again due to taste symmetry violations. More data are required to study this issue further.
Net baryon number fluctuations are consistent with HRG model results for temperatures below MeV, although statistical errors on our data set make the quantification of possible deviations from the HRG result in this temperature range difficult. For larger values of the temperature, the estimates and errors in them grow progressively smaller. Our linear extrapolations suggest that in the temperature interval , deviations from the HRG model calculations are at most 10%. A confirmation of this in our data, through the inclusion of quadratic corrections, however, requires better statistics. Data obtained with the stout action (21) also suggest that the stays close to the HRG model result in this temperature range.
In order to reflect the influence of systematic effects on our continuum extrapolation we varied fit ranges and distribution of knots in the smooth spline interpolation. Moreover, in order to account for possible underestimates of errors on the individual data points, we divided our data samples into two independent subsets consisting of even and odd Tvalues by rankorder. We used the differences in these fits and the fit to the full data sample as an additional error on our spline interpolations.
To summarize, our continuum extrapolated values of are given in Table 5. Extrapolations of used the Ansatz, Eq. (12) truncated at , i.e., we used the quadratic extrapolations, whereas for the exponential was truncated at , i.e., we used the linear extrapolations. For we also find that the quadratic extrapolated values agree with the exponential Ansatz. Our extrapolated results are in good agreement with the recently published analysis using the stout action (21).
[MeV]  

150  0.0790( 57)  0.3736(126)  0.2358(225)  0.0415( 60)  0.0187( 29)  0.0929( 61) 
155  0.1026( 77)  0.3998(158)  0.2770(260)  0.0587( 69)  0.0221( 24)  0.1031( 68) 
160  0.1247( 98)  0.4258(198)  0.3208(289)  0.0767( 75)  0.0242( 23)  0.1173( 70) 
165  0.1470(109)  0.4513(232)  0.3645(322)  0.0949( 81)  0.0259( 21)  0.1332( 78) 
170  0.1662(111)  0.4734(240)  0.4067(334)  0.1115( 84)  0.0269( 21)  0.1491( 84) 
175  0.1849( 99)  0.4964(222)  0.4567(314)  0.1304( 79)  0.0270( 19)  0.1655( 83) 
180  0.2024( 74)  0.5118(179)  0.5013(259)  0.1497( 66)  0.0263( 17)  0.1782( 77) 
185  0.2182( 56)  0.5239(144)  0.5449(215)  0.1682( 57)  0.0249( 19)  0.1892( 76) 
190  0.2326( 52)  0.5379(115)  0.5888(181)  0.1860( 56)  0.0233( 22)  0.2005( 72) 
195  0.2445( 55)  0.5522( 84)  0.6289(210)  0.2005( 60)  0.0220( 21)  0.2122( 68) 
200  0.2524( 46)  0.5638( 68)  0.6623(222)  0.2116( 52)  0.0205( 21)  0.2235( 70) 
205  0.2570( 41)  0.5717( 57)  0.6884(210)  0.2200( 38)  0.0186( 18)  0.2335( 67) 
210  0.2604( 38)  0.5770( 55)  0.7113(194)  0.2267( 31)  0.0168( 15)  0.2424( 62) 
215  0.2637( 36)  0.5792( 63)  0.7282(155)  0.2335( 32)  0.0152( 14)  0.2486( 51) 
220  0.2676( 35)  0.5809( 84)  0.7438(148)  0.2397( 35)  0.0140( 12)  0.2537( 57) 
225  0.2713( 28)  0.5824( 63)  0.7565(113)  0.2456( 34)  0.0128( 12)  0.2576( 56) 
230  0.2749( 27)  0.5841( 64)  0.7672(116)  0.2511( 33)  0.0119( 12)  0.2606( 56) 
235  0.2784( 23)  0.5855( 56)  0.7760(105)  0.2564( 28)  0.0109( 10)  0.2628( 48) 
240  0.2819( 20)  0.5873( 46)  0.7847( 92)  0.2613( 24)  0.0101( 8)  0.2649( 39) 
245  0.2852( 20)  0.5890( 45)  0.7934( 98)  0.2664( 23)  0.0093( 8)  0.2667( 37) 
250  0.2885( 20)  0.5907( 45)  0.8020( 98)  0.2710( 23)  0.0085( 7)  0.2688( 34) 
Lastly, in Fig. 6 (right), we show the ratio of the net baryon number and electric charge fluctuations, . It approaches the HRG model result from above and starts to agree with it for MeV. The continuum extrapolation here is based on the linear and exponential extrapolations for and , respectively. In the region of interest to heavy ion phenomenology, this ratio varies between (at 160 MeV) and (170 MeV). Thus, fluctuations in net electric charge could be 3–4 times larger than in the net baryon number in the vicinity of the freezeout temperature.
V Correlations
Probes of the structure of QCD at finite temperature include correlations among different conserved charges. These correlations show characteristic changes in the crossover region between the low and high temperature phases of QCD, which are correlated with changes in the relevant degrees of freedom. They also provide insight into the applicability of HRG model calculations at low temperatures. The change in correlations between baryon number and electric charge, , is expected to be particularly striking as one goes from the low to the high temperature phase. At low temperatures this correlation is dominated by the contribution of protons plus antiprotons. Consequently, within the HRG model it rises exponentially with temperature in this region. In the high temperature limit of (2+1)flavor QCD, however, vanishes as the quarks become effectively massless, (), and the weighted sum of the charges of up, down and strange quarks vanishes. Results for shown in Fig. 8 are consistent with this picture.
The correlations of strangeness with the baryon number and the electric charge, and , are sensitive to changes in the strangeness degrees of freedom (22); (35); (36). Results for the temperature dependence of these correlations are shown in Fig. 9. They approach the StefanBoltzmann value, , of a massless three flavor quark gas at high temperatures. As observed in the case of the quadratic fluctuations, on decreasing the temperature towards the transition region, these correlations first overshoot the HRG model result and then approach HRG value from above at about MeV. This overshoot is more pronounced for than for . Also shown in Fig. 9 are continuum extrapolations (bands) which in the baryon sector, i.e., for BS and BQ correlations, are based on fits linear in , whereas in the meson sector, i.e., for QS correlations, quadratic corrections are also taken into account. This reflects the larger sensitivity of the latter to taste violations that also has been observed for the quadratic strangeness and electric charge fluctuations.
In the isospin symmetric case considered in this study, the flavor correlations and are equal. Also, the two correlations and are related to each other through the quadratic strangeness fluctuations, i.e., (36). One can then write the following relationships between the charge correlations and quarkflavor fluctuations:
(13) 
At high temperatures because receives perturbative contributions only at and is therefore small (23); (37). On the other hand, corrections to from the ideal gas limit are dominant as they are .
All three charge correlations show significant deviations from the ideal gas limit even at twice the transition temperature (see data in Tables 2, 3 and 4). These deviations are due to large contributions of flavor fluctuations, such as to discussed above. The leading order perturbative correction can be eliminated by forming suitable ratios that can be used to analyze experimental data on charge fluctuations (22); (35); (36),
(14) 
At high temperature the deviations from the ideal gas value of unity are now due to , for example, . Data for is shown in Fig. 10(left) and, in comparison to the quadratic strangeness fluctuations (Fig. 3) or baryon number strangeness correlations (Fig. 9), the approach to the ideal gas limit is much more rapid. This shows that the flavor correlation is already small for . It is, however, large in the vicinity of the transition temperature.
The behavior of the third ratio that one can analyze, , is somewhat different as leading order perturbative corrections do not cancel completely due to differences in the light and strange quark masses. Consequently, the approach to the ideal gas limit is slower as can be seen in Fig. 10(right).
Vi Discussion and Conclusion
We have analyzed quadratic fluctuations and correlations among conserved charges in (2+1)flavor QCD. We find that as the temperature is decreased from the high temperature phase, the net baryon number fluctuations start to agree with the hadron resonance gas model below MeV while for electric charge fluctuations this happens only below MeV. The fluctuations of net strangeness overshoot the HRG model values at MeV. In the temperature range relevant to the discussion of chemical freezeout in heavy ion collisions, 160–170 MeV, strangeness fluctuations are systematically larger than the HRG model result by about 20%. These detailed differences between QCD calculations and HRG model results should become manifest when experimental data on the probability distributions of net charge fluctuations (2) is analyzed. In fact, quadratic fluctuations characterize the bulk structure of these distributions which, in the Gaussian approximation at , is given by
(15) 
where denotes the temperature at the time of chemical freezeout and is the freezeout volume. The results presented here suggest that the largest deviations from HRG model calculations occur in the probability distributions for electric charge and strangeness fluctuations. We give a summary of the fluctuations in conserved charges for temperatures in the transition region in Table 6.
In the phenomenologically interesting temperature regime , continuum extrapolated results for , are smaller than the HRG model results by about 10–20%. Even though at temperatures below MeV, estimates for electric charge fluctuations have large systematic errors due to the distortion of the light meson spectrum in all staggered formulations, our analysis shows that these effects can be taken into account when performing continuum extrapolations. We show that the resulting continuum estimates in the transition region lie below the HRG estimates even after corrections accounting for the distorted pion spectrum have been applied. Thus, our conclusion is that at the highest RHIC energies and at the LHC, the width of probability distribution for the net electric charge should be narrower than what HRG model calculations would suggest since in Eq. 15 is smaller.
In the case of the net baryon number fluctuations, deviations from HRG model results start to become statistically significant only for MeV and will therefore be hard to quantify. We also find that the ratio of the net baryon number and electric charge fluctuations, presented in Fig. 6, approaches the HRG model result from above and starts to agree with it for MeV. In the transition region this ratio is
(16) 
i.e., fluctuations in net electric charge are expected to be about three to four times larger in the vicinity of the freezeout temperature in heavy ion collisions than net baryon number fluctuations. It is worth noting that a comparison of with the ratio of proton to net charge fluctuations, which is accessible in heavy ion collisions, will allow us to relate fluctuations in the proton number to the fluctuations of the conserved net baryon number (38).
Finally, we point out that the continuum extrapolated results presented here, and summarized in Table 6 in a temperature regime relevant to the freezeout conditions in heavy ion collisions, do not rely on any uncertainties in the determination of the QCD crossover temperature or its characterization through different observables as discussed in (27). Systematic errors on the temperature values listed in the first row of Table 6 can only come from uncertainties in the zerotemperature observable used to determine the temperature scale. We estimate these uncertainties to be less than MeV in our calculation (27).
[MeV]  

155  1.049(79)  0.924(36)  1.240(116)  1.353(159)  0.804(86)  1.139(74) 
160  1.020(80)  0.895(41)  1.235(111)  1.384(135)  0.717(67)  1.144(68) 
165  0.972(72)  0.861(44)  1.212(106)  1.356(116)  0.633(51)  1.150(67) 
170  0.898(60)  0.818(41)  1.171( 96)  1.280( 96)  0.544(42)  1.144(64) 
Acknowledgments
This work has been supported in part by contracts DEAC0298CH10886, DEAC5207NA27344, DEFC06ER41446, DEFG0291ER40628, DE FG0291ER40661, DEKA140102 with the U.S. Department of Energy, and NSF grants NSF1067881, PHY0903571, PHY0857333, PHY0757035, PHY0757333 and PHY0703296, the Bundesministerium für Bildung und Forschung under grant 06BI9001, the Gesellschaft für Schwerionenforschung under grant BILAER and the Deutsche Forschungsgemeinschaft under grant GRK881, and the EU Integrated Infrastructure Initiative HadronPhysics2. The numerical simulations have been performed on BlueGene/L computers at Lawrence Livermore National Laboratory (LLNL), the New York Center for Computational Sciences (NYCCS) at Brookhaven National Laboratory, US Teragrid (Texas Advanced Computing Center), and on clusters of the USQCD collaboration in JLab and FNAL. We thank Anton Andronic, Peter BraunMunzinger and Krzysztof Redlich for discussions and for information on their current version of the hadron resonance gas data table.
Appendix: Choice of temperature scale and the hadron spectrum
In this appendix we discuss the effect of setting the scale using different observables with massdimension one calculated in zerotemperature simulations. We denote lattice observables measured in units of the lattice spacing and calculated at a given by . In the limit , these approach their physical value as
(17) 
where is the QCD scale, is the function, and only the leading correction , relevant to our treelevel improved staggered formulation, has been retained.
Consider using the observable to define the temperature scale . On a lattice of temporal size , this temperature is given in terms of as
(18) 
The ratio of any two such temperature scales is then given by
(19)  
Using Eqs. (18) and (19) we can express the observable in terms of the temperature scale obtained from as
(20) 
This shows that if an observable of interest () has a cutoff dependence similar to that of observable () used to determine the temperature scale , i.e., , then , as an estimate of , has small cutoff effects. Since all the need not be small, improving the scaling behavior of one observable does not, in general, imply improvement in all observables.
In the low temperature regime of QCD, the relevant degrees of freedom are hadrons with masses . If the hadron resonance gas is a good approximation in this regime, continuum extrapolations of lattice data are better controlled if a temperature scale is chosen such that all the lattice estimates of have small cutoff dependences. We find that using a hadronic observable such as improves the scaling behavior of the susceptibilities and correlations between charges as shown in Fig. 4.
The pion sector is, however, different and has enhanced cutoff effects due to taste symmetry breaking. One does not, therefore, expect to absorb all these effects with a choice of the temperature scale. For this reason we had to modify the HRG analysis to compare with lattice data for charge fluctuations which are dominated by contributions from the pions.
Footnotes
 Preliminary results of this work had been presented at Quark Matter 2011 (24) and PANIC 2011 (25).
 In this summary table there are a few three starred resonances listed which do not have a known spin assignment. For these we use the minimal degeneracy factors. We also checked that the inclusion of some known heavier resonances as well as the inclusion of charmed hadrons does not alter the picture presented here. Moreover, we checked the stability of HRG results by reducing the mass cutoff from 2.5 GeV to 2.0 GeV. This alters the relevant observables discussed here by at most 2% at MeV.
References

H. Caines [STAR Collaboration],
[arXiv:0906.0305 [nuclex]];
G. Odyniec [STAR Collaboration], J. Phys. G 35, 104164 (2008).  M. M. Aggarwal et al. [STAR Collaboration], Phys. Rev. Lett. 105, 022302 (2010).
 M. Asakawa and K. Yazaki, Nucl. Phys. A 504, 668 (1989).

J. Berges and K. Rajagopal,
Nucl. Phys. B 538, 215 (1999);
A. M. Halasz, A. D. Jackson, R. E. Shrock, M. A. Stephanov, and J. J. M. Verbaarschot, Phys. Rev. D 58, 096007 (1998).  F. Karsch, arXiv:1202.4173 [heplat], presented at the International Conference ”Critical Point and Onset of Deconfinement  CPOD 2011”, Wuhan, November 711, 2011.
 F. Karsch and K. Redlich, Phys. Lett. B 695, 136 (2011).
 B. Friman, F. Karsch, K. Redlich, and V. Skokov, Eur. Phys. J. C 71, 1694 (2011), [arXiv:1103.3511 [hepph]].
 P. BraunMunzinger, B. Friman, F. Karsch, K. Redlich, and V. Skokov, [arXiv:1107.4267 [hepph]].
 P. BraunMunzinger, K. Redlich, and J. Stachel, In Hwa, R.C. (ed.) et al.: Quark gluon plasma 491599, [nuclth/0304013].
 A. Andronic, P. BraunMunzinger, J. Stachel, and M. Winn, [arXiv:1201.0693 [nuclth]].
 A. Andronic, P. BraunMunzinger, and J. Stachel, Phys. Lett. B 673, 142 (2009), [Erratumibid. B 678, 516 (2009)] [arXiv:0812.1186 [nuclth]].
 R. Dashen, S.K. Ma and H.J. Bernstein, Phys. Rev. 187, 345 (1969).
 R. Venugopalan and M. Prakash, Nucl. Phys. A 546, 718 (1992).
 F. Becattini, J. Phys.: Conf. Ser. 5, 175 (2005).
 J. Cleymans, H. Oeschler, K. Redlich and S. Wheaton, Phys. Rev. C 73, 034905 (2006).
 S. Ejiri, F. Karsch, and K. Redlich, Phys. Lett. B 633, 275 (2006).
 S. A. Gottlieb, W. Liu, D. Toussaint, R. L. Renken, and R. L. Sugar, Phys. Rev. Lett. 59, 2247 (1987).
 R. V. Gavai, S. Gupta, and P. Majumdar, Phys. Rev. D 65, 054506 (2002), [heplat/0110032].
 C. Bernard et al. [ MILC Collaboration ], Phys. Rev. D 71, 034504 (2005). [heplat/0405029].
 M. Cheng, P. Hegde, C. Jung, F. Karsch, O. Kaczmarek, E. Laermann, R. D. Mawhinney, C. Miao, P. Petreczky, C. Schmidt, and W. Soeldner, Phys. Rev. D 79, 074505 (2009).
 S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. Szabo, [arXiv:1112.4416 [heplat]].
 V. Koch, A. Majumder, and J. Randrup, Phys. Rev. Lett. 95, 182301 (2005).
 C. R. Allton, M. Doring, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann, and K. Redlich, Phys. Rev. D 71, 054508 (2005), [heplat/0501030].
 S. Mukherjee, J. Phys. G G 38, 124022 (2011), [arXiv:1107.0765 [nuclth]].
 P. Hegde (for the HotQCD Collaboration), arXiv:1110.5932 [heplat].
 E. Follana et al. [HPQCD collaboration and UKQCD collaboration], Phys. Rev. D 75 054502 (2007), [arXiv:heplat/0610092].
 A. Bazavov et al. [HotQCD Collaboration], arXiv:1111.1710 [heplat].
 K. Nakamura et al. (Particle Data Group), J. Phys. G 37, 075021 (2010).
 See for instance: J. Kapusta and C. Gale, Finite Temperature Field Theory, Cambridge University Press, 2006.
 A. Bazavov, et al. [ MILC Collaboration ], PoS LATTICE2010, 074 (2010), [arXiv:1012.0868 [heplat]].
 Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg, and K. K. Szabo, JHEP 0906, 088 (2009).
 F. Karsch, K. Redlich, and A. Tawfik, Phys. Lett. B 571, 67 (2003).
 U. M. Heller, F. Karsch, and B. Sturm, Phys. Rev. D 60, 114502 (1999), [heplat/9901010].
 P. Huovinen and P. Petreczky, Nucl. Phys. A 837, 26 (2010).
 S. Jeon and V. Koch, Phys. Rev. Lett. 85, 2076 (2000).
 A. Majumder and B. Muller, Phys. Rev. C 74, 054901 (2006), [nuclth/0605079].
 J.P. Blaizot, E. Iancu, and A. Rebhan, Phys. Lett. B 523, 143 (2001).
 M. Kitazawa and M. Asakawa, Phys. Rev. C 85, 021901 (2012), [arXiv:1107.2755 [nuclth]].