Distribution of Energy-Momentum Tensor around a Static Quark in the Deconfined Phase of SU(3) Yang-Mills Theory
Ryosuke Yanagihara, Masakiyo Kitazawa, Masayuki Asakawa, Tetsuo Hatsuda
JJ-PARC-TH-0229, RIKEN-iTHEMS-Report-20
Distribution of Energy-Momentum Tensor around a Static Quarkin the Deconfined Phase of SU(3) Yang-Mills Theory
Ryosuke Yanagihara, ∗ Masakiyo Kitazawa,
1, 2, † Masayuki Asakawa, ‡ and Tetsuo Hatsuda § Department of Physics, Osaka University, Toyonaka, Osaka 560-0043, Japan J-PARC Branch, KEK Theory Center, Institute of Particle and Nuclear Studies,KEK, 203-1, Shirakata, Tokai, Ibaraki, 319-1106, Japan RIKEN Interdisciplinary Theoretical and Mathematical Sciences Program (iTHEMS), RIKEN, Wako 351-0198, Japan
Energy momentum tensor (EMT) characterizes the response of the vacuum as well as the thermalmedium under the color electromagnetic fields. We define the EMT by means of the gradient flowformalism and study its spatial distribution around a static quark in the deconfined phase of SU(3)Yang-Mills theory on the lattice. Although no significant difference can be seen between the EMTdistributions in the radial and transverse directions except for the sign, the temporal component issubstantially different from the spatial ones near the critical temperature T c . This is in contrastto the prediction of the leading-order thermal perturbation theory. The lattice data of the EMTdistribution also indicate the thermal screening at long distance and the perturbative behavior atshort distance. I. INTRODUCTION
To study complex quantum systems such as the Yang-Mills (YM) theory, it is customary to introduce testprobe(s) and analyze the response. The Wilson loop isone of such probes whose measurement in YM theoryprovides information on the static quark–anti-quark sys-tem that is closely related to the confinement propertyin YM vacuum [1]. Thanks to the recent developmentof the gradient-flow method [2–4] and its application tothe energy-momentum tensor (EMT) T µν ( x ) [5–8], it be-came possible to study the gauge-invariant structure ofthe flux tube between the quark and anti-quark in theconfining phase through the spatial distribution of EMTunder the Wilson loop [9, 10].The purpose of the present paper is to extend theabove idea and to explore the EMT distribution arounda static quark in YM theory. As a first step, we con-sider the deconfined phase above the critical temperature T c of the SU(3) YM theory in the range of temperature . ≤ T /T c ≤ . and measure the EMT distributionaround the Polyakov loop. The EMT with the gradi-ent flow has been used to study thermodynamics of YMtheory [11–16] and of QCD [17, 18]. However, the ob-servables in these studies are limited to global quantitiessuch as the pressure, energy density, entropy density, andthe specific heat. On the other hand, we focus on the lo-cal observable in this study and examine the followingquestions: (i) How are the energy density and the stresstensor distributed around the static quark?, (ii) How arethe distributions modified as a function of temperature?,and (iii) How can one extract parameters such as therunning coupling and the Debye screening mass from thedistributions? ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] The organization of the present paper is as follows. InSec. II, we briefly review the definition of EMT and itsproperty in the spherical coordinate system. In Sec. III,we introduce the EMT operator on the lattice and itscorrelation with the Polyakov loop operator. In Sec. IV,we discuss the numerical procedure and lattice setup toanalyze the EMT operator around a static quark on thelattice. Numerical results and their physical implicationsare given in Sec. V. Sec. VI is devoted to the summaryand conclusion. In Appendix A, we discuss the procedureto make the tree-level improvement of the correlation be-tween the EMT and the Polyakov loop on the lattice. InAppendix B, the leading-order perturbative analysis ofthe correlation is presented using the high temperatureeffective field theory.
II. EMT AROUND A STATIC CHARGE
The stress tensor σ ij ( i, j = 1 , , is related to thespatial component of EMT, T ij , as [19] σ ij = −T ij . (1)The force per unit area F i acting on a surface with thenormal vector n i is given by the stress tensor as F i = σ ij n j = −T ij n j . (2)The local principal axes n ( k ) j and the corresponding eigen-values λ k of the local stress tensor are obtained by solvingthe eigenvalue problem: T ij n ( k ) j = λ k n ( k ) i ( k = 1 , , . (3)The strength of the force per unit area along n ( k ) i is givenby the absolute values of the eigenvalue λ k . Neighboringvolume elements separated by a surface with the normalvector n ( k ) i pull (push) each other for λ k < ( λ k > )across the wall. Note that the three principal axes n ( k ) i a r X i v : . [ h e p - l a t ] O c t source 𝑟𝜃 FIG. 1. Stress acting on infinitely small volume elementunder the existence of a single static charge (source). In thecase of the classical electromagnetism, a small volume elementat a distance of r from the charge is pulled along the radialdirection from the neighboring volume elements while it ispushed in the transverse direction. are orthogonal with one another because σ ij is a sym-metric tensor.For a system with a single static source, it is convenientto use the spherical coordinate system ( r, θ, ϕ ) with theradial coordinate r = | x | , and the polar and azimuthalangles θ and ϕ . The spherical symmetry allows us to di-agonalize the static EMT in Euclidean spacetime T µν ( x ) in this coordinate system as T γγ (cid:48) ( x ) = diag( T ( r ) , T rr ( r ) , T θθ ( r )) (4)where γ, γ (cid:48) = 4 , r, θ . Due to the spherical symmetry, theazimuthal component degenerates with the polar compo-nent, T ϕϕ ( r ) = T θθ ( r ) , so that only independent compo-nents are given in Eq. (4).In the Abelian case, EMT is given by the Maxwellstress-energy tensor [19], T Maxwell µν = F µρ F νρ − δ µν F ρσ F ρσ with the field strength F µν . When a staticcharge is placed at the origin, the EMT is denoted by T Maxwell γγ (cid:48) = 12 diag( − (cid:126)E , − (cid:126)E , (cid:126)E ) , (5)with E i ( x ) being the electric field. The spatial structureof Eq. (5) is illustrated in Fig. 1 where the neighboringvolume elements around the static electric charge pull(push) each other along the radial (angular) direction.In a static system, the force acting on a volume elementthrough its surface should be balanced. This propertyis guaranteed by the momentum conservation ∂ i T ij = 0 together with the Gauss theorem. III. EMT FOR SU(3) YANG-MILLS THEORYON THE LATTICEA. YM gradient flow
We consider the pure SU(3) YM gauge theory in thefour-dimensional Euclidean space defined by the action, S YM = 14 g (cid:90) d x G aµν ( x ) G aµν ( x ) . (6)Here g is a bare gauge coupling and G aµν ( x ) is the fieldstrength composed of the fundamental gauge field A aµ ( x ) .The YM gradient flow evolves the gauge field along thefictitious fifth dimension t introduced in addition to theordinary four Euclidean dimensions x through the flowequation [2–4], dA aµ ( t, x ) dt = − g δS YM ( t ) δA aµ ( t, x ) . (7)The flowed YM action S YM ( t ) in the (4 + 1) -dimensionalcoordinate is constructed by substituting the flowedgauge field A aµ ( t, x ) in Eq. (6) with an initial condition, A aµ ( t = 0 , x ) = A aµ ( x ) .An important feature of the gradient flow for t > isthat any composite operators composed of flowed gaugefields are UV finite even at equal spacetime point [4, 7].This is a consequence of the smoothing of the gauge fieldsin the four dimensional Euclidean space within the range ∼ √ t . In addition, in the small t limit, composite localoperators are represented by the local operators of theordinary gauge theory at t = 0 . These properties leadus to the renormalized EMT operator defined with thesmall t expansion [5]: T R µν ( x ) = lim t → T µν ( t, x ) , (8) T µν ( t, x ) = c ( t ) U µν ( t, x )+4 c ( t ) δ µν [ E ( t, x ) − (cid:104) E ( t, x ) (cid:105) ] , (9)where (cid:104) E ( t, x ) (cid:105) is the vacuum expectation value of E ( t, x ) . The dimension-four gauge-invariant operatorson the right hand side of Eq. (9) are given by [5] E ( t, x ) = 14 G aµν ( t, x ) G aµν ( t, x ) , (10) U µν ( t, x ) = G aµρ ( t, x ) G aνρ ( t, x ) − δ µν E ( t, x ) , (11)where G aµν ( t, x ) is the field strength composed of theflowed gauge field. Because of the vacuum subtrac-tion in Eq. (9), (cid:104)T R µν ( x ) (cid:105) vanishes. The coefficients c ( t ) and c ( t ) have been calculated perturbatively inRefs. [5, 8, 14] for small t . We use two-loop perturbativecoefficients [8, 14] for the construction of EMT through-out this study. B. EMT around a static heavy quark
To describe a static quark Q on the lattice, we intro-duce the Polyakov loop at the origin, Ω( ) . Then theexpectation value of Eq. (9) around Q is given by (cid:104)T µν ( t, x ) (cid:105) Q = (cid:104)T µν ( t, x )TrΩ( ) (cid:105)(cid:104) TrΩ( ) (cid:105) − (cid:104)T µν ( t, x ) (cid:105) . (12)We note that Eq. (12) is well-defined only when the Z symmetry in SU(3) YM theory is spontaneously broken:In the Z unbroken phase, both numerator and denomi-nator of the first term on the right hand side vanish ex-actly. This is the reason why we focus on the system inthe Z broken phase above T c in this paper. In practice,we choose the state with the Polyakov loop being realamong the three equivalent Z states in the deconfinedphase.The renormalized EMT distribution around Q is ob-tained after taking the double extrapolation, (cid:104)T R µν ( x ) (cid:105) Q = lim t → lim a → (cid:104)T µν ( t, x ) (cid:105) Q . (13)In our actual analysis, we extract the renormalized EMTdistribution by fitting the lattice data with the followingfunctional form [12, 13]: (cid:104)T µν ( t, x ) (cid:105) Q = (cid:104)T R µν ( x ) (cid:105) Q + b µν ( t ) a + c µν t + d µν t , (14)where the contributions from discretization effects ( b µν )as well as the dimension-six and -eight operators ( c µν and d µν ) are considered.To perform the double extrapolation reliably, thesmearing radius ρ ≡ √ t needs to be larger than thelattice spacing to suppress the discretization error. Atthe same time, ρ should be smaller than half the tempo-ral size / T with temperature T as well as the distancefrom the source ( r ) to avoid the overlap of operators.Therefore we require a/ (cid:46) ρ (cid:46) min (cid:16) r, T (cid:17) . (15)The lattice data to be fitted by Eq. (14) should be withinthis window. As will be discussed in Sec. IV C, we im-pose more stringent conditions for the range of t in ournumerical analysis. IV. LATTICE SETUPA. Gauge configurations
Numerical simulations in SU(3) YM theory wereperformed on the four dimensional Euclidean lat-tice with the Wilson gauge action and the periodicboundary conditions at four different temperatures . T c , . T c , . T c , and . T c . The simulation pa-rameters for each T are summarized in Table I. The in-verse coupling β = 6 /g is related to the lattice spacing a determined by the reference scale w [12, 20]. Thespatial and temporal lattice sizes, N s and N τ together TABLE I. Simulation parameters for the four temperatures:The spatial lattice size N s , the temporal lattice size N τ , β =6 /g , the lattice spacing a . N conf represents the number ofconfigurations. T /T c N s N τ β a [fm] N conf with the number of configurations N conf are also summa-rized in Table I. All lattices have the same aspect ratio N s /N τ = 4 .The gauge configurations are generated by the pseudo-heat-bath method followed by five over-relaxations. Eachmeasurement is separated by 200 sweeps. Statisticalerrors are estimated by the jackknife method with 20jackknife bins. We employ the Wilson gauge action for S YM ( t ) in the flow equation Eq. (7) and the clover typerepresentation for the field strength G µν ( t, x ) . The nu-merical solution of the gradient flow equation is obtainedby the third order Runge-Kutta method.In order to suppress the statistical noise, we apply themulti-hit procedure in the measurement of the Polyakovloop by replacing every temporal link by its thermal aver-age [21]. The choice of the temporal argument x of EMTin Eq. (12) is arbitrary. Therefore, we average EMT overthe temporal direction to reduce the statistical error. B. Discretization effect
The EMT in the spherical coordinate system on thelattice reads T (cid:88) x (cid:104)T γγ (cid:48) ( t, x , x ) (cid:105) Q = diag( (cid:104)T ( t, r ) (cid:105) Q , (cid:104)T rr ( t, r ) (cid:105) Q , (cid:104)T θθ ( t, r ) (cid:105) Q ) . (16)The behavior of the EMT distribution close to the source Q is affected by the violation of rotational symmetry ow-ing to lattice discretization. As an example, we showin Fig. 2 the distribution of − r (cid:104)T rr ( t/a . = 1 . , r ) (cid:105) Q as a function of rT at T /T c = 1 . , where a . is thelattice spacing of the finest lattice at this temperature.The figure shows the oscillating behavior of the numer-ical results becomes more prominent on coarser lattices.In this study, we use the lattice data only for N τ ≥ forthe continuum extrapolation to suppress the discretiza-tion errors. In Appendix A, we consider an alternativeanalysis that performs the tree-level improvement of thenumerical results and uses them for the continuum ex-trapolation with the N τ = 10 data. As discussed there,we confirm that the results in both cases are the samewithin the errors. C. Double extrapolation
The double extrapolation Eq. (13) consists of twosteps: (I) the continuum ( a → ) extrapolation, and (II) t → extrapolation. In this subsection, we demonstratethese procedures by using the lattice data at T /T c = 1 . as an example.In Fig. 3, we show −(cid:104)T rr ( t, rT ) (cid:105) Q /T at rT =0 . , . , . as a function of /N τ = ( aT ) for fourvalues of tT . To obtain the EMT values at a given r and t , we first perform the linear interpolation of the lat-tice data along the r direction and then interpolate alongthe t direction by the cubic spline method for each a .In Fig. 3, fitting results of the data at N τ = 12 – ac-cording to Eq. (14) at fixed t are shown by the solid lines,while the results of the continuum limit are shown by thefilled squares on the vertical dotted line at /N τ = 0 . Wethen take t → extrapolation by fitting the continuum-extrapolated results for different t with Eq. (14) at a = 0 .This fit has to be carried out within the range of t satis-fying Eq. (15). We employ tT = 0 . correspondingto t/a = 1 . of the finest lattice data as the lower boundof the fitting window: This choice satisfies Eq. (15) forall the lattices. The upper bound of the fitting windowis taken to be tT = 0 . , since the thermodynamic rT − r › T rr ( t / a . = . , r ) fi Q T/T c = 1 . N τ =10 N τ =12 N τ =14 N τ =16 N τ =18 FIG. 2. Distribution of − r (cid:104)T rr ( t/a . = 1 . , r ) (cid:105) Q as func-tions of rT at T /T c = 1 . . /N τ − › T rr ( t , r T = . ) fi Q / T (a) tT =0 . tT =0 . tT =0 . tT =0 . / o N t =10 /N τ − › T rr ( t , r T = . ) fi Q / T (b) tT =0 . tT =0 . tT =0 . tT =0 . / o N t =10 /N τ − › T rr ( t , r T = . ) fi Q / T (c) tT =0 . tT =0 . tT =0 . tT =0 . / o N t =10 FIG. 3. Color open symbols represent −(cid:104)T rr ( t, rT ) (cid:105) Q /T for various t as functions of /N τ = a T at T /T c = 1 . .The continuum extrapolation for each t is shown by the solidlines, with the extrapolated results represented by the filledsymbols. Panels (a), (b) and (c) show the results at rT =0 . , . , . , respectively. quantities show a linear behavior below this value [12].We consider the following three ranges within . ≤ tT ≤ . to estimate the systematic un-certainty from the fitting ranges [12]: Range-1 : 0 . ≤ tT ≤ . , Range-2 : 0 . ≤ tT ≤ . , Range-3 : 0 . ≤ tT ≤ . . Range-1 is the most conservative window, while Range-2(Range-3) is the extension of Range-1 towards the smaller(larger) values of t . We employ the result of Range-1 asa central value and use the Range-2 and Range-3 for an − › T ( t , r T = . ) fi Q / T (a1) Range - - - tT − › T ( t , r T = . ) fi Q / T (a2) continuum N τ =18 N τ =16 N τ =14 N τ =12Range - - - − › T rr ( t , r T = . ) fi Q / T (b1) tT − › T rr ( t , r T = . ) fi Q / T (b2) continuum N τ =18 N τ =16 N τ =14 N τ =12Range - - - › T θθ ( t , r T = . ) fi Q / T (c1) tT › T θθ ( t , r T = . ) fi Q / T (c2) continuum N τ =18 N τ =16 N τ =14 N τ =12Range - - - FIG. 4. Each component of EMT at rT = 0 . (top) and rT = 0 . (bottom) as functions of tT . Color open symbols denote (cid:104)T γγ (cid:48) ( t, rT ) (cid:105) Q /T for each a as functions of tT . The black solid line with the gray error band is the continuum-extrapolatedresult. The dotted lines show the fitted results of the continuum result with Range-1, 2, and 3. The black symbols are theresults of the t → extrapolation for these fitting ranges. estimate of the systematic error. In the following, all theresults after the double extrapolation contain both thestatistical and systematic errors.In Fig. 4, the open symbols with statistical errors rep-resent (cid:104)T γγ (cid:48) ( t, rT ) (cid:105) Q /T at rT = 0 . (upper) and . (lower) for each a . The results of the continuum limit aredenoted by the black solid lines with the gray statisticalerror band for . ≤ tT ≤ . . Range-1 is high-lighted by the yellow band. The figure also shows thefitted results for Ranges-1, 2, and 3 by the dotted lines.The final results of the t → limit for each range areshown by the open black symbols around tT = 0 : Theyagree with each other within the statistical errors, whichsuggests that the systematic uncertainty from the choiceof the fitting range is not significant. V. RESULTS OF EMT DISTRIBUTIONS
Before entering into the detailed discussions on the spa-tial distribution of EMT, we first show the result of thestress distribution at
T /T c = 2 . on a two-dimensionalplane including the static source in Fig. 5. The same re-sult is later shown in Fig. 8 in a different form. In Fig. 5,the red and blue arrows represent the principal directionsof the stress tensor along the radial and transverse direc-tions, respectively. The length of each arrow representssquare root of the eigenvalue corresponding to each prin-cipal axis. This figure is to be compared with Fig. 1 inRef. [9]. x [fm] y [ f m ] λ k < λ k > FIG. 5. Stress distribution around a static quark at theorigin at
T /T c = 2 . . The red and blue arrows are theprincipal directions along the radial and transverse directions,respectively. The length of each arrow represents square rootof the eigenvalue. A. Channel dependence
In Fig. 6, we show the dimensionless EMT, −(cid:104)T
R44 ( r ) (cid:105) Q /T , −(cid:104)T R rr ( r ) (cid:105) Q /T , and (cid:104)T R θθ ( r ) (cid:105) Q /T , asfunctions of the dimensionless length rT . The errorbands include both the statistical and systematic errors,where the latter is estimated from the three fitting rangesfor the t → extrapolation. Since the thermal expecta-tion value (cid:104)T µν ( t, x ) (cid:105) is subtracted as in Eq. (12), we have (cid:104)T R µν ( r ) (cid:105) Q → in the r → ∞ limit.We find that −(cid:104)T R44 ( r ) (cid:105) Q , −(cid:104)T R rr ( r ) (cid:105) Q , and (cid:104)T R θθ ( r ) (cid:105) Q rT (a) T/T c = 1 . − › T R44 ( r ) fi Q /T − › T R rr ( r ) fi Q /T › T R θθ ( r ) fi Q /T rT (b) T/T c = 1 . − › T R44 ( r ) fi Q /T − › T R rr ( r ) fi Q /T › T R θθ ( r ) fi Q /T rT (c) T/T c = 2 . − › T R44 ( r ) fi Q /T − › T R rr ( r ) fi Q /T › T R θθ ( r ) fi Q /T rT (d) T/T c = 2 . − › T R44 ( r ) fi Q /T − › T R rr ( r ) fi Q /T › T R θθ ( r ) fi Q /T FIG. 6. EMT distribution ( −(cid:104)T R44 ( r ) (cid:105) Q , −(cid:104)T R rr ( r ) (cid:105) Q , (cid:104)T R θθ ( r ) (cid:105) Q ) as functions of rT after the double extrapolation: (a) T =1 . T c , (b) T = 1 . T c , (c) T = 2 . T c , and (d) T = 2 . T c . are all positive for rT (cid:46) and decrease rapidly withincreasing r . These signs are the same as those of theMaxwell stress tensor in Eq. (5). Individual signs phys-ically mean that a volume element has a positive local-ized energy density and receives a pulling (pushing) forcealong the longitudinal (transverse) direction; see Figs. 1and 5.Figure 6 indicates that the absolute values of the spa-tial components |(cid:104)T R rr ( r ) (cid:105) Q | and |(cid:104)T R θθ ( r ) (cid:105) Q | are degener-ated within the error for all temperatures. On the otherhand, |(cid:104)T R44 ( r ) (cid:105) Q | is larger than the spatial componentsespecially at lower temperature. This is in contrast to thedegenerate magnitude of all components in the Maxwellstress Eq. (5) and is also different from the leading-orderthermal perturbation theory (Appendix B). B. Temperature dependence
Shown in Fig. 7 is the temperature dependence of thespatial distribution of EMT with respect to the physicaldistance r [fm]; (a) −(cid:104)T R44 ( r ) (cid:105) Q , (b) −(cid:104)T R rr ( r ) (cid:105) Q , and (c) (cid:104)T R θθ ( r ) (cid:105) Q . Also, shown in Fig. 7(d) is the distribution ofthe trace of EMT given by ∆ Q ( r ) ≡ −(cid:104)T R µµ ( r ) (cid:105) Q = −(cid:104)T R44 ( r ) + T R rr ( r ) + 2 T R θθ ( r ) (cid:105) Q . (17)Figure 7 tells us that the EMT distributions have small T dependence at short distances, r (cid:46) . fm. On the other hand, for large distances, sizable T dependence can beseen despite the growth of the errors at high T .To make these features more explicit, we plot the sameresults with a dimensionless normalization r (cid:104)T R γγ (cid:48) ( r ) (cid:105) Q as a function of r in Fig. 8. The figure shows that the T dependence is suppressed for rT (cid:46) . and all resultsapproach a single line, while the result tends to be moresuppressed for rT (cid:38) . compared with this universal be-havior as temperature is raised. This result is reasonableas the T dependence of r (cid:104)T R γγ (cid:48) ( r ) (cid:105) Q would be suppressedfor r (cid:46) (2 πT ) − .At distance (2 πT ) − (cid:28) r (cid:28) Λ − with Λ being thelambda parameter, the behavior of (cid:104)T R γγ (cid:48) ( r ) (cid:105) Q shouldbe described by the perturbation theory in electrostaticQCD (EQCD). In the leading order of EQCD in thisregime, we have the following ratio (Appendix B), (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∆ Q ( r ) (cid:104)T R44 ,rr,θθ ( r ) (cid:105) Q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 112 π α s + O ( g ) , (18)which is independent of r and T and is given only by afunction of α s .Shown in Fig. 9 is the r dependence of ∆ Q ( r ) / (cid:104)T R44 ,rr,θθ ( r ) (cid:105) Q as a function of r at T /T c = 2 . .From this result and Eq. (18), we obtain, at r = 0 . fm, α s = 0 . from − ∆ Q ( r ) / (cid:104)T R44 ( r ) (cid:105) Q , α s = 0 . from − ∆ Q ( r ) / (cid:104)T R rr ( r ) (cid:105) Q , and α s = 0 . from ∆ Q ( r ) / (cid:104)T R θθ ( r ) (cid:105) Q . Although these values are channeldependent, indicating the existence of non-negligiblehigher order contributions, it is notable that they areconsistent with that obtained from the similar analysis r [fm] − › T R ( r ) fi Q [ G e V / f m ] (a) . T c . T c . T c . T c r [fm] − › T R rr ( r ) fi Q [ G e V / f m ] (b) . T c . T c . T c . T c r [fm] › T R θθ ( r ) fi Q [ G e V / f m ] (c) . T c . T c . T c . T c r [fm] ∆ Q ( r ) [ G e V / f m ] (d) . T c . T c . T c . T c FIG. 7. EMT distribution ( −(cid:104)T R44 ( r ) (cid:105) Q , −(cid:104)T R rr ( r ) (cid:105) Q , (cid:104)T R θθ ( r ) (cid:105) Q ) and ∆ Q ( r ) as func-tions of r [fm] at each temperature: (a) −(cid:104)T R44 ( r ) (cid:105) Q , (b) −(cid:104)T R rr ( r ) (cid:105) Q , (c) (cid:104)T R θθ ( r ) (cid:105) Q , and (d) ∆ Q ( r ) . of the Polyakov loop correlations at r = 0 . fm [22, 23].Higher order α s corrections and the thermal correctionsfor EMT around a static charge to be compared withour lattice data is under way [24].Let us now turn to the long-distance region in Fig. 8.Owing to the large errors in this region, it is not possibleto extract the thermal screening of the form exp( − m D r ) with m D being the Debye screening mass. Nevertheless,Figs. 8(a-d) indicate that the EMT distributions decrease r [fm] − r › T R ( r ) fi Q (a) . T c . T c . T c . T c r [fm] − r › T R rr ( r ) fi Q (b) . T c . T c . T c . T c r [fm] r › T R θθ ( r ) fi Q (c) . T c . T c . T c . T c r [fm] r ∆ Q ( r ) (d) . T c . T c . T c . T c FIG. 8. EMT distribution r ( −(cid:104)T R44 ( r ) (cid:105) Q , −(cid:104)T R rr ( r ) (cid:105) Q , (cid:104)T R θθ ( r ) (cid:105) Q ) and r ∆ Q ( r ) asfunctions of r [fm] at each temperature: (a) − r (cid:104)T R44 ( r ) (cid:105) Q ,(b) − r (cid:104)T R rr ( r ) (cid:105) Q , (c) r (cid:104)T R θθ ( r ) (cid:105) Q , and (d) r ∆ Q ( r ) . faster than /r at long distances, and the tendency isstronger at high temperatures. To draw a definite con-clusion, however, higher statistical data are necessary. r [fm] T/T c = 2 . − ∆ Q ( r ) / › T R44 ( r ) fi Q − ∆ Q ( r ) / › T R rr ( r ) fi Q ∆ Q ( r ) / › T R θθ ( r ) fi Q FIG. 9. Ratios ∆ Q ( r ) / (cid:104)T R44 ,rr,θθ ( r ) (cid:105) Q as functions of r at T /T c = 2 . . VI. SUMMARY AND CONCLUDINGREMARKS
In the present paper, we have studied, for the firsttime, the EMT distribution around a static quark at fi-nite temperature above T c of the SU(3) YM theory onthe lattice. The YM gradient flow plays crucial roles todefine the EMT on the lattice and to explore its spatialstructure.The main results of this paper can be summarized asfollows.As shown in Fig. 6, we found no significant differencebetween the absolute magnitude of EMT along the ra-dial direction and that of the transverse direction for alltemperatures above T c . This seems to be in accordancewith the leading-order thermal perturbation theory inQCD, which predicts the same magnitude for all princi-pal components of EMT. However, we found a substantialdifference between the EMT distribution in the tempo-ral direction and that of the spatial directions, especiallynear T c . This indicates that there is indeed a genuinenon-Abelian effect present at finite temperature, so thatprecise comparison with the higher-order thermal QCDcalculation would be called for.As shown in Figs. 7 and 8, all the EMT distributionshave small T dependence at short distances, r (cid:46) . fm.Also the EMT distributions decrease faster than /r atlong distances, and the tendency is stronger at high tem-peratures. However, owing to the large statistical errors,we could not extract the values of the thermal Debyescreening. By using the fact that the EMT distribu-tions are T independent at short distances, we attemptedto extract the strong coupling constant from the ratiosbetween the different components of EMT. The result, α s (0 . (cid:39) . − . , is consistent with that obtainedfrom the similar analysis for the Q ¯ Q free energy at finite T .We have some important issues to be studied further:Going beyond the leading-order thermal QCD calcula-tion for the EMT [24] is necessary to understand thelattice results presented in this paper. At the same time, increasing the statistics of lattice data is necessary to ex-tract, e.g. the screening mass from the long range partof the EMT distribution.There are also several interesting future problems.First of all, the extension to full QCD is an importantnext step. Since the Z symmetry is explicitly broken bydynamical fermions, the present method can be applieddirectly to a single static quark Q , a static diquark QQ and Q ¯ Q both at low and high temperatures. In particu-lar, the single quark system in QCD at zero temperaturecorresponds to a heavy-light meson [25]. Secondly, theEMT distributions of the QQQ system will provide newinsight into the flux tube formation in baryons [26] aswell as the “gravitational” baryon structure [27–30] atzero and non-zero temperatures.
ACKNOWLEDGMENT
The authors thank T. Iritani for fruitful discussionsin the early stage of this study. They are also grate-ful to M. Berwein for discussions regarding the pertur-bative analysis of the EMT distribution. M. K. thanksF. Karsch for useful discussions. The numerical simula-tion was carried out on OCTOPUS at the CybermediaCenter, Osaka University and Reedbush-U at Informa-tion Technology Center, The University of Tokyo. Thiswork was supported by JSPS Grant-in-Aid for ScientificResearches, 17K05442, 18H03712, 18H05236, 18K03646,19H05598, 20H01903.
Appendix A: Tree-level improvement of the latticeobservables
As shown in Fig. 2, there exists sizable discretizationeffect for the EMT distribution on coarse lattices espe-cially for N τ = 10 . In this Appendix, we attempt toreduce such discretization effects by using the tree-levellattice propagator. Similar idea has been applied to theanalysis of the Polyakov loop correlations in Refs. [22, 31].In calculating the EMT distribution around a staticquark Eq. (16), we need the expectation values ofEqs. (10) and (11) at nonzero flow time t . These quan-tities are constructed from (cid:104) G aµν ( t, x ) G aρσ ( t, x ) (cid:105) Q , wherethe temporal coordinate is suppressed for notational sim-plicity. In the continuum theory at the tree level, thisoperator is calculated to be (cid:104) G aµν ( t, x ) G aρσ ( t, x ) (cid:105) Q = − g N − G µν ( t, x ) G ρσ ( t, x ) , (A1)where g is the gauge coupling and N = 3 is the numberof colors, with G µν ( t, x ) δ ab = (cid:90) /T dτ (cid:104) A a ( τ, ) G bµν ( t, x ) (cid:105) . (A2) rT − r › T rr ( t / a . = . , r ) fi Q (a) w /N t =10w / o N t =10 N τ =10 N τ =12 N τ =14 N τ =16 N τ =18 rT − r › T rr ( t / a . = . , r ) fi Q (b) w /N t =10w / o N t =10 N τ =10 N τ =12 N τ =14 N τ =16 N τ =18 rT − r › T rr ( t / a . = . , r ) fi Q (c) w /N t =10w / o N t =10 N τ =10 N τ =12 N τ =14 N τ =16 N τ =18 rT − r › T rr ( t / a . = . , r ) fi Q (d) w /N t =10w / o N t =10 N τ =10 N τ =12 N τ =14 N τ =16 N τ =18 FIG. 10. Distribution of − r (cid:104)T rr ( t/a . = 1 . , r ) (cid:105) Q asfunctions of rT . (a) Same figure as Fig. 2 shown for a com-parison. (b,c,d) Tree-level improved results with Eqs. (A7)and (A10) for three choices of U tl γγ (cid:48) : (b) U tl γγ (cid:48) = (cid:104) U (cid:105) Q , (c) U tl γγ (cid:48) = (cid:104) U rr (cid:105) Q , (d) U tl γγ (cid:48) = (cid:104) U θθ (cid:105) Q . By selecting appropriate gauge fixing conditions for thegauge action and the gradient flow equation, one obtains G ( t, x ) = G ( t, x ) = G ( t, x ) = 0 and G i ( t, x ) = ∂ i D ( t, x ) , (A3) D ( t, x ) = (cid:90) d p (2 π ) e i p · x e − tp p = 14 π | x | erf (cid:16) | x |√ t (cid:17) . (A4)Next, in lattice gauge theory the propagator corre-sponding to Eq. (A4) with the Wilson gauge action for the gauge action and the flow equation reads [32, 33] D ( t, x n ) = (cid:90) π − π d p (2 π ) e i p · x n e − t (cid:80) i ˆ p i (cid:80) i ˆ p i , (A5)with x n = a n = a ( n x , n y , n z ) and ˆ p i = (2 /a ) sin( p i / a ) .When the clover-leaf operator for the discretized repre-sentation of G aµν ( t, x n ) is employed, the discretized rep-resentation of G i ( t, x n ) is given by [34] G lat i ( t, x n ) = 12 a (cid:0) D ( t, x n +ˆ i ) − D ( t, x n − ˆ i ) (cid:1) . (A6)Using Eq. (A3) and (A6), the tree-level improvementsof Eqs. (10) and (11) denoted by the superscript ‘imp’may be written as (cid:104) E ( t, x n ) (cid:105) imp Q = c ( t, x n ) (cid:104) E ( t, x n ) (cid:105) Q , (A7) (cid:104) U γγ (cid:48) ( t, x n ) (cid:105) imp Q = c ( t, x n ) (cid:104) U γγ (cid:48) ( t, x n ) (cid:105) Q , (A8)where the correction factor c ( t, x n ) is defined by c ( t, x n ) = 13 (cid:88) i =1 (cid:16) G i ( t, x n ) G lat i ( t, x n ) (cid:17) . (A9)In Eq. (A9), the average over i is taken because gener-ally the ratio G i ( t, x n ) / G lat i ( t, x n ) at a lattice site x n depends on i . However, in our particular choice of dis-cretization, i.e. the Wilson gauge actions and the clover-leaf operator, it is easily shown that G i ( t, x n ) / G lat i ( t, x n ) does not depend on i . In this special case the averageover i in Eq. (A9) is redundant. When this property isviolated, the improvement of (cid:104) U γγ (cid:48) ( t, x n ) (cid:105) Q may be re-placed by the one that depends on γ and γ (cid:48) in place ofEq. (A8).There is one more subtle issue about Eq. (A8). Atthe tree level, one easily finds that the matrix elementsof U , U rr and − U θθ are the same, while the actual lat-tice data do not necessarily satisfy such relation as shownin the main text. In our tree-level improvement, there-fore, we decompose our lattice data into tree-like partand the rest, (cid:104) U γγ (cid:48) ( t, x n ) (cid:105) Q = U tl γγ (cid:48) ( t, x n ) + δU γγ (cid:48) ( t, x n ) ,where the tree-like part U tl γγ (cid:48) ( t, x n ) satisfies U tl44 ( t, x n ) = U tl rr ( t, x n ) = − U tl θθ ( t, x n ) . Then we apply our tree-levelimprovement only to the first term: (cid:104) U γγ (cid:48) ( t, x n ) (cid:105) imp Q = c ( t, x n ) U tl γγ (cid:48) ( t, x n ) + δU γγ (cid:48) ( t, x n ) . (A10)Below, we consider three choices of U tl γγ (cid:48) ( t, x n ) to es-timate the systematic uncertainty of this procedure: U tl γγ (cid:48) = (cid:104) U (cid:105) Q , (cid:104) U rr (cid:105) Q , and −(cid:104) U θθ (cid:105) Q . Corresponding re-sults for (cid:104)T R rr ( t, r ) (cid:105) as example are shown in Fig. 10(b),(c), and (d), respectively, together with the case withoutthe correction, Fig. 10(a). Colored open symbols repre-sent the data at each N τ and the gray (yellow) shade isthe continuum result with (without) the data at N τ = 10 .0 /N τ − › T rr ( t T = . , r T = . ) fi Q / T (a) default U tl γγ = › U fi Q U tl γγ = › U rr fi Q U tl γγ = − › U θθ fi Q /N τ − › T rr ( t T = . , r T = . ) fi Q / T (b) default U tl γγ = › U fi Q U tl γγ = › U rr fi Q U tl γγ = − › U θθ fi Q /N τ − › T rr ( t T = . , r T = . ) fi Q / T (c) default U tl γγ = › U fi Q U tl γγ = › U rr fi Q U tl γγ = − › U θθ fi Q /N τ − › T rr ( t T = . , r T = . ) fi Q / T (d) default U tl γγ = › U fi Q U tl γγ = › U rr fi Q U tl γγ = − › U θθ fi Q FIG. 11. Comparison of the different prescriptions inthe continuum extrapolation of −(cid:104)T rr ( t, rT = 0 . (cid:105) Q /T at T /T c = 1 . : (a) tT = 0 . , (b) tT = 0 . , (c) tT = 0 . , and (d) tT = 0 . . The figures show that the tree-level improvement sup-presses the discretization effect at short distances in threecases, especially for N τ = 10 .Let us now compare the continuum extrapolation us-ing the data including N τ = 10 with the tree-level im-provement and that using the data without N τ = 10 andwithout the tree-level improvement. Shown in Fig. 11 FIG. 12. Diagram contributing to the leading-order calcula-tion of the correlation function between the Polyakov loop andthe EMT operator. The vertical line represents the Polyakovloop and two wavy lines exchanged gluons. The symbol ⊗ corresponds to the EMT operator. is such a comparison for −(cid:104)T rr ( t, rT = 0 . (cid:105) Q /T at T /T c = 1 . as functions of /N τ . Colored open trian-gles represent the data without the tree-level improve-ment (red) and with the tree-level improvement for threedifferent prescriptions (blue, orange, and green). Filledsymbols at /N τ = 0 are the continuum extrapolation:The red squares are continuum results without N τ = 10 data discussed in the main text, while the diamonds arethe continuum results with N τ = 10 data after the tree-level improvement. Taking into the uncertainly asso-ciated with the different prescriptions for the tree-levelimprovement, the default results without N τ = 10 in themain text are found to be consistent with the improvedresults including N τ = 10 . Appendix B: Leading order perturbative analysis ofEMT around a static charge
Let us consider the SU( N ) Yang-Mills system at hightemperature where g (2 πT ) (cid:28) . Then, the effective the-ory valid at the length scale of R (cid:29) (2 πT ) − is the di-mensionally reduced electrostatic QCD (EQCD) in threedimensions (see, e.g., [35–38]) S EQCD = (cid:90) d x (cid:104)
12 Tr G + Tr( Dϕ ) + m D Tr ϕ + δ L EQCD (cid:105) . (B1)Here ( A i , ϕ ) = ( A ai t a , ϕ a t a ) = ( A i , A ) / ( g √ T ) , G ij = ∂ i A j − ∂ j A i + ig E [ A i , A j ] , and D i ϕ = ∂ i ϕ + ig E [ A i , ϕ ] with Tr( t a t b ) = δ ab . The higher dimensional operators aredenoted by δ L EQCD . The effective coupling and the Debyescreening mass in the leading-order (LO) read g E = g √ T and m D = ( N/ gT ) , respectively.Under the “Feynman static gauge” ( ∂ A = 0 for thetemporal component and the Feynman gauge for the spa-tial component A i ) [36, 37], the tree-level propagators1read (cid:104) ϕ a ( ) ϕ b ( x ) (cid:105) = δ ab e − m D | x | π | x | , (B2) (cid:104)A ai ( ) A bj ( x ) (cid:105) = δ ab δ ij π | x | , (B3) (cid:104)A ai ( ) ϕ b ( x ) (cid:105) = 0 , (B4)where a and b are color indices. Moreover, the Polyakovloop operator Ω( x ) is written as Ω( x ) = P e − ig (cid:82) /T dτA ( x ,τ ) = e − igϕ ( x ) / √ T . (B5)The leading order (LO) contribution to the connectedcorrelation between the Polyakov loop and the EMTstems from the two gluon exchange of O ( g ) and is dia-grammatically shown in Fig. 12. Since the Polyakov loopoperator has only the scalar component ϕ ( x ) , the termswhich survive in the LO are the connected diagrams with G i , i.e., (cid:104) ( G ai ) ( x ) (cid:105) Q = (cid:104) ( G ai ) ( x )TrΩ( ) (cid:105) c (cid:104) TrΩ( ) (cid:105) , (B6)where the suffix c implies the connect correlation.By expanding ( G ai ) and TrΩ up to O ( ϕ ) for fixed i = 1 , , , we obtain (cid:104) ( G ai ) ( x ) (cid:105) Q = − N g (cid:104) ϕ a ( ) ∂ i ϕ b ( x ) (cid:105)(cid:104) ϕ b ( ) ∂ i ϕ a ( x ) (cid:105) + O ( g )= − N − N g (cid:26) ∂ i (cid:18) e − m D | x | π | x | (cid:19)(cid:27) + O ( g )= − C F π α s x i ( m D | x | + 1) | x | e − m D | x | + O ( g ) , (B7)where α s = g / π and C F = ( N − / N .Picking up the contributions of ( G a i ) ( x ) in each com-ponent of the EMT, we obtain the following perturbativeestimate for r ≡ | x | (cid:29) (2 πT ) − up to O ( g ) , (cid:104)T ( x ) (cid:105) Q = (cid:104)T rr ( x ) (cid:105) Q = −(cid:104)T θθ ( x ) (cid:105) Q = − C F π α s ( m D r + 1) r e − m D r + O ( g ) , (B8)Simplest way to show the above relation is to choose x = ( r, , , so that T rr ( r, ,
0) = T ( r, , and T θθ ( r, ,
0) = T ( r, , .Although one finds that the EMT trace (cid:104)T µµ ( x ) (cid:105) Q van-ishes at O ( g ) , one can utilize the following trace anomalyto evaluate its O ( g ) contribution: T µµ = β g G aµν G aµν , (B9)where the Yang-Mills beta function reads β = − β g − β g + · · · , with β = (11 / C A / (4 π ) , β =(34 / C A / (4 π ) , and C A = N . By using the right handside of the formula and follow the same procedure asabove, we find (cid:104)T µµ ( x ) (cid:105) Q = − C A C F (4 π ) α s ( m D r + 1) r e − m D r + O ( g ) . (B10)This is indeed O ( g ) higher than Eq. (B8). [1] G. S. Bali, Phys. Rept. , 1 (2001).[2] R. Narayanan and H. Neuberger, JHEP , 064 (2006)[hep-th/0601210].[3] M. Lüscher, JHEP , 071 (2010) [arXiv:1006.4518[hep-lat]].[4] M. Lüscher and P. Weisz, JHEP , 051 (2011)[arXiv:1101.0963 [hep-th]].[5] H. Suzuki, PTEP , 083B03 (2013) [Erratum: PTEP , 079201 (2015)] [arXiv:1304.0533 [hep-lat]].[6] H. Makino and H. Suzuki, PTEP , 063B02 (2014)[Erratum: PTEP , 079202 (2015)] [arXiv:1403.4772[hep-lat]].[7] K. Hieda, H. Makino and H. Suzuki, Nucl. Phys. B ,23-51 (2017) [arXiv:1604.06200 [hep-lat]].[8] R. V. Harlander, Y. Kluth and F. Lange, Eur. Phys. J. C , 944 (2018) [arXiv:1808.09837 [hep-lat]].[9] R. Yanagihara, T. Iritani, M. Kitazawa, M. Asakawaand T. Hatsuda, Phys. Lett. B , 210 (2019)[arXiv:1803.05656 [hep-lat]].[10] R. Yanagihara and M. Kitazawa, PTEP ,093B02 (2019) [Erratum: PTEP , 079201 (2020)][arXiv:1905.10056 [hep-ph]].[11] M. Asakawa et al. [FlowQCD], Phys. Rev. D ,011501 (2014) [Erratum: Phys. Rev. , 059902 (2015)][arXiv:1312.7492 [hep-lat]].[12] M. Kitazawa, T. Iritani, M. Asakawa, T. Hatsuda,and H. Suzuki, Phys. Rev. D , 114512 (2016)[arXiv:1610.07810 [hep-lat]].[13] M. Kitazawa, T. Iritani, M. Asakawa, and T. Hatsuda,Phys. Rev. D , 111502 (2017) [arXiv:1708.01415 [hep- lat]].[14] T. Iritani, M. Kitazawa, H. Suzuki and H. Takaura,PTEP , 023B02 (2019) [arXiv:1812.06444 [hep-lat]].[15] T. Hirakida, E. Itou and H. Kouno, PTEP , 033B01(2019) [arXiv:1805.07106 [hep-lat]].[16] M. Kitazawa, S. Mogliacci, I. Kolbé and W. Horowitz,Phys. Rev. D , 094507 (2019) [arXiv:1904.00241 [hep-lat]].[17] Y. Taniguchi, S. Ejiri, R. Iwami, K. Kanaya, M. Ki-tazawa, H. Suzuki, T. Umeda, and N. Wakabayashi,Phys. Rev. D , 014509 (2017) [arXiv:1609.01417 [hep-lat]].[18] Y. Taniguchi et al. [WHOT-QCD], Phys. Rev. D ,014510 (2020) [arXiv:2005.00251 [hep-lat]].[19] L. D. Landau and E. M. Lifshitz, “The Classical The-ory of Fields” (fourth Edition) (Butterworth-Heinemann,1980).[20] S. Borsanyi, S. Dürr, Z. Fodor et al. , JHEP , 010(2012) [arXiv:1203.4469 [hep-lat]].[21] G. Parisi, R. Petronzio, and F. Rapuano, Phys. Lett. , 418 (1983).[22] O. Kaczmarek, F. Karsch, F. Zantow and P. Petreczky,Phys. Rev. D , 074505 (2004) [arXiv:hep-lat/0406036[hep-lat]].[23] O. Kaczmarek and F. Zantow, Phys. Rev. D , 114510(2005) [arXiv:hep-lat/0503017 [hep-lat]].[24] M. Berwein, in preparation.[25] See, for example, L. Müller, O. Philipsen, C. Reisinger and M. Wagner, Phys. Rev. D , no.5, 054503 (2019)[arXiv:1907.01482 [hep-lat]].[26] T. T. Takahashi, H. Suganuma, Y. Nemoto, andH. Matsufuru, Phys. Rev. D , 114509 (2002) [hep-lat/0204011].[27] S. Kumano, Q. T. Song and O. V. Teryaev, Phys. Rev.D , 014020 (2018) [arXiv:1711.08088 [hep-ph]].[28] M. V. Polyakov and P. Schweitzer, Int. J. Mod. Phys. A , 1830025 (2018) [arXiv:1805.06596 [hep-ph]].[29] V. D. Burkert, L. Elouadrhiri and F. X. Girod, Nature , 396 (2018).[30] P. E. Shanahan and W. Detmold, Phys. Rev. Lett. ,072003 (2019) [arXiv:1810.07589 [nucl-th]].[31] S. Necco and R. Sommer, Nucl. Phys. B , 328-346(2002) [arXiv:hep-lat/0108008 [hep-lat]].[32] Z. Fodor, K. Holland, J. Kuti, S. Mondal, D. Nogradiand C. H. Wong, JHEP , 018 (2014) [arXiv:1406.0827[hep-lat]].[33] L. Altenkort, A. M. Eller, O. Kaczmarek, L. Mazur,G. D. Moore and H. T. Shu, [arXiv:2009.13553 [hep-lat]].[34] P. Fritzsch and A. Ramos, JHEP , 008 (2013)[arXiv:1301.4388 [hep-lat]].[35] T. Appelquist and R. D. Pisarski, Phys. Rev. D , 2305(1981)[36] E. D’Hoker, Nucl. Phys. B , 401-428 (1982)[37] S. Nadkarni, Phys. Rev. D , 917 (1983)[38] E. Braaten and A. Nieto, Phys. Rev. Lett.74