Distribution of Energy-Momentum Tensor around a Static Quark in the Deconfined Phase of SU(3) Yang-Mills Theory

Energy momentum tensor (EMT) characterizes the response of the vacuum as well as the thermal medium under the color electromagnetic fields. We define the EMT by means of the gradient flow formalism 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 EMT distributions in the radial and transverse directions except for the sign, the temporal component is substantially different from the spatial ones near the critical temperature $T_c$. This is in contrast to the prediction of the leading-order thermal perturbation theory. The lattice data of the EMT distribution also indicate the thermal screening at long distance and the perturbative behavior at short distance.


I. INTRODUCTION
To study complex quantum systems such as the Yang-Mills (YM) theory, it is customary to introduce test probe(s) and analyze the response. The Wilson loop is one of such probes whose measurement in YM theory provides information on the static quark-anti-quark system that is closely related to the confinement property in YM vacuum [1]. Thanks to the recent development of the gradient-flow method [2][3][4] and its application to the energy-momentum tensor (EMT) T µν (x) [5][6][7][8], it became possible to study the gauge-invariant structure of the flux tube between the quark and anti-quark in the confining phase through the spatial distribution of EMT under the Wilson loop [9,10].
The purpose of the present paper is to extend the above idea and to explore the EMT distribution around a static quark in YM theory. As a first step, we consider the deconfined phase above the critical temperature T c of the SU(3) YM theory in the range of temperature 1.2 ≤ T /T c ≤ 2.6 and measure the EMT distribution around the Polyakov loop. The EMT with the gradient flow has been used to study thermodynamics of YM theory [11][12][13][14][15][16] and of QCD [17,18]. However, the observables in these studies are limited to global quantities such as the pressure, energy density, entropy density, and the specific heat. On the other hand, we focus on the local observable in this study and examine the following questions: (i) How are the energy density and the stress tensor distributed around the static quark?, (ii) How are the distributions modified as a function of temperature?, and (iii) How can one extract parameters such as the running coupling and the Debye screening mass from the distributions?
The organization of the present paper is as follows. In Sec. II, we briefly review the definition of EMT and its property in the spherical coordinate system. In Sec. III, we introduce the EMT operator on the lattice and its correlation with the Polyakov loop operator. In Sec. IV, we discuss the numerical procedure and lattice setup to analyze the EMT operator around a static quark on the lattice. Numerical results and their physical implications are given in Sec. V. Sec. VI is devoted to the summary and conclusion. In Appendix A, we discuss the procedure to make the tree-level improvement of the correlation between the EMT and the Polyakov loop on the lattice. In Appendix B, the leading-order perturbative analysis of the correlation is presented using the high temperature effective field theory.

II. EMT AROUND A STATIC CHARGE
The stress tensor σ ij (i, j = 1, 2, 3) is related to the spatial component of EMT, T ij , as [19] σ ij = −T ij . (1) The force per unit area F i acting on a surface with the normal vector n i is given by the stress tensor as The local principal axes n (k) j and the corresponding eigenvalues λ k of the local stress tensor are obtained by solving the eigenvalue problem: The strength of the force per unit area along n (k) i is given by the absolute values of the eigenvalue λ k . Neighboring volume elements separated by a surface with the normal vector n (k) i pull (push) each other for λ k < 0 (λ k > 0) across the wall. Note that the three principal axes n Stress acting on infinitely small volume element under the existence of a single static charge (source). In the case of the classical electromagnetism, a small volume element at a distance of r from the charge is pulled along the radial direction from the neighboring volume elements while it is pushed in the transverse direction.
are orthogonal with one another because σ ij is a symmetric tensor.
For a system with a single static source, it is convenient to use the spherical coordinate system (r, θ, ϕ) with the radial coordinate r = |x|, and the polar and azimuthal angles θ and ϕ. The spherical symmetry allows us to diagonalize the static EMT in Euclidean spacetime T µν (x) in this coordinate system as where γ, γ = 4, r, θ. Due to the spherical symmetry, the azimuthal component degenerates with the polar component, T ϕϕ (r) = T θθ (r), so that only independent components are given in Eq. (4).
In the Abelian case, EMT is given by the Maxwell stress-energy tensor [19], T Maxwell µν = F µρ F νρ − 1 4 δ µν F ρσ F ρσ with the field strength F µν . When a static charge is placed at the origin, the EMT is denoted by with E i (x) being the electric field. The spatial structure of Eq. (5) is illustrated in Fig. 1 where the neighboring volume 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 element through its surface should be balanced. This property is guaranteed by the momentum conservation ∂ i T ij = 0 together with the Gauss theorem. We consider the pure SU(3) YM gauge theory in the four-dimensional Euclidean space defined by the action, Here g 0 is a bare gauge coupling and G a µν (x) is the field strength composed of the fundamental gauge field A a µ (x). The YM gradient flow evolves the gauge field along the fictitious fifth dimension t introduced in addition to the ordinary four Euclidean dimensions x through the flow equation [2][3][4], The flowed YM action S YM (t) in the (4 + 1)-dimensional coordinate is constructed by substituting the flowed gauge 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 > 0 is that any composite operators composed of flowed gauge fields are UV finite even at equal spacetime point [4,7]. This is a consequence of the smoothing of the gauge fields in the four dimensional Euclidean space within the range ∼ √ 2t. In addition, in the small t limit, composite local operators are represented by the local operators of the ordinary gauge theory at t = 0. These properties lead us to the renormalized EMT operator defined with the small t expansion [5]: where E(t, x) 0 is the vacuum expectation value of E(t, x). The dimension-four gauge-invariant operators on the right hand side of Eq. (9) are given by [5] where G a µν (t, x) is the field strength composed of the flowed gauge field. Because of the vacuum subtraction in Eq. (9), T R µν (x) 0 vanishes. The coefficients c 1 (t) and c 2 (t) have been calculated perturbatively in Refs. [5,8,14] for small t. We use two-loop perturbative coefficients [8,14] for the construction of EMT throughout this study.

B. EMT around a static heavy quark
To describe a static quark Q on the lattice, we introduce the Polyakov loop at the origin, Ω(0). Then the expectation value of Eq. (9) around Q is given by We note that Eq. (12) is well-defined only when the Z 3 symmetry in SU(3) YM theory is spontaneously broken: In the Z 3 unbroken phase, both numerator and denominator of the first term on the right hand side vanish exactly. This is the reason why we focus on the system in the Z 3 broken phase above T c in this paper. In practice, we choose the state with the Polyakov loop being real among the three equivalent Z 3 states in the deconfined phase.
The renormalized EMT distribution around Q is obtained after taking the double extrapolation, In our actual analysis, we extract the renormalized EMT distribution by fitting the lattice data with the following functional form [12,13]: 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, the smearing radius ρ ≡ √ 2t needs to be larger than the lattice spacing to suppress the discretization error. At the same time, ρ should be smaller than half the temporal size 1/2T with temperature T as well as the distance from the source (r) to avoid the overlap of operators. Therefore we require a/2 ρ min r, 1 2T .
The lattice data to be fitted by Eq. (14) should be within this window. As will be discussed in Sec. IV C, we impose more stringent conditions for the range of t in our numerical analysis.

A. Gauge configurations
Numerical simulations in SU(3) YM theory were performed on the four dimensional Euclidean lattice with the Wilson gauge action and the periodic boundary conditions at four different temperatures 1.20T c , 1.44T c , 2.00T c , and 2.60T c . The simulation parameters for each T are summarized in Table I. The inverse coupling β = 6/g 2 0 is related to the lattice spacing a determined by the reference scale w 0 [12,20]. The spatial and temporal lattice sizes, N s and N τ together with the number of configurations N conf are also summarized in Table I. All lattices have the same aspect ratio The gauge configurations are generated by the pseudoheat-bath method followed by five over-relaxations. Each measurement is separated by 200 sweeps. Statistical errors are estimated by the jackknife method with 20 jackknife bins. We employ the Wilson gauge action for S YM (t) in the flow equation Eq. (7) and the clover type representation for the field strength G µν (t, x). The numerical solution of the gradient flow equation is obtained by the third order Runge-Kutta method.
In order to suppress the statistical noise, we apply the multi-hit procedure in the measurement of the Polyakov loop by replacing every temporal link by its thermal average [21]. The choice of the temporal argument x 4 of EMT in Eq. (12) is arbitrary. Therefore, we average EMT over the temporal direction to reduce the statistical error.

B. Discretization effect
The EMT in the spherical coordinate system on the lattice reads The behavior of the EMT distribution close to the source Q is affected by the violation of rotational symmetry owing to lattice discretization. As an example, we show in Fig. 2 the distribution of −r 4 T rr (t/a 2 6.910 = 1.3, r) Q as a function of rT at T /T c = 1.44, where a 6.910 is the lattice spacing of the finest lattice at this temperature. The figure shows the oscillating behavior of the numerical results becomes more prominent on coarser lattices. In this study, we use the lattice data only for N τ ≥ 12 for the continuum extrapolation to suppress the discretization errors. In Appendix A, we consider an alternative analysis that performs the tree-level improvement of the numerical results and uses them for the continuum extrapolation with the N τ = 10 data. As discussed there, we confirm that the results in both cases are the same within the errors.

C. Double extrapolation
The double extrapolation Eq. (13) consists of two steps: (I) the continuum (a → 0) extrapolation, and (II) t → 0 extrapolation. In this subsection, we demonstrate these procedures by using the lattice data at T /T c = 1.44 as an example.
In Fig. 3, we show − T rr (t, rT ) Q /T 4 at rT = 0.40, 0.48, 0.60 as a function of 1/N 2 τ = (aT ) 2 for four values of tT 2 . To obtain the EMT values at a given r and t, we first perform the linear interpolation of the lattice data along the r direction and then interpolate along the t direction by the cubic spline method for each a.
In Fig. 3, fitting results of the data at N τ = 12-18 according to Eq. (14) at fixed t are shown by the solid lines, while the results of the continuum limit are shown by the filled squares on the vertical dotted line at 1/N 2 τ = 0. We then take t → 0 extrapolation by fitting the continuumextrapolated results for different t with Eq. (14) at a = 0. This fit has to be carried out within the range of t satisfying Eq. (15). We employ tT 2 = 0.00401 corresponding to t/a 2 = 1. quantities show a linear behavior below this value [12].
We consider the following three ranges within 0.00401 ≤ tT 2 ≤ 0.015 to estimate the systematic uncertainty from the fitting ranges [12]: Range-1 is the most conservative window, while Range-2 (Range-3) is the extension of Range-1 towards the smaller estimate of the systematic error. In the following, all the results after the double extrapolation contain both the statistical and systematic errors.
In Fig. 4, the open symbols with statistical errors represent T γγ (t, rT ) Q /T 4 at rT = 0.40 (upper) and 0.60 (lower) for each a. The results of the continuum limit are denoted by the black solid lines with the gray statistical error band for 0.00401 ≤ tT 2 ≤ 0.015. Range-1 is highlighted by the yellow band. The figure also shows the fitted results for Ranges-1, 2, and 3 by the dotted lines. The final results of the t → 0 limit for each range are shown by the open black symbols around tT 2 = 0: They agree with each other within the statistical errors, which suggests that the systematic uncertainty from the choice of the fitting range is not significant.

V. RESULTS OF EMT DISTRIBUTIONS
Before entering into the detailed discussions on the spatial distribution of EMT, we first show the result of the stress distribution at T /T c = 2.60 on a two-dimensional plane including the static source in Fig. 5. The same result is later shown in Fig. 8 in a different form. In Fig. 5, the red and blue arrows represent the principal directions of the stress tensor along the radial and transverse directions, respectively. The length of each arrow represents square root of the eigenvalue corresponding to each principal axis. This figure is to be compared with Fig. 1 in Ref. [9].

A. Channel dependence
In Fig. 6, we show the dimensionless EMT, − T R 44 (r) Q /T 4 , − T R rr (r) Q /T 4 , and T R θθ (r) Q /T 4 , as functions of the dimensionless length rT . The error bands include both the statistical and systematic errors, where the latter is estimated from the three fitting ranges for the t → 0 extrapolation. Since the thermal expectation value T µν (t, x) is subtracted as in Eq. (12), we have T R µν (r) Q → 0 in the r → ∞ limit. We find that − T R 44 (r) Q , − T R rr (r) Q , and T R θθ (r) Q are all positive for rT 1 and decrease rapidly with increasing r. These signs are the same as those of the Maxwell stress tensor in Eq. (5). Individual signs physically mean that a volume element has a positive localized energy density and receives a pulling (pushing) force along the longitudinal (transverse) direction; see Figs. 1 and 5. Figure 6 indicates that the absolute values of the spatial components | T R rr (r) Q | and | T R θθ (r) Q | are degenerated within the error for all temperatures. On the other hand, | T R 44 (r) Q | is larger than the spatial components especially at lower temperature. This is in contrast to the degenerate magnitude of all components in the Maxwell stress Eq. (5) and is also different from the leading-order thermal perturbation theory (Appendix B). Fig. 7 is the temperature dependence of the spatial distribution of EMT with respect to the physical distance r [fm]; (a) − T R 44 (r) Q , (b) − T R rr (r) Q , and (c) T R θθ (r) Q . Also, shown in Fig. 7(d) is the distribution of the trace of EMT given by Figure 7 tells us that the EMT distributions have small T dependence at short distances, r 0.2 fm. On the other hand, for large distances, sizable T dependence can be seen despite the growth of the errors at high T .

Shown in
To make these features more explicit, we plot the same results with a dimensionless normalization r 4 T R γγ (r) Q as a function of r in Fig. 8. The figure shows that the T dependence is suppressed for rT 0.3 and all results approach a single line, while the result tends to be more suppressed for rT 0.3 compared with this universal behavior as temperature is raised. This result is reasonable as the T dependence of r 4 T R γγ (r) Q would be suppressed for r (2πT ) −1 .
At distance (2πT ) −1 r Λ −1 with Λ being the lambda parameter, the behavior of T R γγ (r) Q should be described by the perturbation theory in electrostatic QCD (EQCD). In the leading order of EQCD in this regime, we have the following ratio (Appendix B), which is independent of r and T and is given only by a function of α s . Shown in Fig. 9 is the r dependence of ∆ Q (r)/ T R 44,rr,θθ (r) Q as a function of r at T /T c = 2.60. (d) of the Polyakov loop correlations at r = 0.1 fm [22,23]. Higher order α s corrections and the thermal corrections for EMT around a static charge to be compared with our 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 possible to extract the thermal screening of the form exp(−2m D r) with m D being the Debye screening mass. Nevertheless, Figs. 8(a-d)   faster than 1/r 4 at long distances, and the tendency is stronger at high temperatures. To draw a definite conclusion, however, higher statistical data are necessary.

VI. SUMMARY AND CONCLUDING REMARKS
In the present paper, we have studied, for the first time, the EMT distribution around a static quark at finite temperature above T c of the SU(3) YM theory on the lattice. The YM gradient flow plays crucial roles to define the EMT on the lattice and to explore its spatial structure.
The main results of this paper can be summarized as follows.
As shown in Fig. 6, we found no significant difference between the absolute magnitude of EMT along the radial direction and that of the transverse direction for all temperatures above T c . This seems to be in accordance with the leading-order thermal perturbation theory in QCD, which predicts the same magnitude for all principal components of EMT. However, we found a substantial difference between the EMT distribution in the temporal direction and that of the spatial directions, especially near T c . This indicates that there is indeed a genuine non-Abelian effect present at finite temperature, so that precise comparison with the higher-order thermal QCD calculation would be called for.
As shown in Figs. 7 and 8, all the EMT distributions have small T dependence at short distances, r 0.2 fm. Also the EMT distributions decrease faster than 1/r 4 at long distances, and the tendency is stronger at high temperatures. However, owing to the large statistical errors, we could not extract the values of the thermal Debye screening. By using the fact that the EMT distributions are T independent at short distances, we attempted to extract the strong coupling constant from the ratios between the different components of EMT. The result, α s (0.1 fm) 0.28 − 0.32, is consistent with that obtained from the similar analysis for the QQ free energy at finite T .
We have some important issues to be studied further: Going beyond the leading-order thermal QCD calculation for the EMT [24] is necessary to understand the lattice results presented in this paper. At the same time, increasing the statistics of lattice data is necessary to extract, e.g. the screening mass from the long range part of the EMT distribution.
There are also several interesting future problems. First of all, the extension to full QCD is an important next step. Since the Z 3 symmetry is explicitly broken by dynamical fermions, the present method can be applied directly to a single static quark Q, a static diquark QQ and QQ both at low and high temperatures. In particular, the single quark system in QCD at zero temperature corresponds to a heavy-light meson [25]. Secondly, the EMT distributions of the QQQ system will provide new insight into the flux tube formation in baryons [26] as well as the "gravitational" baryon structure [27][28][29][30] at zero and non-zero temperatures.

ACKNOWLEDGMENT
The authors thank T. Iritani for fruitful discussions in the early stage of this study. They are also grateful to M. Berwein for discussions regarding the perturbative analysis of the EMT distribution. M. K. thanks F. Karsch for useful discussions. The numerical simulation was carried out on OCTOPUS at the Cybermedia Center, Osaka University and Reedbush-U at Information Technology Center, The University of Tokyo. This work was supported by JSPS Grant-in-Aid for Scientific Researches, 17K05442, 18H03712, 18H05236, 18K03646, 19H05598, 20H01903. As shown in Fig. 2, there exists sizable discretization effect for the EMT distribution on coarse lattices especially for N τ = 10. In this Appendix, we attempt to reduce such discretization effects by using the tree-level lattice propagator. Similar idea has been applied to the analysis of the Polyakov loop correlations in Refs. [22,31].
In calculating the EMT distribution around a static quark Eq. (16), we need the expectation values of Eqs. (10) and (11) at nonzero flow time t. These quantities are constructed from G a µν (t, x)G a ρσ (t, x) Q , where the temporal coordinate is suppressed for notational simplicity. In the continuum theory at the tree level, this operator is calculated to be where g is the gauge coupling and N = 3 is the number of colors, with By selecting appropriate gauge fixing conditions for the gauge action and the gradient flow equation, one obtains Next, in lattice gauge theory the propagator corresponding to Eq. (A4) with the Wilson gauge action for the gauge action and the flow equation reads [32,33] with x n = an = a(n x , n y , n z ) andp i = (2/a) sin(p i /2a). When the clover-leaf operator for the discretized representation of G a µν (t, x n ) is employed, the discretized representation of G i4 (t, x n ) is given by [34] G lat Using Eq. (A3) and (A6), the tree-level improvements of Eqs. (10) and (11) denoted by the superscript 'imp' may be written as where the correction factor c(t, x n ) is defined by In Eq. (A9), the average over i is taken because generally the ratio G i4 (t, x n )/G lat i4 (t, x n ) at a lattice site x n depends on i. However, in our particular choice of discretization, i.e. the Wilson gauge actions and the cloverleaf operator, it is easily shown that G i4 (t, x n )/G lat i4 (t, x n ) does not depend on i. In this special case the average over i in Eq. (A9) is redundant. When this property is violated, the improvement of U γγ (t, x n ) Q may be replaced by the one that depends on γ and γ in place of Eq. (A8).
There is one more subtle issue about Eq. (A8). At the tree level, one easily finds that the matrix elements of U 44 , U rr and −U θθ are the same, while the actual lattice data do not necessarily satisfy such relation as shown in the main text. In our tree-level improvement, therefore, we decompose our lattice data into tree-like part and the rest, U γγ (t, Then we apply our tree-level improvement only to the first term: Below, we consider three choices of U tl γγ (t, x n ) to estimate the systematic uncertainty of this procedure: U tl γγ = U 44 Q , U rr Q , and − U θθ Q . Corresponding results for T R rr (t, r) as example are shown in Fig. 10(b), (c), and (d), respectively, together with the case without the correction, Fig. 10(a). Colored open symbols represent the data at each N τ and the gray (yellow) shade is the continuum result with (without) the data at N τ = 10. The figures show that the tree-level improvement suppresses the discretization effect at short distances in three cases, especially for N τ = 10.
Let us now compare the continuum extrapolation using the data including N τ = 10 with the tree-level improvement and that using the data without N τ = 10 and without the tree-level improvement. Shown in Fig. 11  is such a comparison for − T rr (t, rT = 0.40) Q /T 4 at T /T c = 1.44 as functions of 1/N 2 τ . Colored open triangles represent the data without the tree-level improvement (red) and with the tree-level improvement for three different prescriptions (blue, orange, and green). Filled symbols at 1/N 2 τ = 0 are the continuum extrapolation: The red squares are continuum results without N τ = 10 data discussed in the main text, while the diamonds are the continuum results with N τ = 10 data after the treelevel improvement. Taking into the uncertainly associated with the different prescriptions for the tree-level improvement, the default results without N τ = 10 in the main text are found to be consistent with the improved results including N τ = 10. read ϕ a (0)ϕ b (x) = δ ab e −m D |x| 4π|x| , (B2) where a and b are color indices. Moreover, the Polyakov loop operator Ω(x) is written as The leading order (LO) contribution to the connected correlation between the Polyakov loop and the EMT stems from the two gluon exchange of O(g 2 ) and is diagrammatically shown in Fig. 12. Since the Polyakov loop operator has only the scalar component ϕ(x), the terms which survive in the LO are the connected diagrams with G 2 4i , i.e., where the suffix c implies the connect correlation.