Distribution of Stress Tensor around Static Quark--Anti-Quark from Yang-Mills Gradient Flow

The spatial distribution of the stress tensor around the quark--anti-quark ($Q\bar{Q}$) pair in SU(3) lattice gauge theory is studied. The Yang-Mills gradient flow plays a crucial role to make the stress tensor well-defined and derivable from the numerical simulations on the lattice. The resultant stress tensor with a decomposition into local principal axes shows, for the first time, the detailed structure of the flux tube along the longitudinal and transverse directions in a gauge invariant manner. The linear confining behavior of the $Q\bar{Q}$ potential at long distances is derived directly from the integral of the local stress tensor.

The energy-momentum tensor (EMT), T µν (x), in classical and quantum field theories is a special quantity among other observables in the sense that it relates the local properties and the global behaviors of the system in a gauge invariant manner. A classic example is the Maxwell stress-tensor in electromagnetism: σ Maxwell [1]. It describes the local response under external charges, and its integration on the surface surrounding a charge gives the Coulomb force acting on the charge. In quantum Yang-Mills (YM) theory, the EMT is even more important than in the Abelian case, since it provides gauge-invariant and non-perturbative information.
The purpose of this Letter is to explore novel aspects of EMT in YM theory at zero temperature under the presence of static quark (Q) and anti-quark (Q) charges separated by a distance R. In such a setup, the YM field strength is believed to be squeezed into a quasi-one-dimensional flux tube [2] and gives rise to the linear confining potential at large R (see the reviews [3][4][5] and references Email address: yanagihara@kern.phys.sci.osaka-u.ac.jp (Ryosuke Yanagihara) therein). Although the action density, the color electric field and the plaquettes have been employed before to probe such a flux tube [6][7][8][9][10][11][12][13], the present Letter is a first attempt to provide gauge invariant EMT distribution around the QQ pair in three spatial dimensions. The fundamental theoretical tool to make this analysis possible is the YM gradient flow [14][15][16], which was recently put in practice to treat T µν (x) [17,18] and has been applied extensively to the equation of state of SU(3) YM theory at finite temperature [19][20][21][22].
Before going into the details of our lattice study, let us first discuss the general feature of T µν (x) in the Euclidean spacetime with µ, ν = 1, 2, 3, 4. The local energy density and the stress tensor read respectively as The force per unit area F i , which induces the momentum flow through a given surface element with the normal vector n i , is given by [1] Then the local principal axes and the corresponding eigenvalues of the local stress can be obtained by diagonalizing T ij : where the strengths of the force per unit area along the principal axes are given by the absolute values of the eigenvalues, λ k . The force acting on a test charge is obtained by the surface integral F i = − S T ij dS j , where S is a surface surrounding the charge with the surface vector S j oriented outward from S. In quantum YM theory, obtaining T µν (x) nonperturbatively around static QQ on the lattice requires us to go through the following steps.
The first step is to start with the YM gradient flow equation [15], with the fictitious 5-th coordinate t. The YM action S YM (t) is composed of A µ (t, x), whose initial condition at t = 0 is the ordinary gauge field A µ (x) in the four dimensional Euclidean space. The gradient flow for positive t smooths the gauge field with the radius √ 2t. Then the renormalized EMT operator is defined as [17] Here E(t, x) = (1/4)G a µν (t, x)G a µν (t, x) and U µν (t, x) = G a µρ (t, x)G a νρ (t, x) − δ µν E(t, x) with the field strength G a µν (t, x) composed of the flowed gauge field A µ (t, x). The vacuum expectation value T µν (t, x) 0 is normalized to be zero due to the subtraction of E(t, x) 0 . We use the perturbative coefficients for α U (t) and α E (t) [17] in the following analysis. Thermodynamic quantities in SU(3) YM theory have been shown to be accurately obtained with this EMT operator with smaller statistics than with the previous methods [19][20][21].
The second step is to prepare a static QQ system on the lattice. We use the rectangular Wilson loop W (R, T ) with static color charges at R ± = (0, 0, ±R/2) and in the temporal interval [−T /2, T /2]. Then the expectation value of T µν (t, x) around the QQ is obtained by [23] T µν (t, x) QQ = lim where T → ∞ is to pick up the ground state of QQ.
The measurements of T µν (t, x) for different values of t are made at the mid temporal plane The final step is to obtain the renormalized EMT distribution around QQ from the lattice data by taking double limit [19,20], In lattice simulations we measure T µν (t, x) lat QQ at finite t and a, and make an extrapolation to (t, a) = (0, 0) according to the formula [20,22], where b µν (t) and c µν are contributions from lattice discretization effects and the dimension six operators, respectively. The numerical simulations in SU(3) YM theory are performed on the four-dimensional Euclidean lattice with the Wilson gauge action and the periodic boundary condition. Shown in Table 1 are five different inverse couplings β = 6/g 2 0 and corresponding lattice spacings a determined from the w 0scale [20,25]. The lattice size N 4 size , and the number of gauge configurations N conf are also summarized in the table. Gauge configurations are generated by the same procedure as in [20] with the separation of 200 (100) sweeps on the 64 4 (48 4 ) lattice. 1 Statistical errors are estimated by the jackknife method with 20 jackknife bins at which the errors saturate. In the flow equation, the Wilson gauge action is used for S YM (t), while the clover-type representation is adopted for G a µν (t, x) in T µν (t, x). Other than the gradient flow for T µν (t, x) described above, we adopt the standard APE smearing for each spatial link along the Wilson loop [26] with the same smearing parameter as in [27] to enhance the coupling of W (R, T ) to the QQ ground state. We keep a √ N APE which is proportional to the transverse size of the spatial links [14,27] to be approximately constant by changing the iteration number N APE . Also, to reduce the statistical noise, we adopt the standard multi-hit procedure by replacing each temporal link by its mean-field value [7,28]. We consider three QQ distances (R = 0.46, 0.69, 0.92 fm), which are comparable to the typical scale of strong interaction. These values as well as the corresponding dimensionless distances R/a are summarized in Table 1. While the largest R is half the spatial lattice extent aN size for the two finest lattice spacings, effects of the periodic boundary are known to be well suppressed even with this setting [7]. A measure of the ground state saturation in the QQ system reads with the ground-state potential V (R) and the ground-state overlap C 0 (R) obtained at large T [7].
Using the data at a = 0.038 fm with N APE = 160, we found |1 − P (R, T )| < 0.5% as long as T > 0.19 fm for all R in Table 1. By keeping T 0.46 fm to extract observables as shown in the last column of Table 1, the ground state saturation of the Wilson loop is, therefore, secured in our simulations. 2 To avoid the artifact due to finite a and the oversmearing of the gradient flow [20,22], we need to choose an appropriate window of t satisfying the condition a/2 ρ L. Here ρ ≡ √ 2t is the flow radius, and L ≡ min(| x − R + |, | x − R − |, T /2) is the minimal distance between x µ = ( x, 0) and the Wilson loop.
Before taking the double limit in Eq. (9), we illustrate a qualitative feature of the distribution of T ij 2 An alternative estimate of the ground state saturation is obtained through the excitation energy of a bosonic string: ∆En = πn/R (n =, 1, 2, · · · ) [29]. By taking n = 2, which corresponds to the excitation with the same symmetry (Σ + g ) as the ground state [30], the excited state is suppressed at least by a factor exp(−(2π/R)T ) = exp(−π) 4% for T = 0.46 fm and R = 0.92 fm, even if C n=2 (R) ∼ C n=0 (R).
around QQ in YM theory at fixed a and t by considering the case a = 0.029 fm with R = 0.69 fm, T = 0.46 fm, and t/a 2 = 2.0. Shown in Fig. 1 are the two eigenvectors in Eq. (4) along with the principal axes of the local stress. The other eigenvector is perpendicular to the y-z plane. The eigenvector with negative (positive) eigenvalue is denoted by the red outward (blue inward) arrow with its length proportional to |λ k |: Neighbouring volume elements are pulling (pushing) with each other along the direction of red (blue) arrow according to Eq. (3). The spatial regions near Q andQ, which would suffer from over-smearing, are excluded in the figure. Spatial structure of the flux tube is clearly revealed through the stress tensor in Fig. 1 (a) in a gauge invariant way. This is in contrast to the same plot of the principal axes of T ij for opposite charges in classical electrodynamics given in Fig. 1 (b).
Let us now turn to the mid-plane between the QQ pair and extract the stress-tensor distribution at x = (x, y, 0) by taking the double limit Eq. (9). In Fig. 2 (a), we show an example of the T dependence of T zz (t, 0) lat QQ , which gives the largest eigenvalue on the mid-plane, for several values of t with a = 0.029 fm (≡ā). The figure indicates that significant T -dependence arises owing to oversmearing of the gradient flow for t/ā 2 ∼ 6 (i.e., ρ ∼ √ 12ā = 0.10 fm) already around T = 16ā = 0.46 fm. This is so for the same t/ā 2 in all other cases in Table 1. On the other hand, under-smearing of the gradient flow on the coarsest lattice becomes significant for t/ā 2 = 2 (ρ = a = 0.058 fm). Therefore, in the following analysis, we focus on the data in the interval 2 ≤ t/ā 2 ≤ 5 (2ā ≤ ρ ≤ √ 10ā), which satisfies a/2 ρ L with margin.
In Fig. 2 (b), we show T zz (t, 0) lat QQ as a function of the dimensionless ratio a  gles for different values of t. The continuum limit (a → 0) with fixed t is taken by using these data together with the formula Eq. (10). The results are shown by the filled black squares with error bars.
In Fig. 2 (c), the open symbols with error bars correspond to the original data for various values of a and t/ā 2 . The result of the continuum limit in the interval 2 ≤ t/ā 2 ≤ 5 is denoted by the black solid line with the shaded error band. The t → 0 limit is carried out by using the values in the continuum limit according to Eq. (10) with a = 0. As a most conservative range for the extrapolation, we take 3 ≤ t/ā 2 ≤ 4 (Range-1). Also, to estimate the systematic errors from the extrapolation, we consider two different ranges by changing the upper and lower limits: 2 ≤ t/ā 2 ≤ 4 (Range-2) and 3 ≤ t/ā 2 ≤ 5 (Range-3). The resulting values of T R zz (0) QQ after the double limit are shown by the filled black symbols at t = 0. The dashed line corresponds to the extrapolation with Range-1. The double limit for general x = (x, y, 0) on the mid-plane can be carried out essentially through the same procedure with a few extra steps. First of all, the cylindrical coordinate system c = (r, θ, z) with r = x 2 + y 2 and 0 ≤ θ < 2π is useful for the present QQ system. On the mid-plane we have T cc (t, x) lat QQ = diag( T rr (t, r) lat QQ , T θθ (t, r) lat QQ , T zz (t, r) lat QQ ). Next, we need data at the same r for different a to take the continuum limit.
We consider the values of r at which the lattice data are available on the finest lattice. To obtain the data at these r on lattices with different a, we interpolate the lattice data T cc (t, r) lat QQ and T 44 (t, r) lat QQ with the commonly used functions to parametrize the transverse profile of the flux tube: f Bessel (r) = A 0 K 0 √ Br 2 + C with the 0th-order modified Bessel function K 0 (x) [31] and f exp (r) = (A 0 + A 1 r 2 )e (−2 √ r 2 +B 2 +2B)/C [13]. Once it is done, the t → 0 limit is taken in the same way as explained above.
In Fig. 3, we show the r dependence of the resulting EMT (the stress tensor − T R cc (r) QQ and the energy density − T R 44 (r) QQ ). From the figure, one finds several noticeable features: (i) Approximate degeneracy between temporal and longitudinal components is found for a wide range of r: T R 44 (r) QQ T R zz (r) QQ < 0. This feature is compatible with the leadingorder prediction of the worldsheet theory of QCD string [11]. We also find T R rr (r) QQ T R θθ (r) QQ > 0, which does not have simple interpretation except at r = 0.
(ii) The scale symmetry broken in the YM vacuum (the trace anomaly) is partially restored inside the flux tube, which arises in the numerical results, T R µµ (r) QQ = T R 44 (r) + T R zz (r) + T R rr (r) + T R θθ (r) QQ < 0. This is in sharp contrast to the case of classical electrodynamics; T 44 (r) = T zz (r) = −T rr (r) = −T θθ (r) and T µµ (r) = 0 for all r.
(iii) Each component of EMT at r = 0 decreases as R becomes larger, while the transverse radius of the flux tube, typically about 0.2 fm, seems to increase for large R [10,13,23], although the statistics is not enough to discuss the radius quantitatively.   Finally, we consider a non-trivial relation between the force acting on the charge located at z > 0 evaluated by the QQ potential through F pot = −dV (R)/dR (13) and the force evaluated by the surface integral of the stress-tensor surrounding the charge, For F pot , we fit the numerical data of V (R) obtained from the Wilson loop at a = 0.038 fm by the Cornell parametrization, V Cornell (R) = −A/R+ σR + B. Note that V (R) at this lattice spacing is shown to be already close to the continuum limit [3].
For F stress , we take the mid-plane for the surface integral: F stress = 2π ∞ 0 T zz (r) QQ rdr. Here T zz (r) QQ is obtained by fitting Fig. 3 with either f Bessel (r) or f exp (r). In Fig. 4, −F pot and −F stress thus obtained are shown by the solid line and the horizontal bars, respectively. −F pot increases as A/R 2 at short distance and approaches the string tension σ at large distance. For −F stress , we take into account not only the statistical error but also the systematic errors from the double limit and the fitting in terms of f Bessel,exp (r). The agreement between the two quantities within the errors is a first numerical evidence that the "action-at-a-distance" QQ force can be described by the local properties of the stress tensor in YM theory.
In this Letter, we have performed a first study on the spatial distribution of EMT around the QQ system in SU(3) lattice gauge theory. The EMT operator defined through the YM gradient flow plays a crucial role here. The transverse structure of the stress-tensor distribution in the mid-plane is analyzed in detail by taking the continuum limit and the zero flow-time limit successively. The linear confining behavior of the QQ potential at long distances can be shown to be reproduced by the surface integral of the stress tensor. Further details of the stress-tensor distribution, not only in the transverse direction but also in the longitudinal direction, and R dependence of the transverse radius [10,13,23], will be reported in a forthcoming publication [32]. There are also interesting future problems to be studied on the basis of the formalism presented in this Letter: Those include the applications to the QQQ system [27] and to the QQ system at finite temperature [33] as well as the generalization to full QCD [34] with the QCD flow equation [35,36].
The authors thank K. Kondo, F. Negro, A. Shibata, and H. Suzuki for discussions. Numerical simulation was carried out on IBM System Blue Gene Solution at KEK under its Large-Scale Simulation Program (No. 16/17-07). This work was supported by JSPS Grant-in-Aid for Scientific Researches, 17K05442, 25287066 and 18H05236. T.H. is grateful to the Aspen Center for Physics, supported in part by NSF Grants PHY1607611.