Finite temperature properties of QCD with two flavors and three, four and five colors
aa r X i v : . [ h e p - l a t ] F e b Finite temperature properties of QCD with two flavors and 3, 4and 5 colors
Thomas DeGrand Department of Physics, University of Colorado, Boulder, CO 80309 USA ∗ (Dated: February 3, 2021) Abstract
I present a numerical study of the crossover between the low temperature chirally broken phaseand the high temperature chirally restored phase in SU ( N c ) gauge theory with N c = 3 − N f = 2 degenerate fermion flavors. Fermion masses span a range of intermediate values(represented by the squared ratio of pseudoscalar to vector meson masses ( m P S /m V ) ∼ .
25 to0.63). Observables include the temperature dependent chiral condensate and screening masses. Ateach fermion mass these quantities show nearly identical temperature dependence across N c . ∗ Electronic address: [email protected] . INTRODUCTION AND MOTIVATION The limit of QCD when the number of colors N c is taken large has a long history as asource of insight about N c = 3 QCD itself [1–3] . Many of the applications of large N c ideasto phenomenology actually involve nonperturbative quantities (masses or decay constants),although the predictions are often based on counting the color weight of Feynman diagrams.To what extent are these predictions true? Checking them requires a lattice simulation, andthere is a small literature of lattice calculations away from N c = 3 to provide such tests.(See Refs. [4–6] for a selection of reviews.)At a qualitative level, lattice calculations confirm large N c intuition rather nicely: mesonmasses, baryon masses, and decay constants scale as N pc f ( m q ) where p is a characteristicpower and m q is the fermion mass. The interplay of small dynamical fermion mass and large N c is less well explored, but simple matrix elements (decay constants, the kaon weak matrixelement calculations of Ref. [7]) mostly scale as expected. So do some chiral observableswhich are governed by the pseudoscalar decay constant (scaling as m P S /f P S ∝ m P S /N c ) orby the condensate (as in the topological susceptibility χ T ∝ m q Σ ∝ m q N c [8]).The subject of this paper is the crossover temperature for the transition from the lowtemperature confining and chirally broken phase to the high temperature deconfined andchirally restored phase, for QCD with N c = 3, 4 and 5 colors and N f = 2 flavors of degeneratemass fermions in the fundamental representation. The reason why this project might beinteresting is that there is not a single large N c story for what happens, there are at leastthree possibilities.The first possibility comes from the naive large N c expectation that gluonic degreesof freedom dominate fermionic ones as N c → ∞ . There should be a finite temperatureconfinement - deconfinement transition which converges to the pure glue one in the large N c limit. This is a first order transition with T c ∼
320 MeV. In N c = 3 the pure gauge transitionis first order and as the fermion mass falls from infinity the transition becomes a crossover.In this scenario the large N c transition temperature would remain roughly constant across N c as the fermion mass fell from infinity, and the critical point where the first order regionends would move to smaller fermion mass as N c rises.The second scenario assumes naive chiral symmetry breaking dominance. Even QCDat large N c has an SU ( N f ) × SU ( N f ) symmetry which (modulo issues with the eta-prime[9]) undergoes spontaneous symmetry breaking to vectorial SU ( N f ). The Pisarski -Wilczekanalysis [10] approximates the Goldstone sector as a linear sigma model. For N f = 2 thesystem is expected to have a second order transition at zero fermion mass, with O (4) criticalexponents. Second order transitions are unstable under perturbation, so the transitionbecomes a crossover away from m q = 0. All QCDs with any N c should then share a commonbehavior at small fermion mass. But where is the crossover temperature T c ? In a linearsigma model, shouldn’t it scale like the pseudoscalar decay constant (which itself scales as √ N c ) [11]? Is T c ∝ √ N c ?The final scenario predates QCD. Confining theories are expected to show an exponen-tially growing spectrum of resonances with mass, forming a Hagedorn spectrum [12]. Thetower of resonances implies a limiting temperature T and (as first stated explicitly byCabibbo and Parisi [13], as far as I can tell) this implies a crossover temperature T c ∼ T . InQCD, the Hagedorn temperature is about 160 MeV. The extension of the story for large N c and nonzero N f is that the spectrum of resonances is basically identical across N c . Mesonstates dominate baryon ones up to a few GeV. If two theories have the same spectrum, then2hey will have the same critical properties. Notice that large N c with nonzero N f is differentfrom quenched QCD: the latter case is has a much sparser spectrum below 2 GeV. Thereare only glueballs in contrast to all the excited states labeled (for example) by quark modelcounting. The prediction of the third scenario is that any N c = 3 with N f = 2 will quali-tatively resemble N c = 3, N f = 2. This argument has been used to justify the observationthat the deconfinement transition in pure gauge SU ( N c ) is nearly independent of N c [14, 15],since the glueball spectrum is also nearly N c independent.The first scenario was already unlikely given that the deconfinement critical point for N c = 3 is already at a very high mass. Ref [16] observes it at a pseudoscalar mass ofabout 4 GeV. And indeed, there is no evidence in any of the simulations reported herefor anything other than a smooth crossover. (This makes the identification of a particularcrossover temperature problematic.) As to the other scenarios, my results indicate that thetemperature dependence of observables computed at common values of the fermion massshow essentially identical behavior, including inflection points at finite temperature, whichis nearly independent of N c .To illustrate this statement, the temperature - dependent and zero - temperature sub-tracted condensate, defined below in Eq. 12, is shown in Fig. 1. The naive N c scaling of thecondensate is divided out, and data is presented for three values of the ratio ( m P S /m V ) .This quantity is close to zero at low temperature and becomes negative at high temperatureas (speaking loosely) chiral symmetry is restored and the finite temperature condensate fallsto zero. The different plotting symbols label different numbers of colors. This picture illus-trates the smooth crossover from broken to restored chiral symmetry with nearly identicaltemperature variation across N c .A natural question to ask is, can one identify a feature in the data which serves as amarker for a crossover temperature? The short answer is No. The quantities I have studiedhave a smooth temperature dependence with no sharp features. I was able to identify twoquantities which could serve as markers: there is a peak in the derivative d Σ( T ) /dT , andthere is a transition region where screening masses cross over from near independence withtemperature to strong temperature dependence. Both features are broad.The outline of the rest of the paper is as follows: Technical aspects of the calculation(lattice action, data sets, a bit about data analysis methodology) are described in Sec. II.This is all completely conventional. Then the various observables are described and re-sults are presented for them: the temperature dependent condensate in Sec. III, screeningmasses in Sec. IV. I mention the the Polyakov line susceptibility in Sec. V. Conclusions aresummarized in Sec. VI. II. TECHNICAL ASPECTS OF THE CALCULATIONA. Simulation methodology, lattice actions, data sets
The simulations involved two degenerate flavors of Wilson–clover fermions. The gaugeaction is the usual Wilson plaquette action. The fermion action uses gauge connectionsdefined as normalized hypercubic (nHYP) smeared links [17–20]. The bare gauge coupling g is set by the simulation parameter β = 2 N c /g . The bare quark mass m q is introducedvia the hopping parameter κ = (2 m q a + 8) − . The clover coefficient is fixed to its tree levelvalue, c SW = 1. Gauge-field updates used the Hybrid Monte Carlo (HMC) algorithm [21–23]with a multi-level Omelyan integrator [24] and multiple integration time steps [25] with one3 IG. 1: The temperature-dependent condensate, rescaled by 3 /N c , as a function of temperature, inappropriate units of t . Squares, octagons, and diamonds label N c = 3, 4, and 5. (a) ( m P S /m V ) ∼ .
63; (b) ( m P S /m V ) ∼ .
5; (c) ( m P S /m V ) ∼ . level of mass preconditioning for the fermions [26]. Lattices used for analysis are spaced aminimum of 10 HMC time units apart. All data sets are based on a single stream for eachset of bare parameters.The only new feature to report in these simulations is a first order bulk transition instrong coupling for SU (5). Its location depends, of course, on the particular form of thebare action. These transitions are well known features of the pure gauge systems and theone I found has a small κ limit which appears to approach the value of gauge coupling whereRef. [14] observed a pure gauge transition. I only mention this in passing; data sets used todo physics were selected to avoid this transition.The choice of Wilson-clover fermions in a finite temperature study is not optimal, sincethe best signals for finite temperature behavior with dynamical fermions are ones which aresensitive to chiral symmetry. The reason I chose this discretization is that I already had a4arge collection of zero temperature simulations which could be used to find lines of constantphysics in bare parameter space, and because I was only looking for gross finite temperaturefeatures. B. Fixing the lattice spacing
The lattice spacing is set by the Wilson flow parameter t [27, 28], and quantities will bepresented as dimensionless ones by a rescaling by an appropriate power of t .The determination of t is done on zero temperature lattices, in the standard way, fromthe observable E ( t ) extracted from the field strength tensor, t h E ( t ) i = C ( N c ) . (1) C ( N c ) is chosen to match what most other large N c simulations take, C ( N c ) = C (cid:18) N c − N c (cid:19) , (2)with C = 0 . SU (3). The extraction of t from the data is identicalto what was done in Ref. [8] and details may be found there.I use values of t computed at each bare parameter value, that is, there is no extrapolationor interpolation in fermion mass. Let’s recall a few useful numbers which will place results incontext. The critical temperature for pure gauge systems was published in Refs. [14, 15]. Theauthors of these papers quote T c / √ σ where σ is the string tension. I converted their numbersto √ t T c using the Sommer parameter r ∼ .
49 fm [29], r √ σ = 1 .
175 for quenched SU (3)and SU (5) from Ref. [30], and the quenched N c = 3 √ t = 0 . √ t T c = 0 . N c = 3, 4, 5.The nominal “physical point” value (with 2 + 1 flavors) of fermions with T c = 150 MeV, √ t = 0 .
147 fm of Ref. [34] from chiral observables, is √ t T c = 0 . C. Data sets
At each N c the bare parameter space is (at least) three dimensional, involving a baregauge coupling, a bare fermion mass expressed through κ in units of the lattice spacing as am q , and the temperature, T = 1 / ( aN t ) where N t is the length of the lattice in the Euclideantemporal direction. The issue for a study like this is that no single simulation (at any oneset of bare parameters) is interesting in itself: all features are broad. It is hard to avoidgenerating data sets at many parameters and easy to get lost wading through them.The literature suggests at least two ways to proceed.The first approach is to map out a crossover line (or phase boundary, if it exists) as afunction of bare parameters at fixed N t values. Then zero temperature simulations are donealong the crossover lines to determine physical quantities such as the lattice spacing and( m P S /m V ) . (Examples of this approach I found useful were Refs. [35–37].) The issue withthis approach is that since there is no true phase transition, the crossover region is broad,and to present a result such as √ t T c versus ( m P S /m V ) involves collecting a large numberof zero temperature data sets followed by interpolation of their output onto a surface in bare5oupling constant space. The surface’s own bare parameters are poorly defined because the“transition” is just a crossover. This issue is compounded of course by the need to work atseveral values of N c . The preliminary version of this project, Ref. [38] used this approach,but I found it to be unwieldy.Instead, I took an approach inspired by Refs. [39, 40]. I started with zero temperaturedata at selected bare parameters and then varied the temperature by varying N t holdingall other bare parameters (and hence, the lattice spacing, t and ( m P S /m V ) ) fixed. Adisadvantage of this approach is that it is unlikely that one will obtain data precisely at acrossover temperature. This can be compensated for, to some extent, by combining datasets from several values of the lattice spacing (but the same ( m P S /m V ) values) to fill in thecurves. This of course introduces another disadvantage: the data sets from different valuesof the lattice spacing have different lattice artifacts. The naive approach of just plottingthem all together ignores lattice artifacts: they will appear in the scatter between points ona graph. An advantage of the approach is that it allows one to use a “natural” observable– the temperature -dependent condensate – to compare crossover behavior across N c . Atthe end, I have eight bare parameter sets per N c (divided into three different ( m P S /m V ) values), five values of N t per bare parameter, – 40 sets per N c , 120 data sets in all.For zero temperature data sets I took 16 ×
32 volumes. Results for many of the parametersets have been published before (see Refs.[8, 20]) and additional sets were generated to giveseveral lattice spacings at three matching points in the fermion mass. The squared ratio ofthe pseudoscalar to vector masses ( m P S /m V ) will be taken as a proxy for a fermion mass.and data was collected at r = ( m P S /m V ) ∼ .
63, 0.5 and about 0.25 (in the latter case,ranging from 0.22-0.27). Useful results from the zero temperature data sets are summarizedin Table I.The calculation of a temperature dependent condensate described below in Sec. III re-quires data sets at fixed spatial volume (for each lattice spacing), again 16 sites. I varied N t from 4 to 12 at r = 0 .
63 and r = 0 .
5, and 6 to 16 at r = 0 .
25 since the crossovertemperature appeared to fall with decreasing r . Almost all of these finite temperature datasets are 1600 trajectories long after equilibration, again saving lattices spaced 10 trajectoriesapart for further analysis. The results for screening masses and related quantities are basedon subsets of these data sets since these results were not used for any precise calculations. D. Data analysis
Most results are global observables, a single quantity averaged over the simulation run.For most observables, there is not really enough data for a full autocorrelation analysis;instead, errors are estimated from a jackknife analysis. For an observable Q I compute anaverage h Q i and a susceptibility and χ Q = h Q i − h Q i ; the uncertainty of each comes froma jackknife. I varied the size of the jackknife (dropping n J successive measurements). Theuncertainty in this quantity has a fractional error from loss of statistics,∆(∆ Q )∆ Q = r n (3)where n = N/n J for N measurements. Plots of h ∆ Q i show a rise as n J increases followedby saturation within the uncertainty of h ∆ Q i when n J reaches the autocorrelation time.Generally, I observed that a jackknife dropping one or two successive lattices (spaced10 molecular dynamics time steps apart) was sufficient to remove time autocorrelations in6bservables related to the temperature dependent condensate. Volume - averaged Polyakovlines were measured every trajectory in the normal course of data collection. There, Eq. 3indicated that typically jackknifes with n J = 5 −
20 (the higher number coming closer tothe crossover temperature) sufficed to handle time autocorrelations.
III. THE TEMPERATURE DEPENDENT CONDENSATE AND RELATEDQUANTITIES
With Wilson fermions the chiral condensate is a bit awkward to measure [39, 41]. Thebare condensate for Wilson fermions at bare mass m has an expansion (cid:10) ¯ ψ ψ (cid:11) = c + c ( m − m ) + c ( m − m ) + . . . (4)where m is the bare mass at which the axial Ward identity fermion mass (defined below)vanishes, m is the bare mass at the simulation point, and the c i ’s contain cutoff (diver-gent) behavior, c ∼ /a , c ∼ /a , c ∼ /a . The divergent pieces are independent oftemperature, and so the difference (cid:10) ¯ ψψ (cid:11) sub = (cid:10) ¯ ψψ (cid:11) T − (cid:10) ¯ ψψ (cid:11) T =0 (5)is finite and sensible. Inspired by the Gell-Mann, Oakes, Renner relation, a potential defi-nition for (cid:10) ¯ ψψ (cid:11) sub is (cid:10) ¯ ψψ (cid:11) sub ∝ m AW I [ Z d x h | P ( x, t ) P (0 , | i T − Z d x h | P ( x, t ) P (0 , | i T =0 ] (6)where P ( x, t ) = ¯ ψ ( x, t ) γ ψ ( x, t ) is the pseudoscalar current. The first term on the righthand side of Eq. 6 is evaluated on an N s × N t lattice where T = 1 / ( aN T ) and the secondterm is evaluated on an N s × N t lattice where N t ≫ N s .With the “130 MeV” definition of f P S the Gell-Mann, Oakes, Renner relation betweenthe condensate Σ and other observables isΣ = m P S f P S m q . (7)In this convention, with the axial current A aµ = ¯ ψγ µ γ ( τ a / ψ , and the pseudoscalar density P a = ¯ ψγ ( τ a / ψ , the vacuum to pseudoscalar matrix elements are h | A | P S i = m P S f P S and h | P | P S i = m P S f P S / (2 m q ) . The partial conservation of axial current relation is ∂ µ A µ ( x, t ) = 2 m q P ( x, t ) . (8)Matrix elements of this relation define m q to be the Axial Ward Identity (AWI) fermionmass, through the two-point functions ∂ t X x h A a ( x , t ) O a i = 2 m q X x h P a ( x , t ) O a i , (9)where O a can be any convenient source. 7he massive correlator in finite volume with periodic temporal boundary conditions intemporal length N t , saturated by a single state of mass m P S , is C ( t ) = X x h P ( x, t ) P (0 , i = | h | P | P S i | cosh( m P S ( N t / − t ))2 m P S sinh( m P S N t / . (10)Thus its integral over the simulation volume is Z N t dtC ( t ) = | h | P | P S i | m P S = ( m P S f P S m q ) m P S = Σ m q . (11)In the passage from lattice to continuum regularization there is a factor of Z A on the righthand side of Eq. 11 where Z A = (1 − (3 κ ) / (4 κ c )) z A is the tadpole improved Z − factor; z A = 1 + cα π ∼ z A is close to unity for the lattice action used here and so we omit itfrom further discussion.The authors of Refs. [39, 40] wrote in the days before the use of t and so they presentedplots of ( m q (cid:10) ¯ ψψ (cid:11) sub ) /m P S as a dimensionless observable. I instead will look at the quantity3 N c t / Σ( T ) = 3 N c t / × m q (∆ P P ( T ) − ∆ P P ( T = 0)) (12)where (explicitly showing the conversion from the lattice quantity computed with cloverfermions to a continuum one) ∆ P P ( T ) = ˆ∆ P P ( N t )(1 − κ κ c ) . (13)and the lattice quantity measured with the usual convention for the definition of lattice fieldvariables is ˆ∆ P P ( N t ) = N t X t =0 X x h P ( x, t ) P (0 , i . (14)The factor of t in Eq. 12 renders the observable dimensionless and the overall factor of3 /N c is included so that plots can show collapse to a common curve across N c when thecondensate scales proportional to N c as expected by large N c counting.Data for the integrated pseudoscalar correlator is recorded in Tables II, III and IV.Fig. 1 shows the rescaled temperature dependent condensate from Eq. 12 at three values of r = ( m P S /m V ) = 0 .
63, 0.5 and about 0.25. These plots are sufficient to show that the finitetemperature behavior of the systems is reasonably independent of N c . However, the curvesare too featureless to identify an inflection point as a signal for a crossover temperature.One expectation would be that the curves in Fig. 1 would show a sigmoid behavior withΣ( T ) zero at low temperature since (loosely speaking) the condensate is unchanged as thetemperature rises, then a fall as the finite temperature condensate goes to zero, followed bya plateau at a constant value, basically the negative of the zero temperature condensate.One could identify a crossover temperature as the midpoint on the sigmoid.The issue with doing this is that T = 1 / ( aN t ) so going to high temperature at fixed a means going to smaller N t , and at some point N t is so small that lattice artifacts mustappear. An alternative is to do a direct measurement of the condensate at T = 0 andcompare it to Fig. 1. There are modern ways based on Ref. [41], but maybe a quicker8 IG. 2: Condensate from the GMOR relation versus ( m P S /m V ) . (a) SU (3), (b) SU (4), (c) SU (5). The plotting symbols label different β values. For SU (3), squares, octagons, diamonds andcrosses are β = 5 .
4, 5.3, 5.45 and 5.25. For SU (4), squares, octagons, diamonds and crosses are β = 10 .
2, 10.3, 10.1 and 10.0. For SU (5), squares, octagons, diamonds and crosses are β = 16 . (though dirtier) way is to use the Gell Mann, Oakes, Renner relation, Eq. 7, using a singleelimination jackknife from separate fits to the AWI quark mass, the decay constant, and thepseudoscalar mass. This will allow for a couple of checks: first, does the zero temperaturecondensate lie on a scaling curve as seen for the subset of the data used in Ref. [20]? Second,will it let us bracket the crossover region in temperature? Fig. 2 shows Σ( m ) plotted as afunction of ( m P S /m V ) . The data exhibit reasonable scaling properties with lattice spacing.Comparing Figs. 1 and 2, the zero temperature condensate appears (rather noisily) to be N c t / Σ ∼ . − .
04 ar r = 0 .
63, 0.02 at r = 0 . r = 0 .
25. The results inFig. 1 seem to be plausible.It would be better to have an observable with a peak, so I looked at two more quantities9
IG. 3: ∆Σ( T ) / ∆ T , rescaled by 3 /N c , as a function of temperature, in appropriate units of t .Squares, octagons, and diamonds label N c = 3, 4, and 5. (a) ( m P S /m V ) ∼ .
63; (b) ( m P S /m V ) ∼ .
5; (c) ( m P S /m V ) ∼ . related to the condensate. The first is the temperature derivative of Σ( T ), just taken fromthe finite difference ∆Σ(( T m )∆ T = Σ( T ) − Σ( T ) T − T (15)where T m = ( T + T ) /
2. The difference in Eq. 15 can just be taken from the integratedpseudoscalar correlator at each value of N t , without doing the T = 0 subtraction. Of course,only data sets at the same bare parameters can be used. Fig. 3 shows this (rescaled by N c and the appropriate power of t ). The differences are (for example for r = 0 . N t = 6 − √ t T ∼ .
18, pretty much independent of( m P S /m V ) . Tis corresponds to a temperature of about 225 MeV, an intermediate value be-10ween the quenched and physical transition points. The figure certainly shows no differencebetween N c = 3, 4 and 5.Attempts to refine this statement came to naught. They were mostly based on doing fitsto an arbitrary peaked function, a Gaussian, y ( x ) = C exp( −
12 ( x − x ) σ ) (16)with x = √ t T . Fits to all data sets with the same N c value (for each choice of ( m P S /m V ) )generally had poor chi-squared, probably due to lattice artifacts: the data span a wide rangeof t values. Fits to a single set of bare parameters fared better, though the issue is thatthere are only four values of N t . Some of the data sets (especially the ones at coarse latticespacing) do not themselves include the peak and then of course the fit fails immediately.I also attempted to compute a susceptibility from the time histories of ∆ P P ( T ). This wasunsuccessful and I just mention it in passing. It is just χ P P = (cid:10) (∆ P P ( T )) (cid:11) − h ∆ P P ( T ) i (17)with an error taken from a jackknife. To convert to a dimensionless parameter, I rescale χ in lattice units by the square of the factor used for Σ( T ):( 3 N c ) t χ = [ 3 N c m q t / ] χ P P (18)This result is shown in Fig. 4.(I believe that this is not the same quantity as the “disconnected chiral susceptibility”measured in large scale SU (3) simulations. There, the condensate is measured directly fromΣ f = TV ∂∂m f ln Z (19)and a susceptibility could come from one more derivative, χ ∝ TV ∂∂m f Σ . (20)Eq. 19 says that the condensate is the one point function Σ ∝ h ¯ qq i = (cid:10) Tr D − ii (cid:11) where D isthe Dirac operator and χ is a difference of (cid:10) (Tr D − ii ) (cid:11) and (cid:10) Tr D − ii (cid:11) .)Fig. 4 is obviously noisy. The issue is (I think) that in the simulation runs with largeuncertainties, the time history of ˆ∆ P P ( N t ).shows occasional spikes above a flat background,reminiscent of the “exceptional configurations” of olden times. These are associated withfluctuations in the correlator C ( t ) when it flattens across t . One could speculate as to thecause (the log normal distribution of C ( t ) [42, 43]?) but for now I will just move on. IV. SCREENING MASSES
Meson screening masses in the scalar, pseudoscalar, vector, and axial vector channelsare taken from two-point correlation functions extending in a spatial lattice direction. Thetemperature - dependent pseudoscalar decay constant is also extracted from these spatial11
IG. 4: The susceptibility of the temperature dependent condensate, rescaled by (3 /N c ) , as afunction of temperature, in appropriate units of t . Squares, octagons, and diamonds label N c = 3,4, and 5. (a) ( m P S /m V ) ∼ .
63; (b) ( m P S /m V ) ∼ .
5; (c) ( m P S /m V ) ∼ . correlators. Propagators are constructed with composite boundary conditions to double theeffective length of the lattice [44–46].In the low temperature, chirally broken phase, the spectroscopy of screening massesshould qualitatively resemble ordinary T = 0 spectroscopy with a light pion and no degen-eracies in the spectrum. When chiral symmetry is restored parity partners (the pion andthe scalar mesons, the vector and axial vector mesons) should become degenerate, and allfour states should become degenerate when U (1) A is restored. A naive expectation for ascreening mass is that it behaves something like m H = 4 "(cid:18) πN t (cid:19) + m q (21)where π/N t is the lowest nonzero Matsubara frequency associated with antiperiodic bound-12 IG. 5: Pseudoscalar and scalar meson screening masses ( √ t m P S and √ t m S ) versus √ t T . (a) m P S /m V = 0 .
63, (b) m P S /m V = 0 . m P S /m V = 0 .
25. Symbols are squares for SU (3),octagons for SU (4) and diamonds for SU (5) for the pseudoscalar and fancy squares, fancy crossesand fancy diamonds for the scalar, for SU (3), SU (4) and SU (5), respectively. The line is just m H = 2 πT . ary conditions in a lattice of temporal length N t . Since 1 /N t = aT , this gives m H = 2 πT athigh temperature.Results for screening masses (scaled by √ t ) are displayed in Fig. 5 for the pseudoscalarand scalar states and Fig. 6 for the vector an axial vector states, They show the expectedbehavior. (The diagonal lines are just y = 2 πT , the small m q (or large T ) limit of Eq. 21.)The very dirty signals for the a (scalar) and a (axial vector) mesons in the chirally brokenphase are also expected. There, the pseudoscalar meson is light and the noise to signal ratio σ ( t ) /C ( t ) ∼ exp(( m H − m π )) t is exponentially bad [47, 48]; in the chirally restored phaseall states are massive and the noise to signal ratio improves.13 IG. 6: Vector and axial vector meson screening mass ( √ t m V and √ t m A ) versus √ t T . (a) m P S /m V = 0 .
63, (b) m P S /m V = 0 . m P S /m V = 0 .
25. Symbols are squares for SU (3),octagons for SU (4) and diamonds for SU (5) for the vector and fancy squares, fancy crosses andfancy diamonds for the axial vector, for SU (3), SU (4) and SU (5). The line is just m H = 2 πT . The pseudoscalar and vector screening masses make a clear transition from a temperatureindependent value at low T to linear behavior at high T . This suggests that linear fits O = c + c T should show better quality (smaller chi-squared) when fits keeping only datain one phase are included, and c would be much larger in the higher T phase. The fit woulddeteriorate when points in both phases were included. This would determine a crossovertemperature.The transition can be seen in individual data sets, and no fit is needed to see it, just atable of the mass values. In most sets the change in slope occurs in the middle of the dataset, with (typically) two or three masses with nearly the same value at low temperatureand the remaining higher temperature masses rising linearly with temperature. Data set14y data set, one can identify a √ t T low and a √ t T high where masses remain unchangedfor T ≤ T low and rise linearly for √ t T ≥ √ t T high . Generally, √ t T low and √ t T high areconsistent between the pseudoscalar and vector mass sets. The individual data sets give arelatively wide interval between √ t T low and √ t T high simply because T = 1 / ( aN t ) and thevalues of N t are small.Fits to all mass values at each N c and ( m P S /m V ) value also show jumps in chi-squaredwhen a fit including high temperature points extends too low, or fits to low temperaturepoints extends too high. The issue is that even with data in one phase, the chi-squaredtended to be unacceptably large due to the a dependent variation in the masses. The bestone can say is that the crossover region is in the range √ t T c ∼ . − .
24, which is notinconsistent with the location of the broad peak in d Σ /dT .The last screening quantity to present is the pseudoscalar decay constant. Here, thefollowing alternative seems to show a sharper result than plots of f P S versus temperatureat fixed r . In the chirally broken phase, f P S is nonvanishing at zero fermion mass (due tospontaneous breaking of chiral symmetry) and rises modestly with fermion mass (due to toexplicit breaking of chiral symmetry from the fermion mass). In the chirally restored phasethere is only explicit symmetry breaking, and f P S should be proportional to m q . A plotof f P S /m q will show collapse to a common value when that situation occurs. Fig. 7 showsthat behavior. I have scaled f P S by p /N c . Here the different plotting symbols in eachpanel represent different values of ( m P S /m V ) with the diamonds representing the lightestfermion mass.Notice that the smaller mass f P S /m q points (diamonds) seem to fall onto the commontemperature line at larger √ t T than the higher mass ones (octagons and squares) do –perhaps a suggestion that the crossover moves to lower temperature as the fermion massfalls. V. POLYAKOV LINE AND ITS SUSCEPTIBILITY
It is well known that the Polyakov line is not a sensitive observable at small fermionmass, but I present one picture for completeness. Following Ref. [15], I define the volumeaveraged Polyakov line as l P = 1 N s N c X x Tr N t − Y t =0 U ( x, t ) (22)and then the susceptibility is χ P = N s ( (cid:10) | l P | (cid:11) − h| l P |i ) . (23)I measure these quantities from time histories of the average Polyakov line taking a jackknifeaverage over the simulation run. I have only measured the original Polyakov line, not anysmoothed one.The Polyakov line susceptibility is shown in Fig. 8. The bare Polyakov line is not a scalingquantity, so the results for different N c values should not coincide. The only comment tomake about these noisy figures is that the magnitude of the susceptibility falls with thefermion mass, as seen for N c = 3 [49].Note that the shoulder in the Polyakov line susceptibility appears at √ t T ∼ .
15 or so,about where the screening masses begin to take their high temperature functional form.15
IG. 7: Pseudoscalar decay constant divided by the AWI fermion mass ( f P S /m q ) versus √ t T . (a) SU (3), (b) SU (4), (c) SU (5). The plotting symbols correspond to ( m P S /m V ) = 0 .
63 for squares,0.5 for octagons, and 0.25 for diamonds.
VI. CONCLUSIONS
The question this study was designed to address was whether the finite temperaturecrossover behavior of SU ( N c ) gauge theories with N f = 2 flavors of fermions showed differentbehavior as the number of colors was varied. These simulations, with low statistics andcarried out on small volumes, studied the temperature dependence of the condensate (asdetermined from the volume integral of the pseudoscalar correlator) and of screening massesand the pseudoscalar decay constant. A smooth crossover from a low temperature confiningand chirally broken phase to a high temperature chirally restored and deconfined phase isobserved. This crossover behavior is of course not surprising given our extensive knowledge16 IG. 8: Polyakov line susceptibility as a function of temperature, scaled by √ t .of t . Squares,octagons, and diamonds label N c = 3, 4, and 5. Lines connect simulation results from the samebare parameters. (a) ( m P S /m V ) ∼ .
63; (b) ( m P S /m V ) ∼ .
5; (c) ( m P S /m V ) ∼ . of N c = 3: the one new result is that the temperature dependence of these quantities (whencompared at fixed values of m P S /m V ) shows no observable dependence on N c .It was not possible to determine a crossover temperature (to the extent that such aquantity makes sense when the crossover is smooth) but, with some plausible assumptions,it appears to be someplace between the known SU (3) result with physical quark masses( √ t T ∼ .
11) and the large N c pure gauge result ( √ t T c ∼ . SU (3) results for T c /m V as a function of m P S /m V . With a nominal crossover temperature √ t T c ∼ .
18 and m V taken from Table I, T c /m V ∼ .
20, 0.23, 0.25 at r = 0 .
63, 0.5, 0.25, inreasonably good agreement with Refs. [35, 36]. Ref. [37] shows a plot of r T c versus r m P S where r is the Sommer parameter [29], r / √ t ∼
3. This conversion gives r T c ∼ .
54 at17 m P S = 2 .
2, 1.6 and 1.1 at r = 0 .
63, 0.5, 0.25, to be compared with r T c = 0 . − . r m P S values.Of the three scenarios described in the Introduction, the first two seem to be disfavored.No evidence was seen for any first order behavior at large mass at any N c studied. Ofcourse, it could be that N c = 5 is still not “large N c ” from the point of view of thermo-dynamics. Scenario two (or at least the scenario that T c ∝ f P S ) s disfavored because thethe crossover temperature shows no obvious N c separation despite the (known) variation ofthe pseudoscalar decay constant with respect to N c . As far as the third scenario, latticesimulations show that the lowest part of the meson spectrum show little N c dependence.The third scenario assumes that the two results I have presented – a common spectrumand a common crossover behavior – are correlated. Of course, the scenario asks that thecorrelation persists high in the meson spectrum (that the density of states is as given byHagedorn) and lattice simulations say nothing about that.The results in this paper are rather poor quality, but at the mass values studied, every-thing is smooth. The situation at very low mass and very high mass remains open. I thinkthat to do more work in this area, one should not continue to use Wilson-clover fermions –at least, not in pilot tests. The importance of chiral symmetry, the need to subtract resultsfrom finite and zero temperature simulations, and the fierce scaling of the cost with simu-lation volume of thermodynamical observables argue in favor of using some modern versionof staggered fermions for dedicated simulations.The obvious next target in future studies of large N c systems at finite temperature couldbe a more direct attack on thermodynamic observables – the internal energy, pressure,speed of sound, and related quantities. The pure gauge calculations of Panero [50] are aninspiration, and they have been extensively cited in the plenomenological and gravitationalduality literature of QCD thermodynamics. But we know that the critical behavior of N c = 3 with dynamical fermions is different from the behavior of pure gauge systems. Doesthis matter, for the questions which interest researchers in these areas? If so, numericalsimulation which includes dynamical fermions might be a worthwhile project. Acknowledgments
Daniel Hackett participated in the earliest stages of this project. I am grateful to RobPisarski for a conversation about large N c expectations for QCD thermodynamics. Mycomputer code is based on the publicly available package of the MILC collaboration [51].The version I use was originally developed by Y. Shamir and B. Svetitsky. This material isbased upon work supported by the U.S. Department of Energy, Office of Science, Office ofHigh Energy Physics under Award Number DE-SC-0010005. Some of the computations forthis work were also carried out with resources provided by the USQCD Collaboration, whichis funded by the Office of Science of the U.S. Department of Energy using the resources of theFermi National Accelerator Laboratory (Fermilab), a U.S. Department of Energy, Office ofScience, HEP User Facility. Fermilab is managed by Fermi Research Alliance, LLC (FRA),acting under Contract No. DE- AC02-07CH11359. [1] G. ’t Hooft, Nucl. Phys. B , 461 (1974). doi:10.1016/0550-3213(74)90154-0 c β κ κ c am q a m P S a m V ( m P S /m V ) t /a N r ∼ .
633 5.3 0.1250 0.12923 0.153 0.757(1) 0.925(5) 0.670(7) 0.871(3) 303 5.4 0.1250 0.12838 0.107 0.563(1) 0.707(1) 0.634(3) 1.657(3) 4003 5.45 0.1250 0.12795 0.093 0.497(2) 0.636(5) 0.611(11) 2.284(21) 304 10.1 0.1240 0.12794 0.150 0.716(2) 0.881(3) 0.661(6) 1.189(7) 304 10.2 0.1245 0.12792 0.108 0.556(1) 0.701(1) 0.629(3) 1.966(3) 1904 10.3 0.1240 0.12759 0.112 0.537(2) 0.657(3) 0.668(8) 2.727(18) 505 16.3 0.1230 0.12795 0.162 0.726(1) 0.875(3) 0.688(5) 1.467(7) 305 16.4 0.1240 0.12785 0.119 0.582(1) 0.725(1) 0.644(3) 2.030(2) 1905 16.6 0.1240 0.12740 0.106 0.514(2) 0.638(2) 0.649(6) 2.942(14) 50 r ∼ .
53 5.25 0.1280 0.12964 0.080 0.545(2) 0.773(4) 0.497(6) 0.863(2) 303 5.4 0.1265 0.12838 0.058 0.395(1) 0.563(2) 0.492(4) 2.019(5) 4003 5.45 0.1265 0.12795 0.044 0.331(2) 0.486(6) 0.464(13) 2.747(17) 914 10.0 0.1270 0.12926 0.104 0.623(3) 0.860(6) 0.525(9) 0.855(2) 404 10.2 0.1262 0.12792 0.054 0.377(1) 0.561(2) 0.452(4) 2.270(4) 1904 10.3 0.1260 0.12759 0.049 0.343(3) 0.491(5) 0.488(13) 3.106(17) 905 16.2 0.1260 0.12853 0.087 0.531(1) 0.757(4) 0.492(6) 1.290(7) 305 16.4 0.1258 0.12785 0.063 0.404(1) 0.592(2) 0.466(4) 2.272(5) 1905 16.6 0.1252 0.12740 0.069 0.406(2) 0.549(3) 0.547(8) 3.108(19) 40 r ∼ .
253 5.4 0.1276 0.12838 0.021 0.234(2) 0.444(8) 0.278(11) 2.412(7) 4003 5.45 0.1273 0.12795 0.018 0.228(4) 0.437(6) 0.272(12) 3.164(23) 904 10.2 0.1272 0.12792 0.022 0.238(2) 0.472(8) 0.254(10) 2.520(20) 1014 10.3 0.1270 0.12759 0.017 0.201(3) 0.432(6) 0.216(9) 3.375(21) 905 16.4 0.1270 0.12785 0.025 0.248(1) 0.493(3) 0.253(4) 2.483(6) 2105 16.6 0.1268 0.12740 0.019 0.208(2) 0.435(14) 0.229(16) 3.565(23) 50TABLE I: Zero temperature data sets from 16 ×
32 volumes. (The SU (5) β = 16 . κ = 0 . × κ c is thevalue of the hopping parameter where the Axial Ward Identity fermion mass vanishes. It is neededto perform tadpole renormalization of condensate - related quantities.[2] G. ’t Hooft, Nucl. Phys. B , 461 (1974). doi:10.1016/0550-3213(74)90088-1[3] E. Witten, Nucl. Phys. B , 57-115 (1979) doi:10.1016/0550-3213(79)90232-3[4] B. Lucini and M. Panero, Phys. Rept. , 93 (2013) doi:10.1016/j.physrep.2013.01.001[arXiv:1210.4997 [hep-th]].[5] M. Garcia Perez, arXiv:2001.10859 [hep-lat].[6] P. Hern´andez and F. Romero-L´opez, [arXiv:2012.03331 [hep-lat]].[7] A. Donini, P. Hern´andez, C. Pena and F. Romero-L´opez, Phys. Rev. D , no.11, 114511(2016) doi:10.1103/PhysRevD.94.114511 [arXiv:1607.03262 [hep-ph]].[8] T. DeGrand, Phys. Rev. D , no.11, 114509 (2020) doi:10.1103/PhysRevD.101.114509[arXiv:2004.09649 [hep-lat]].[9] R. Kaiser and H. Leutwyler, Eur. Phys. J. C , 623 (2000) doi:10.1007/s100520000499 c β κ N t = 4 N t = 6 N t = 8 N t = 12 N t = 323 5.3 0.1250 16.64(2) 18.11(6) 18.34(5) 18.61(12) 18.35(9)3 5.4 0.1250 16.49(2) 17.30(4) 18.06(5) 18.02(8) 18.11(8)3 5.5 0.1250 16.45(2) 17.14(3) 17.61(3) 17.86(5) 17.89(9)4 10.1 0.1240 22.03(2) 23.44(4) 24.08(5) 24.01(9) 23.88(8)4 10.2 0.1245 21.93(2) 22.91(2) 23.72(4) 23.73(6) 23.94(6)4 10.3 0.1240 21.84(1) 22.68(2) 23.24(3) 23.36(4) 23.42(5)5 16.3 0.1230 27.40(2) 28.58(3) 29.35(3) 29.52(8) 29.41(13)5 16.4 0.1240 27.41(2) 28.60(3) 29.54(4) 29.54(5) 29.54(7)5 16.6 0.1240 27.26(1) 28.28(2) 28.82(2) 29.17(4) 29.21(4)TABLE II: The (bare, lattice regulated) quantity ˆ∆ P P ( N t ) as defined in Eq. 14 for r ∼ . × N t . N c β κ N t = 4 N t = 6 N t = 8 N t = 12 N t = 323 5.25 0.1280 16.90(9) 18.94(9) 20.12(16) 19.94(12) 19.96(18)3 5.4 0.1265 16.58(2) 17.43(8) 17.96(5) 18.78(11) 18.72(16)3 5.45 0.1265 16.49(2) 17.26(2) 17.65(3) 18.36(8) 18.45(8)4 10.0 0.1270 22.44(9) 25.88(10) 26.22(9) 26.16(14) 26.00(9)4 10.2 0.1262 22.10(2) 23.08(2) 23.96(5) 25.12(12) 24.68(8)4 10.3 0.1260 21.99(1) 22.92(2) 23.41(3) 24.35(7) 24.34(8)5 16.2 0.1260 27.79(2) 29.87(7) 31.51(7) 31.51(10) 31.77(23)5 16.4 0.1258 27.61(2) 28.82(2) 30.52(5) 30.54(6) 30.66(4)5 16.6 0.1252 27.40(1) 28.50(2) 29.01(2) 29.95(6) 29.79(5)TABLE III: The (bare, lattice regulated) quantity ˆ∆ P P ( N t ) at r ∼ .
5. Lattice volumes are16 × N t .[hep-ph/0007101].[10] R. D. Pisarski and F. Wilczek, Phys. Rev. D , 338 (1984). doi:10.1103/PhysRevD.29.338[11] Thanks to Rob Pisarski for discussions about this point.[12] R. Hagedorn, Nuovo Cim. A , 1027 (1968). doi:10.1007/BF02751614[13] N. Cabibbo and G. Parisi, Phys. Lett. , 67 (1975). doi:10.1016/0370-2693(75)90158-6[14] B. Lucini, M. Teper and U. Wenger, JHEP , 033 (2005) doi:10.1088/1126- N c β κ N t = 6 N t = 8 N t = 12 N t = 16 N t = 323 5.4 0.1276 17.43(3) 18.07(6) 19.87(15) 20.05(18) 20.80(27)3 5.5 0.1273 17.32(2) 17.69(3) 19.08(11) 19.69(21) 19.48(15)4 10.2 0.1272 23.25(2) 23.88(3) 26.69(14) 26.55(17) 27.04(24)4 10.3 0.1270 23.00(2) 23.79(11) 25.60(13) 25.33(13) 25.72(11)5 16.4 0.1270 29.00(3) 30.21(7) 33.10(13) 33.05(17) 32.95(17)5 16.6 0.1268 28.66(2) 29.28(3) 31.59(9) 31.97(23) 32.25(16)TABLE IV: The (bare, lattice regulated) quantity ˆ∆ P P ( N t ) at r ∼ .
25. Lattice volumes are16 × N t . , 279 (2012)doi:10.1016/j.physletb.2012.04.070 [arXiv:1202.6684 [hep-lat]].[16] F. Cuteri, O. Philipsen, A. Sch¨on and A. Sciarra, [arXiv:2009.14033 [hep-lat]].[17] A. Hasenfratz and F. Knechtli, Phys. Rev. D , 034504 (2001).doi:10.1103/PhysRevD.64.034504 [hep-lat/0103029].[18] A. Hasenfratz, R. Hoffmann and S. Schaefer, JHEP , 029 (2007). doi:10.1088/1126-6708/2007/05/029 [hep-lat/0702028].[19] T. DeGrand, Y. Shamir and B. Svetitsky, Phys. Rev. D , 074506 (2012).doi:10.1103/PhysRevD.85.074506 [arXiv:1202.2675 [hep-lat]].[20] T. DeGrand and Y. Liu, Phys. Rev. D , no. 3, 034506 (2016) Erratum: [Phys. Rev. D , no. 1, 019902 (2017)] doi:10.1103/PhysRevD.95.019902, 10.1103/PhysRevD.94.034506[arXiv:1606.01277 [hep-lat]].[21] S. Duane and J. B. Kogut, Nucl. Phys. B , 398 (1986). doi:10.1016/0550-3213(86)90606-1[22] S. Duane and J. B. Kogut, Phys. Rev. Lett. , 2774 (1985). doi:10.1103/PhysRevLett.55.2774[23] S. A. Gottlieb, W. Liu, D. Toussaint, R. L. Renken and R. L. Sugar, Phys. Rev. D , 2531(1987). doi:10.1103/PhysRevD.35.2531[24] T. Takaishi and P. de Forcrand, Phys. Rev. E , 036706 (2006).doi:10.1103/PhysRevE.73.036706 [hep-lat/0505020].[25] C. Urbach, K. Jansen, A. Shindler and U. Wenger, Comput. Phys. Commun. , 87 (2006).doi:10.1016/j.cpc.2005.08.006 [hep-lat/0506011].[26] M. Hasenbusch, Phys. Lett. B , 177 (2001). doi:10.1016/S0370-2693(01)01102-9[hep-lat/0107019].[27] R. Narayanan and H. Neuberger, JHEP , 064 (2006) doi:10.1088/1126-6708/2006/03/064[hep-th/0601210]. M. Luscher, Commun. Math. Phys. , 899 (2010) doi:10.1007/s00220-009-0953-7 [arXiv:0907.5491 [hep-lat]].[28] M. L¨uscher, JHEP , 071 (2010) Erratum: [JHEP , 092 (2014)]doi:10.1007/JHEP08(2010)071, 10.1007/JHEP03(2014)092 [arXiv:1006.4518 [hep-lat]].[29] R. Sommer, Nucl. Phys. B , 839 (1994). doi:10.1016/0550-3213(94)90473-1[hep-lat/9310022].[30] T. DeGrand, Phys. Rev. D , 034508 (2012) doi:10.1103/PhysRevD.86.034508[arXiv:1205.0235 [hep-lat]].[31] S. Lottini [ALPHA Collaboration], PoS LATTICE , 315 (2014) doi:10.22323/1.187.0315[arXiv:1311.3081 [hep-lat]].[32] M. Bruno et al. [ALPHA Collaboration], PoS LATTICE , 321 (2014)doi:10.22323/1.187.0321 [arXiv:1311.5585 [hep-lat]].[33] R. Sommer, PoS LATTICE , 015 (2014) doi:10.22323/1.187.0015 [arXiv:1401.3270 [hep-lat]].[34] R. A. Soltz, C. DeTar, F. Karsch, S. Mukherjee and P. Vranas, Ann. Rev. Nucl. Part. Sci. ,379 (2015) doi:10.1146/annurev-nucl-102014-022157 [arXiv:1502.02296 [hep-lat]].[35] A. Ali Khan et al. [CP-PACS Collaboration], Phys. Rev. D , 074510 (2001)doi:10.1103/PhysRevD.64.074510 [hep-lat/0103028].[36] S. Ejiri et al. [WHOT-QCD Collaboration], Phys. Rev. D , 014508 (2010)doi:10.1103/PhysRevD.82.014508 [arXiv:0909.2121 [hep-lat]].[37] V. G. Bornyakov, R. Horsley, S. M. Morozov, Y. Nakamura, M. I. Polikarpov,P. E. L. Rakow, G. Schierholz and T. Suzuki, Phys. Rev. D , 014504 (2010) oi:10.1103/PhysRevD.82.014504 [arXiv:0910.2392 [hep-lat]].[38] T. DeGrand, D. C. Hackett and E. T. Neil, PoS LATTICE2018 , 175 (2018)doi:10.22323/1.334.0175 [arXiv:1809.00073 [hep-lat]].[39] S. Borsanyi et al. , JHEP , 126 (2012) doi:10.1007/JHEP08(2012)126 [arXiv:1205.0440[hep-lat]].[40] S. Borsanyi et al. , Phys. Rev. D , no. 1, 014505 (2015) doi:10.1103/PhysRevD.92.014505[arXiv:1504.03676 [hep-lat]].[41] L. Giusti, F. Rapuano, M. Talevi and A. Vladikas, Nucl. Phys. B , 249-277 (1999)doi:10.1016/S0550-3213(98)00659-2 [arXiv:hep-lat/9807014 [hep-lat]].[42] M. G. Endres, D. B. Kaplan, J. W. Lee and A. N. Nicholson, PoS LATTICE2011 , 017 (2011)doi:10.22323/1.139.0017 [arXiv:1112.4023 [hep-lat]].[43] T. DeGrand, Phys. Rev. D , 014512 (2012) doi:10.1103/PhysRevD.86.014512[arXiv:1204.4664 [hep-lat]].[44] T. Blum et al. [RBC], Phys. Rev. D , 114506 (2003) doi:10.1103/PhysRevD.68.114506[arXiv:hep-lat/0110075 [hep-lat]].[45] Y. Aoki, T. Blum, N. H. Christ, C. Dawson, T. Izubuchi, R. D. Mawhinney,J. Noaki, S. Ohta, K. Orginos and A. Soni, et al. Phys. Rev. D , 094507 (2006)doi:10.1103/PhysRevD.73.094507 [arXiv:hep-lat/0508011 [hep-lat]].[46] T. DeGrand and S. Schaefer, [arXiv:0712.2914 [hep-lat]].[47] G. P. Lepage, “The Analysis Of Algorithms For Lattice Field Theory,” invited lectures at the1989 TASI summer school, Boulder CO, June 4-30, 1989. CLNS-89-971.[48] G. Parisi, Phys. Rept. , 203-211 (1984) doi:10.1016/0370-1573(84)90081-4[49] D. A. Clarke, O. Kaczmarek, F. Karsch and A. Lahiri, PoS LATTICE2019 , 194 (2020)doi:10.22323/1.363.0194 [arXiv:1911.07668 [hep-lat]].[50] M. Panero, Phys. Rev. Lett. , 232001 (2009) doi:10.1103/PhysRevLett.103.232001[arXiv:0907.3719 [hep-lat]].[51] https://github.com/milc-qcd/milc qcd/https://github.com/milc-qcd/milc qcd/