Equation of State for physical quark masses
Abstract
We calculate the QCD equation of state for temperatures corresponding to the transition region with physical mass values for two degenerate light quark flavors and a strange quark using an improved staggered fermion action (p4action) on lattices with temporal extent . We compare our results with previous calculations performed at twice larger values of the light quark masses as well as with results obtained from a resonance gas model calculation. We also discuss the deconfining and chiral aspects of the QCD transition in terms of renormalized Polyakov loop, strangeness fluctuations and subtracted chiral condensate. We show that compared to the calculations performed at twice larger value of the light quark mass the transition region shifts by about MeV toward smaller temperatures.
pacs:
11.15.Ha, 11.10.Wx, 12.38Gc, 12.38.MhI Introduction
First calculations of the equation of state (EoS) of hot strongly interacting matter date back to the early 80s su2 ; su3 . These early purely gluonic calculations have been improved steadily by including contributions from dynamical quark degrees of freedom (for recent reviews see detar_lat08 ; petreczky_sewm06 ; petreczky_qm09 ). In the most recent calculations the equation of state has been evaluated for 2+1 flavor QCD, i. e. in QCD with one strange quark and two light () quarks using various improved staggered fermion actions stout ; MILCeos ; ourEoS ; hotqcd_eos . The most extensive calculations of the EoS have been performed with p4 and asqtad staggered fermion formulations on lattices with temporal extent MILCeos ; ourEoS and hotqcd_eos . These actions improve both the flavor symmetry of the staggered fermions as well as the quark dispersion relations. The latter insures that thermodynamic observables are improved at high temperatures and thus have only a small cutoff dependence in this regime. The stoutlink action, which has been used for the calculation of the EoS on lattices with temporal extent stout , only improves the flavor symmetry of the staggered fermions and therefore has the same large discretization errors at high temperatures as the standard staggered fermion formulation. Indeed, the calculations of the EoS with the stoutlink action on lattices with temporal extent and reflect the expected strong cutoff dependence above the transition temperature stout .
While at high temperatures the masses of the relevant degrees of freedom, quarks and gluons, are small compared to the temperature scale, this is not the case at low temperatures and in the transition region. One thus may expect that at these temperatures thermodynamic observables are more sensitive to the quark masses, which control the mass of the light pseudoscalars and eventually are responsible for the occurrence of a true phase transition in the chiral limit. Calculations with p4 and asqtad actions have sofar been performed using light quark masses () which are one tenth of the strange quark mass () and correspond to a pseudoscalar Goldstone mass^{1}^{1}1In calculations with staggered fermions flavor symmetry is broken at nonvanishing values of the lattice spacing . As a consequence only one of the pseudoscalar mesons has a light mass that is proportional to and vanishes in the chiral limit at fixed . Full chiral symmetry with the correct Goldstone pion multiplet is recovered only for . For an estimate of the remaining flavor symmetry violations in spectrum calculations with the p4action see FlavBreak . of and MeV respectively hotqcd_eos . The calculations with the stoutlink action have been performed at the physical value of the light quark mass.
The purpose of this paper is to investigate the quark mass dependence of the EoS by calculating it with the p4action for physical values of the (degenerate) light quark masses. The calculational procedure used in this work closely follows that used in our previous calculations at ourEoS . The paper is organized as follows. In the next section we discuss the parameters and some technical details of the numerical calculations, including the choice of the quark masses and the determination of the lattice spacing that fixes the temperature scale, . In section III we show the QCD equation of state and compare it with the resonance gas model. In section IV we discuss chiral and deconfining aspects of the QCD transition. Finally section V contains our conclusions. In the appendix we give further numerical details of our calculations.
Ii Zero temperature calculations
We have performed calculations with the p4action for several values of the gauge coupling in the region of the finite temperature crossover for lattices with temporal extent . The values of the quark masses and gauge coupling are given in Table 1. Simulations were performed with the Rational Hybrid MonteCarlo (RHMC) algorithm with Hasenbush preconditioning and different time steps in the molecular dynamic (MD) evolution for gauge fields, strange quarks and light quarks. We used 5 time steps for strange quark updates per each light quark update. The gauge fields were updated 10 times more frequently than the strange quarks. The length of the MD trajectory was and the acceptance rate about . The zero temperature calculations were performed on lattices.
We would like to perform calculations of thermodynamic quantities keeping the physical values of the strange quark mass as well as the light quark masses fixed. As a starting point for our calculations we used the same finetuned values of the bare strange quark mass as in ourEoS but bare light quark masses which are twice smaller than in Ref. ourEoS , i.e. the ratio of strange to light quark mass was chosen to be .
At each beta value we then have determined the static quark potential. The static potential and the corresponding scales and defined as
(1) 
were calculated from smeared Wilson loops as described in Ref. ourTc . Here we used 30 levels of APE smearing of the spatial links to get good signals for the Wilson loop expectation values. The values of and /a are given in Table 1. The comparison with previous results given in Table I of ourEoS reveals that the new values of the scale parameters in lattice units differ from the old ones obtained at by less than 1 %, with deviations which are not systematic. Thus, the smaller values of the light quark masses do not change the lattice spacing in physical units beyond the errors, which are about (c.f. Table 1). As the consequence the values of the temperature in units of have errors of about . To give the temperature in physical units, e.g. MeV, we use the value fm. Thus the temperature scale in physical units has an overall uncertainty of about . However, since the same uncertainty is present in the calculations performed at twice larger quark mass it does not show up in the comparison of thermodynamic quantities obtained at two different quark masses.
To remove the additive divergent constant in the potential, following Ref. ourEoS we normalized it to the string form of the static quark potential, , at distance . This amounts to renormalize the temporal gauge links with . The resulting multiplicative renormalization factors, , are also given in Table 1 and will be used to define renormalized Polyakov loops that are discussed in section IV.
Having determined the scale we extracted pseudoscalar meson masses using wall sources in the calculation of meson propagators. It turned out that the mass in lattice units is, with 1 2 % accuracy, the same as in ourEoS . Thus, a readjustment of the line of constant physical mass corresponding to our new and smaller light quark masses was not necessary. In fact, in the present calculations we find that our quark mass values define a line of constant physics characterized by the following relations
(2) 
Using fm, as determined in Ref.gray , we get MeV, MeV and^{2}^{2}2 A physical value for the mass can be obtained from the relation MeV. MeV. This means that both the light quark masses and the strange quark mass are very close to their physical values. Furthermore, in the entire parameter range covered by our thermodynamic calculations deviations of the meson masses from the above values are less than .
3.430  0.0370  2610  2.6698(101)  1.8378(184)  0.1415(2)  0.4423(3)  0.5992(4)  1.5124(90) 

3.460  0.0313  3220  2.9377(90)  2.0000(202)  0.1274(4)  0.3994(4)  0.5427(2)  1.5219(71) 
3.490  0.0290  3460  3.2188(70)  2.2115(101)  0.1213(15)  0.3775(5)  0.5117(10)  1.5277(26) 
3.500  0.0253  3030  3.3251(114)  2.2550(159)  0.1110(5)  0.3491(4)  0.4743(3)  1.5292(50) 
3.510  0.0260  3040  3.4021(141)  2.3116(105)  0.1121(10)  0.3510(9)  0.4770(6)  1.5289(55) 
3.520  0.0240  2980  3.5167(68)  2.4050(106)  0.1056(5)  0.3322(6)  0.4525(3)  1.5303(20) 
3.530  0.0240  2450  3.5927(101)  2.4723(81)  0.1045(14)  0.3300(17)  0.4500(11)  1.5291(48) 
3.540  0.0240  3360  3.6802(138)  2.5090(94)  0.1040(4)  0.3251(13)  0.4432(11)  1.5272(41) 
3.545  0.0215  3090  3.7513(67)  2.5840(75)  0.0947(7)  0.3061(3)  0.4178(2)  1.5298(23) 
3.560  0.0205  3010  3.8790(67)  2.6732(81)  0.0935(12)  0.2941(13)  0.4024(11)  1.5275(18) 
3.585  0.0192  1400  4.1501(118)  2.8503(106)  0.1049(6)  0.3258(17)  0.4434(11)  1.5275(66) 
3.600  0.0192  4080  4.3033(181)  2.9351(73)  0.0888(20)  0.2731(6)  0.3722(7)  1.5277(41) 
3.630  0.0170  2000  4.6853(349)  3.1858(167)  0.0773(24)  0.2442(5)  0.3343(9)  1.5293(82) 
3.660  0.0170  2850  4.8497(232)  3.3368(138)  0.0757(7)  0.2364(12)  0.3246(8)  1.5175(50) 
Iii Calculation of the thermodynamic quantities
The calculation of the EoS starts with the evaluation of the trace anomaly, i.e. the trace of the energymomentum tensor . It is related to the temperature derivative of the pressure through thermodynamic identities,
(3) 
The trace anomaly can be expressed in terms of the expectation values of quark condensates and the gluon action density
(4)  
(5)  
(6) 
where the zero temperature expectation values are subtracted to render the trace anomaly UV finite. The expectation values of the light and strange quark condensates and the action are defined as
(7)  
(8) 
with being the staggered fermion matrix. Furthermore, and are the nonperturbative beta function and the mass anomalous dimension
(9)  
(10) 
Following Ref. ourEoS the dependence of the Sommer scale is fitted to a renormalization group inspired Ansatz allton
(11) 
where
(12) 
is the twoloop beta function for 3flavor QCD and . From this the nonperturbative beta function can be calculated as
(13) 
Likewise, was obtained from a parametrization of the bare quark mass,
(14) 
with a sixth order rational function as in ourEoS to account for deviations from the leading order scaling relation.
We performed finite temperature calculations on lattices for the 14 parameter sets shown in Table 1. The number of MD trajectories for each finite temperature run is given in Table 2. Using the expectation values of the quark condensates and the gluonic action density as well as the nonperturbative beta functions described above we have calculated the trace anomaly. The numerical results are shown in Figure 1 and are compared to the previous calculation at twice larger quark mass on lattices ourEoS and lattices hotqcd_eos . The differences between and calculations are due to cutoff effects and have been discussed in Ref. hotqcd_eos .
As one can see from the figure the main differences to the results at arise for temperatures MeV. This difference can be understood as resulting from a shift of the transition temperature when the light quark mass is lowered to approximately its physical value. Based on our previous study on coarser lattices we expect that the reduction of the quark mass from to effectively leads to a shift of several observables by a few MeV towards smaller values of the temperature ourTc . As discussed further in section IV this shift indeed amounts to about MeV. At lower temperatures it also is expected that the trace anomaly increases with decreasing quark masses as hadrons become lighter when the quark mass is decreased. While a tendency for such an increase may be indicated by the data at the lowest two temperatures reached in our calculation, this effect is certainly not significant within the current statistical accuracy. The observed consistency within errors between the trace anomaly calculated at and at might possibly also hint towards a distortion of the hadron spectrum due to finite lattice spacing effects. In fact, it is known that for improved staggered fermions the ground state hadron masses approach their continuum limit from above milc01 ; milc04 . For temperatures MeV there is no visible dependence of on the quark mass. This is presumably due to the fact that hadronic degrees of freedom are no longer relevant in this temperature domain; the relevant degrees of freedom have thermal masses of the order of the temperature and are insensitive to the light quark masses already for .
At temperatures below the transition temperature it is expected that thermodynamic quantities are well described by a hadron resonance gas (HRG) model. In fact, the freezeout of hadrons in heavy ion experiments takes place in the transition region and the observed particle abundances are well described by the HRG model cleymans ; pbm . Therefore in Figure 1 we also show the prediction of the HRG model, which includes all the known resonances up to the mass GeV. The lattice data for are below the HRG prediction although the deviations from it are smaller compared to the results obtained at . We mention again the present statistical accuracy and the possibility of discretization effects in the hadron spectrum. In particular, due to taste breaking of staggered fermions pseudoscalar mesons are not degenerate at finite lattice spacing, therefore their contribution to thermodynamic quantities maybe suppressed.
From the trace anomaly the pressure and thus other thermodynamic quantities can be calculated by performing the integration over the temperature
(15) 
Here is an arbitrary temperature value that is usually chosen in the low temperature regime where the pressure and other thermodynamical quantities are suppressed exponentially by Boltzmann factors associated with the lightest hadronic states, i.e. the pions. Energy and entropy () densities are then obtained by combining results for and .
To perform the integration numerically a reliable interpolation of the lattice data on the trace anomaly is needed. These interpolations are shown in Figure 1 as curves. In the region MeV, the curves correspond to exponential fits. Therefore our parametrization ensures that the pressure is zero at . However, the fits give very small values of the trace anomaly already at temperatures around MeV, therefore we could also use MeV as lower integration limit. Above MeV, we divide the data into several intervals and perform quadratic interpolations. In each interval, these quadratic fits have been adjusted to match the value and slope at the boundary with the previous interval. These interpolating curves are then used to calculate the pressure and other thermodynamic quantities using Eqs. (3) and (15). The numerical results for the pressure and energy density are shown in Figure 2. At temperatures between 170 MeV and 200 MeV, energy density and pressure are slightly larger than in earlier calculations performed at hotqcd_eos . For temperatures below MeV the pressure is almost the same.
Because the discretization errors increase as the temperature decreases it is interesting to consider other choices for the normalization point . At sufficiently low temperatures the HRG model provides a fair estimate for the pressure. In particular, at MeV the pressure is not very sensitive on how many resonances above GeV are included in the HRG model. Therefore we also calculated the pressure taking the lower integration limit to be MeV with calculated in the HRG model. The difference in the pressure and energy density calculated this way and using the standard procedure, where , can be used as an estimate of the systematic error. It is shown in Figure 2 as the horizontal band in upper right corner. The thickness of the band indicates the expected size of the systematic error. The systematic error estimated this way is about in the energy density at the highest temperature of MeV, and about for the pressure.
Iv Deconfinement and chiral aspects of the QCD transition
In this section we are going to discuss the deconfinement and chiral aspects of the QCD transition in terms of renormalized Polyakov loop, subtracted chiral condensate and strangeness susceptibility. The QCD transition in terms of these quantities has been studied previously in Refs. ourEoS ; hotqcd_eos ; fodorTc .
In the previous section we have seen that the energy density shows a rapid rise in the temperature interval MeV. This is usually interpreted to be due to deconfinement, i.e. liberation of many new degrees of freedom. For sufficiently large quark mass this transition is known to be a first order transition (see e.g. Ref. stickan ). In the limit of infinitely large quark mass the order parameter for the deconfinement phase transition is the Polyakov loop. After renormalization it can be related to the free energy of a static quark antiquark pair at infinite separation mclerran81 ; okacz02 ; plqcd
(16) 
A rapid change in this quantity is indicative for deconfinement also in the presence of light quarks. The renormalized Polyakov loop is calculated from the bare Polyakov loop by multiplying it by the renormalization constant given in Table 1,
(17) 
In the above formula denotes the temporal link variables. Note that after performing the ensemble average is independent of the space coordinate .
In the opposite limit of zero quark mass one expects a chiral transition and the corresponding order parameter is the quark condensate defined in section III. For a genuine phase transition, i.e. in the chiral limit the quark condensate vanishes at the critical temperature . However, we expect that even for the crossover at finite quark mass the light quark condensate rapidly drops in the transition region, indicating an approximate restoration of the chiral symmetry. At nonvanishing quark mass the quark condensate needs additive and multiplicative renormalization. Therefore, following Ref. ourEoS we introduce the socalled subtracted chiral condensate
(18) 
Subtraction of the strange quark condensate multiplied by the ratio of the light to strange quark mass removes the quadratic divergence proportional to the quark mass.
In Figure 3 we show the renormalized Polyakov loop and the subtracted chiral condensate and compare with previous calculations performed at light quark masses equal to one tenth of the strange quark mass hotqcd_eos . The renormalized Polyakov loop rises in the temperature interval MeV where we also see the rapid increase of the energy density. At the same time the subtracted chiral condensate rapidly drops in the transition region, indicating that the approximate restoration of the chiral symmetry happens in the same temperature interval as deconfinement. Compared to the calculation performed at light quark masses equal to one tenth of the strange quark mass we see a shift of the transition region by roughly MeV. We note that such a shift arises differently in different observables. In the case of the subtracted chiral condensate, for instance, a major ingredient to the ’shift’ is the fact, that at fixed temperature the condensate in the transition region is strongly quark mass dependent and drops proportional to Goldstone .
The fluctuation of strangeness is also indicative of deconfinement. It can be defined as the second derivative of the free energy density with respect to the strange quark chemical potential
(19) 
At low temperatures strangeness is carried by massive hadrons and therefore strangeness fluctuations are suppressed. At high temperatures strangeness is carried by quarks and the effect of the strange quark mass is small. Therefore strangeness fluctuations are not suppressed at high temperatures. As discussed in Ref. hotqcd_eos strangeness fluctuations behave like the energy density in the transition region, i.e. they rapidly rise in a narrow temperature interval. In Fig. 4 we show the strangeness fluctuations calculated at and compare them with previous calculations performed at hotqcd_eos . In the bottom figure we also show the strangeness fluctuation for with a MeV shift of the temperature scale. As one can see this shift accounts for most of the quark mass dependence of the strangeness fluctuations. This is consistent with the conclusion obtained from the quark mass dependence of other thermodynamic observables.
V Conclusion
We have calculated the EoS, renormalized Polyakov loop, subtracted chiral condensate and strangeness fluctuations in (2+1)flavor QCD in the crossover region from low to high temperatures using the improved p4 staggered fermion formulation on lattices with temporal extent at physical values of the light and strange quark masses. We found that thermodynamic quantities below the deconfinement transition are larger compared to the previous calculations performed at twice larger quark mass but fall below the resonance gas model result. The differences in the thermodynamic quantities calculated at and can be well understood in terms of the shift of the transition temperatures towards smaller values when the quark mass is decreased. This conclusion is also supported by the calculation of renormalized Polyakov loop, subtracted chiral condensate and strangeness fluctuations. No additional enhancement of the pressure and the energy density is seen at low temperatures. This and the deviation from the resonance gas model may be a cutoff effect due to taste violations. However, better statistical accuracy and calculations at smaller lattice spacing are needed to quantify this assertion. At temperatures above MeV no quark mass dependence is seen in the equation of state.
Acknowledgments
This work has been supported in part by contracts DEAC0298CH10886 and DEFG0292ER40699 with the U.S. Department of Energy, the Bundesministerium für Bildung und Forschung under grant 06BI401, the Gesellschaft für Schwerionenforschung under grant BILAER, the Helmholtz Alliance HA216/EMMI and the Deutsche Forschungsgemeinschaft under grant GRK 881. Numerical simulations have been performed on the QCDOC computer of the RIKENBNL research center, the DOE funded QCDOC at BNL, the apeNEXT at Bielefeld University and the BlueGene/L at the New York Center for Computational Sciences (NYCCS).
Appendix
In this appendix we present some numerical details of our calculations. In Table 2 we give the expectation values of the gauge action and quark condensates calculated on (zero temperature) and (finite temperature) lattices. In Table 3 we present the numerical values of the trace anomaly, pressure, energy density, bare Polyakov loop and strangeness fluctuations.
[MeV]  # traj.  

3.4300  0.0370  139  7950  4.10978(19)  4.10922(15)  0.06996(19)  0.06788(12)  0.14253(12)  0.1419(8) 
3.4600  0.0313  154  15290  4.04326(16)  4.04242(17)  0.05263(14)  0.04932(9)  0.11670(8)  0.1157(8) 
3.4900  0.0290  170  33710  3.98357(19)  3.98215(11)  0.04018(13)  0.03453(11)  0.10023(10)  0.0984(7) 
3.5000  0.0253  175  32520  3.96333(14)  3.96161(12)  0.03538(12)  0.02783(18)  0.08873(9)  0.0862(9) 
3.5100  0.0260  180  30050  3.94573(19)  3.94322(7)  0.03294(13)  0.02339(20)  0.08733(10)  0.0839(7) 
3.5200  0.0240  185  39990  3.92713(9)  3.92431(10)  0.02956(6)  0.01718(17)  0.08003(5)  0.0756(8) 
3.5300  0.0240  191  70230  3.90994(11)  3.90651(6)  0.02715(9)  0.01250(16)  0.07735(7)  0.0721(7) 
3.5400  0.0240  196  56740  3.89333(10)  3.88947(6)  0.02518(13)  0.00929(8)  0.07507(8)  0.0691(5) 
3.5450  0.0215  198  20700  3.88427(6)  3.88048(9)  0.02328(6)  0.00681(12)  0.06848(3)  0.0618(9) 
3.5600  0.0205  206  11660  3.85944(14)  3.85579(7)  0.02039(11)  0.00459(6)  0.06315(8)  0.0559(8) 
3.5850  0.0192  219  17250  3.81987(14)  3.81650(10)  0.01638(8)  0.00297(2)  0.05583(7)  0.0479(8) 
3.6000  0.0192  227  4260  3.79720(7)  3.79408(17)  0.01478(7)  0.00272(2)  0.05359(4)  0.0460(9) 
3.6300  0.0170  243  10810  3.75307(9)  3.75040(9)  0.01141(8)  0.00202(7)  0.04524(6)  0.0378(4) 
3.6600  0.0170  259  45050  3.71067(8)  3.70856(3)  0.00916(11)  0.00184(2)  0.04202(5)  0.0355(2) 
[MeV]  

139  0.89(30)  0.067  1.09(30)  0.00087(6)   
154  1.34(30)  0.161  1.82(30)  0.0015(5)  0.062(4) 
170  2.34(30)  0.317  3.29(30)  0.0028(6)  0.150(9) 
175  2.87(27)  0.393  4.05(27)  0.0036(9)  0.186(10) 
180  4.23(29)  0.495  5.71(29)  0.0046(9)  0.261(10) 
185  4.87(19)  0.621  6.73(39)  0.0057(6)  0.335(6) 
191  5.99(19)  0.768  8.30(19)  0.0069(7)  0.392(10) 
196  6.84(18)  0.929  9.63(18)  0.0082(4)  0.462(5) 
198  6.78(18)  1.013  9.81(18)  0.0089(9)  0.492(7) 
206  6.73(27)  1.274  10.55(27)  0.0104(11)  0.599(11) 
219  6.48(31)  1.703  11.59(31)  0.0135(17)  0.750(10) 
227  6.15(34)  1.926  11.93(34)  0.0146(24)  0.798(8) 
243  5.49(26)  2.296  12.38(26)  0.0172(11)  0.841(7) 
259  4.58(18)  2.602  12.38(18)  0.0199(5)  0.874(3) 
References
 (1) J. Engels, F. Karsch, H. Satz and I. Montvay, Phys. Lett. B 101, 89 (1981) and Nucl. Phys. B 205, 545 (1982).

(2)
J. B. Kogut, H. Matsuoka, M. Stone, H. W. Wyld, S. H. Shenker, J. Shigemitsu and
D. K. Sinclair,
Phys. Rev. Lett. 51, 869 (1983);
T. Celik, J. Engels and H. Satz, Phys. Lett. B 129, 323 (1983).  (3) C. E. DeTar, PoS LAT2008, 001 (2008).
 (4) P. Petreczky, Nucl. Phys. Proc. Suppl. 140, 78 (2005).
 (5) P. Petreczky, Nucl. Phys. A 830, 11C (2009).
 (6) Y. Aoki, Z. Fodor, S. D. Katz and K. K. Szabo, JHEP 0601, 089 (2006).
 (7) C. Bernard et al., Phys. Rev. D 75, 094505 (2007).
 (8) M. Cheng et al., Phys. Rev. D 77, 014511 (2008).
 (9) A. Bazavov et al., Phys. Rev. D 80, 014504 (2009).
 (10) M. Cheng et al., Eur. Phys. J. C 51, 875 (2007).
 (11) M. Cheng et al., Phys. Rev. D 74, 054507 (2006).
 (12) A. Gray et al., Phys. Rev. D 72, 094507 (2005).
 (13) C. Allton, Nucl. Phys. B [Proc. Suppl.] 53, 867 (1997).
 (14) C. W. Bernard et al., Phys. Rev. D 64, 054506 (2001).
 (15) C. Aubin et al., Phys. Rev. D 70, 094505 (2004).
 (16) J. Cleymans and K. Redlich, Phys. Rev. C 60, 054908 (1999).
 (17) A. Andronic, P. BraunMunzinger and J. Stachel, Nucl. Phys. A 772, 167 (2006).
 (18) Y. Aoki, Z. Fodor, S. D. Katz and K. K. Szabo, Phys. Lett. B 643, 46 (2006); Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg and K. K. Szabo, JHEP 0906, 088 (2009).
 (19) F. Karsch, C. Schmidt and S. Stickan, Comput. Phys. Commun. 147, 451 (2002).
 (20) L.D.McLerran and B.Svetitsky, Phys. Rev. D 24, 450 (1981).
 (21) O. Kaczmarek, F. Karsch, P. Petreczky and F. Zantow, Phys. Lett. B 543, 41 (2002)
 (22) P. Petreczky and K. Petrov, Phys. Rev. D 70, 054503 (2004); O. Kaczmarek and F. Zantow, Phys. Rev. D 71, 114510 (2005).
 (23) S. Ejiri et al., Phys. Rev. D 80, 094505 (2009).