Glueball spectroscopy in lattice QCD using gradient flow

Removing ultraviolet noise from the gauge fields is necessary for glueball spectroscopy in lattice QCD. It is known that the Yang-Mills gradient flow method is an alternative approach instead of link smearing or link fuzzing in various aspects. In this work we study the application of the gradient flow technique to the construction of the extended glueball operators. We examine a simple application of the gradient flow method, which has some problems in glueball mass calculations at large flow time because of its nature of diffusion in space-time. To avoid this problem, the spatial links are evolved by the ``spatial gradient flow'', that is defined to restrict the diffusion to spatial directions only. We test the spatial gradient flow in calculations of glueball two-point functions and Wilson loops as a new smearing method, and then discuss its efficiency in comparison with the original gradient flow method and the conventional method. Furthermore, to demonstrate the feasibility of our proposed method, we determine the masses of the three lowest-lying glueball states, corresponding to the $0^{++}$, $2^{++}$ and $0^{-+}$ glueballs, in the continuum limit in the pure Yang-Mills theory.


I. INTRODUCTION
The existence of composite states consisting solely of gluons, called glueballs, is one of the important predictions of QCD. Since none of them have been identified in experiments as a glueball state, the lattice QCD results play an essential role in studying properties of the glueball states including their masses. However, ultraviolet noise from the gauge fields makes it difficult to calculate the glueball spectrum in lattice QCD. Therefore, noise reduction techniques such as the single-link smearing or the double-link fuzzing procedure plays an increasingly important role in the construction of the extended glueball operators.
Various smearing techniques such as APE smearing [1], HYP-smearing [2], and stout smearing [3] are developed for many purposes, while the fuzzing approach is proposed for a specific purpose that requires a significant improvement of having a better overlap with the glueball ground states [4]. In several previous works on the lattice glueball mass calculations [5][6][7][8], some sophisticated combinations of both the single-link smearing and the double-link fuzzing schemes are conventionally adopted (denoted as the conventional approach).
Recently, it was found that the Yang-Mills gradient flow method [9] is an alternative approach instead of the smearing in various aspects (e.g., the computation of topological charge [10]). Indeed, the gradient flow equation can be regarded as a continuous version of the recursive update procedure in the stout-link smearing with the small smearing parameter [9,10]. Therefore, in this study, we investigate the application of the gradient flow technique to the glueball calculation and also demonstrate its effectiveness in comparison to the conventional approach. This paper is organized as follows. In Sec. II, after a brief introduction of the original Yang-Mills gradient flow, we describe our proposal of "spatial gradient flow" as a new smearing method. In Sec. III, we give a short outline of how to construct glueball two-point functions based on spacelike Wilson loops. In Sec. IV, we first briefly summarize the numerical ensembles used in this study. Then we present results of the static quark potential computed with the spatial gradient flow on every ensemble. Section V gives the features of the spatial gradient flow in glueball spectroscopy. The results of glueball masses obtained by the spatial gradient flow are summarized in Sec. VI. Finally, we close with summary in Sec. VII.

II. CALCULATION METHOD I: GRADIENT FLOW METHOD
In this section, we first provide a brief review of the original Yang-Mills gradient flow, which makes the Wilson flow diffused in the four-dimensional space-time, in Sec. II A. As reported in Ref. [11], a simple application of the gradient flow technique to the glueball calculation has some problem in measuring the glueball mass from the two-point function.
We thus propose the "spatial gradient flow" as a new smearing method, which is described in Sec. II B.

A. Original gradient flow
The Yang-Mills gradient flow on the lattice is a kind of diffusion equation, where the link variables U µ (x) evolve smoothly as a function of fictitious time τ (denoted as flow time) [9].
The associated flow V µ (x, τ ) of the link variables (hereafter called the Wilson flow) driven by the gradient of the action with respect the link variables [9]. The simplest choice of the action of the link variables U µ (x) is the standard Wilson plaquette action, where g 0 is the bare coupling. The flowed link variables V µ (x, τ ) are defined by the following equation with the initial conditions V µ (x, 0) = U µ (x), where S W [V ] denotes the standard Wilson plaquette action in terms of the flowed link variables (see Appendix A for the definition of the link derivative operator ∂ x,µ and the explicit expression of ∂ x,µ S W ).
According to Eq. (2), the link variables are diffused in the four-dimensional space-time, so that the Wilson flow is approximately spread out in a Gaussian distribution with the diffusion radius (or length) of R d = √ 8τ [9]. Although such smearing procedure works well with the longer flow time, too much smearing will destroy or hide the true temporal correlation of the glueball two-point function due to the overlap of two glueball operators given by the Wilson flow as discussed in Ref. [11]. Therefore, the longer flow is not applicable for the glueball spectroscopy to avoid over smearing, that was observed in Ref. [12].
B. Spatial gradient flow as a new smearing method As described in Sec II A, the previous attempt to apply the gradient flow to the glueball spectroscopy is not fully satisfactory [12]. We propose the "spatial gradient flow" as a new smearing method in order to overcome the limited usage of the Wilson flow due to over smearing. The spatial gradient flow is defined to restrict the diffusion to spatial directions only, so that the spatial links U i (x) are evolved into the spatial Wilson flow V i (x, τ ) as the initial conditions of V i (x, 0) = U i (x) in the following gradient flow equation: Here, S sW denotes the spatial part of the standard Wilson plaquette action, where the plaquette values are composed only of the spatial links. The indices i and j run only over spatial directions. Since the spatial Wilson flow is diffused only in threedimensional space, its diffusion radius is given by R d = √ 6τ . We will later show that this new smearing works well even for the glueball spectroscopy without over smearing.

III. CALCULATION METHOD II: GLUEBALL TWO-POINT FUNCTION
We are interested in three lowest-lying glueball states, which carry specific quantum numbers, J P C = 0 ++ (scalar), 0 −+ (pseudoscalar), or 2 ++ (tensor), in this study. In this section, we briefly describe how to construct two-point correlation functions of the desired glueball state having spin J, parity P , and charge-conjugation parity C. and A 2 ), one two-dimensional representation (denoted as E), and two three-dimensional representations (denoted as T 1 and T 2 ) [13].
The inclusion of inversion (X is mapped to −X) results in the symmetry group known as O h , which has 48(=24×2) symmetry operations. The irreducible representations of O h are obtained from those of O by appending the index g (gerade) or u (ungerade), which indicates even or odd parity [14]. For convenience, we use the indices + and −, instead of g and u.
Therefore, ten different irreducible representations of O h are denoted as R P , hereafter. The glueball states are also eigenstates of charge conjugation. Therefore, the quantum number of the lattice glueball state is specified by R P C . The quantum number R P C is expected to have the following correspondence: 0 in the continuum limit [13,14].
in Table I. In the case of SU (N ) with N ≥ 3, the real part of Wilson loops has a charge conjugation parity C = +1, while the imaginary part has C-parity C = −1 [14]. In this study, we restrict ourselves to consider the real part of Wilson loops, since the three lowestlying glueball states carry C = +1.
We take the following procedure to get the irreducible contents of the representation R P with fixed C-parity from the operators O k according to Ref. [14]. First, 48 symmetry operationsŜ i , which are defined in Table II where the characters χ Γ (S i ) of the irreps Γ = A + 1 , A − 1 , E + and T + 2 are listed in Table III.
We next construct two-point correlation functions of glueball states for given irreps Γ as where the tilde over O Γ k implies the vacuum-subtracted operator defined asÕ Γ . We may also consider an N ×N correlation matrix using a set of different shape operators for given irreps Γ [15] as which allow us to perform the variational method [16,17].

A. Gauge ensembles
We perform the pure Yang-Mills lattice simulations using the standard Wilson plaquette action with a fixed physical volume (La ≈ 1.6 fm) at four different gauge couplings (β = 6/g 2 0 = 6.2, 6.4, 6.71 and 6.93). The gauge configurations in each simulation are separated by n update sweeps after n therm sweeps for thermalization as summarized in Table IV. Each sweep consists of one heat bath [18] combined with four over-relaxation [19] steps. The number of configurations analyzed is O(500-4000). All lattice spacings are set by the Sommer scale (r 0 = 0.5 fm) [20,21].
For both original and spatial gradient flows, the forth-order Runge-Kutta scheme is used with an integration step size of = 0.025. The flow time τ is given by n flow × where n flow denotes the number of flow iterations. Our simulation code for the gradient flow had been already checked in Table II of Ref. [22], where the values of the gradient flow reference scale are directly compared with the results given in the original work of Lüscher [9].

B. Scale setting and static quark potential
The static potential V (r) between heavy quark and antiquark pairs, which are separated by relative distance r, is calculated from the Wilson-loop expectation value W (r, t) with spatial extent r and temporal extent t as  .   [20,21]. where the ellipsis denotes some contribution from the excited states.
To determine the static quark potential V (r), let us consider the following quantity: Since the t-dependence of V (r, t) is expected to disappear for sufficiently large t, the static potential V (r) can be determined from a plateau seen in V (r, t) as t increases for fixed r. At glance, spatial gradient flow provides the better behavior, where the plateau starts at earlier t and extends to larger t, comparing with APE smearing. It indicates that the systematic uncertainties stemming from the excited-state contamination are better under control to determine the static potential using the spatial gradient flow method.
Hereafter, we adopt the spatial gradient flow method to evaluate the value of V (r) from the Wilson-loop expectation value, and then aim to determine the Sommer scale from the resulting static potential at each β. In this study, the Wilson loops are restricted to on-axis loops only. To extract the value of V (r) from W (r, t) at fixed r, we use the doubleexponential fit, where the a correlation among W (r, t) at different value of t is taken into account by using a covariance matrix, for our final analysis. To apply a tree-level improvement on the static quark potential, we consider is the (scalar) lattice propagator in three dimensions [21].
In Fig. 3, all data points of V I (R), which are computed at four different lattice spacings, are plotted as a function of R. The vertical and horizontal axes are normalized by the Sommer scale r 0 given in Ref. [21]. For clarity of the figure, the self-energy contribution is subtracted by the value at R = r 0 . Our results of V I (R) obtained by spatial gradient flow exhibit good scaling behavior with the literature values of r 0 /a [21]. We finally determine r 0 from our results of the static quark potentials computed at four gauge couplings. For this purpose, we first adopt the Cornell potential parametrization by fitting the data of V I (R) as with the Coulombic coefficient A, the string tension σ, and a constant V 0 . Finally, the parameters are A and σ can be used to determine the Sommer scale r 0 as for each gauge coupling β. In Table V The origin for the underestimation of r 0 /a is twofold. According to Eq. (12), one is the overestimation of A, while the other is the overestimation of σ. In the string regime (R > 1 fm) [23], the effective string theory predicts that the coefficient A is given by the universal Lüscher constant A = π/12 [24], which is smaller than our estimates of A. Thus, we alternatively choose fits with A fixed at π/12, though our data falls outside of this range. Nevertheless, the obtained results for the string tension σ becomes slightly larger than the values tabulated in Table V tension since the excited-state contaminations are not fully eliminated in our analysis of V (r) especially for large r. We simply use the double-exponential fit to determine V (r) from W (r, t) instead of the variational method that was adopted in Ref. [21].   that if the two-point function C(t) forms a Gaussian shape, with a Gaussian width corresponding to the diffusion radius R d , its effective mass gives rise to a peculiar t-dependence as whose value linearly increases from zero with a coefficient of 1/R 2 d as a function of the time slice t = t + 1 2 . This feature can be observed in the right panel of Fig. 4, where each effective mass 1 approximately starts from zero and linearly raise up to around t ≈ R d with increasing of the time slice t. Furthermore, as expected in Eq. (13) 1 To take into account "the wrap-around effect" due to the periodic boundary condition, the effective masses M eff (t ) are given by a solution of C(t) in the right panels of Fig. 4 and Fig. 5.   In the right panel of Fig. 4 and As will be discussed in Appendix A, the spatial gradient flow is slightly more efficient than the gradient flow, which is diffused in the four-dimensional space-time, in terms of reduction of relative uncertainties. As reported in an earlier work [25], although the cooling method that can smoothen the whole four-dimensional space-time was also used for calculating the string tension and glueball masses, the similar conclusion is made that the results were not better than the conventional approach that can smoothen only the three-dimensional space.
B. Equivalence to the stout smearing We will later show numerical equivalence between the spatial gradient flow and the stout smearing in the glueball calculations. As emphasized in Ref. [3], the stout smearing is a relatively new type of smearing technique, which can keep the differentiability with respect to the link variables during the smearing procedure. This property is maintained by the use of the exponential function in the stout-link smearing algorithm to remain within the SU (3) group. For the gradient flow, the numerical integrations of Eqs. (2) and (3)  We will later rederive the above matching relation in more rigorous manner in Appendix A.
In Fig. 7, we show the effective masses of the A ++ 1 glueball state obtained from the spatial gradient flow and the stout smearing at the same flow time τ that is determined by the matching relation, τ = ρn st , between the two methods. In this work, for the stout smearing, the spatially isotropic three-dimensional parameter set is chosen to be ρ ij = ρ = 0.1 and ρ 4µ = ρ µ4 = 0. The numerical equivalence between the two methods is clearly observed in both the lower and higher-diffusion cases as shown in Fig. 6. It is worth remarking that the values of n st adopted in Fig. 7, is much larger than a typical value of less than ten in the usual usage. Although the usage of the stout smearing with a small value of n st is not effective for the glueball calculations, the almost identical result to the one made by the spatial gradient flow with the diffusion radius ( √ 6τ ) can be obtained by the case if the same amount of the diffusion radius ( √ 6ρn st ) is adopted in the stout smearing.

VI. RESULTS
In this section, we present the results of glueball masses in the four channels (A ++ 1 , A −+ 1 , E ++ and T ++ 2 ) using the spatial gradient flow. To perform the continuum extrapolation, we calculate the glueball masses at four different gauge couplings with a fixed physical volume (La ≈ 1.6 fm). In this section, we rotate the temporal direction using hypercubic symmetry of each gauge configuration and then increase the total number of glueball mass measurements by a factor of four as listed in Table VI. The maximum number of flow iterations corresponds to the diffusion radius R d ≈ 0.5 − 0.6 fm at each ensemble.  Table I. In this sense, the variational analysis [16,17] based on the different shapes is in principle applicable, according to Ref. [15]. However, it is found that the shape-dependence of the resulting two-point functions disappears due to a strong isotropic nature in the spatial gradient flow method as shown in Fig. 8. is given at the different flow time or the different smearing step, to carry out the variational analysis.

B. Variational analysis
As described in the previous subsection, we perform the variational analysis [16,17] with a set of basis operators, which are made of the flowed link variables at the different flow time for a fixed shape k. In this study, we choose the "fish" shaped operator O fish which contains all irreps of our target states (A ++ 1 , For the variational analysis, we construct the N × N correlation matrix of two-point functions of glueball states for given irreps Γ as where the labels of α, β, which run from 1 to N , identify the different flow iterations. The . We next solve the generalized eigenvalue problem, to obtain the nth eigenvalue λ n,Γ (t, t 0 ), where t 0 is a reference time slice, and its eigenvector ω n,β . If only the N lowest states are propagating in the region where t ≥ t 0 , the nth eigenvalue λ n,Γ (t, t 0 ) for n ≤ N is given by a single exponential form with the rest mass of the nth glueball state as which corresponds to the eigenvalue of the transfer matrix between two time slices t and t 0 .
Details of how to practically compute the eigenvalues λ n,Γ (t, t 0 ) are described in Appendix B of Ref. [26]. An effective mass is defined as where λ n,Γ (t, t 0 ) is the nth eigenvalue of the N × N correlation matrix for Γ = A ++ 1 , A −+ 1 , E ++ , T ++ 2 . In this study, we choose N = 6 and the reference time slice as t 0 /a = 0, where the resulting mass is less sensitive to variation of t 0 .
Let us first present the effective masses of glueballs obtained from the variational method using the 6 × 6 correlation matrix constructed by the O fish operator with six different flow iterations. Figure 9 show the effective mass plots of the first two eigenvalues in the A ++

C. Continuum extrapolation and comparison with previous results
It is worth comparing our result obtained from the spatial gradient flow with previous results obtained from both the original gradient flow and a conventional approach. For this purpose, we choose the results obtained in the simulations performed on isotropic lattice with the standard Wilson plaquette action. The previous attempt to apply the gradient flow to the A ++ 1 glueball spectroscopy was done by Chowdhury-Harindranath-Maiti [12] (denoted as CHJ result). Meyer [7] and Athenodorou-Teper [8] (denoted as AT result) performed comprehensive studies of the low-lying glueball states using the conventional approach at several lattice spacings.
In Fig. 13, we show our results of the ground-state glueball masses calculated in the with c 0 being the continuum-extrapolated glueball masses M G (0) in units of r 0 . The fit results using all four data points are shown in Fig. 13 as solid lines. The statistical uncertainty which is estimated by the jackknife analysis is indicated as the gray-shaded band in each panel. The data points are well described by the fitted curves. As mentioned above, our data points calculated at β = 6.2 are slightly overestimated in comparison to the previous works except for the A ++ 1 state, though our continuum results (marked as filled circles) are consistent with the continuum-extrapolated AT results obtained from their high-precision data in all four channels.
Our fit results are compiled in Table VIII. The quoted values of r 0 M G include a systematic error stemming from the continuum-extrapolation fits as the second error, which are evaluated from a difference between the fit results obtained by all four data points and three data points closest to the continuum. Table VIII also includes the previous continuum-limit results for comparison.
In Fig. 13, we finally plot all of data included in Table VIII. The masses are given in terms of the Sommer scale r 0 along the left vertical axis, while the right vertical axis is converted to physical units by r 0 = 0.472(5) fm, 2 which was adopted in Ref. [8].
where the first error is statistical, while the second one represents a systematic error in the continuum extrapolation as explained above.
It is worth recalling that our results and the results of Meyer [7] and AT [8] are obtained from the isotropic lattice simulations, while the results given by Morningstar-Peardon [5] (denoted as MP result) and Chen et al. [6] are obtained from the anisotropic lattice simulations.
The results from the isotropic lattice simulations are systematically underestimated in comparison to those of the anisotropic lattice simulations. This may suggest that there remains some subtlety in taking the continuum limit for the results obtained from the anisotropic lattice simulations. It is beyond the scope of this study, while our aim is rather to demonstrate the feasibility of our proposed approach. Furthermore, as discussed in Appendix B, (bottom-right) channels from this work, Meyer [7], CHJ [12], and AT [8]. Our results are calculated by the spatial gradient flow method. On the other hand, both the Meyer and AT results are given by the conventional approach, while the CHJ result are given by the ordinary gradient flow.
we found that the spatial gradient flow is a few times more effective than the original gradient flow and the conventional approach. We therefore stress that the spatial gradient flow method can efficiently reproduce the recent high precision results of the glueball masses [8] within the statistical uncertainties.

VII. SUMMARY
We have studied the glueball two-point function with two types of the gradient flow method. The original gradient flow, which makes the Wilson flow diffused in the four-  (16) 6.120(52) 6.33 (7) 6.25 (6) dimensional space-time, has some problem in measuring the glueball mass from the twopoint function. It is known to be over smearing due to the overlap of two glueball operators extended in both space-time as reported in the previous study [12]. This particular issue makes the plateau behavior uncertain in the effective mass plot, so that it is difficult to extract the ground-state mass of the glueball with high accuracy.
To avoid over smearing, we propose the spatial gradient flow approach and also apply it to the glueball calculations. Our numerical simulations show that the spatial gradient flow method works well as a noise-reduction technique, meanwhile it has a good property that the plateau behavior in the effective mass plot does not change with variation in flow time for sufficiently large flow time. The latter point gives an advantage for extracting the ground-state mass of the glueball with high accuracy without over smearing.
It is also observed that the spatial gradient flow eliminates dependence of the Wilson loop shapes in the glueball two-point function due to a strong isotropic nature. Therefore, the variational method based on the different shapes is not applicable. Instead, the different diffuseness of the extended operator, which is given at the different flow time, are used for the variational analysis to separate the ground-state contribution from the excited-state contributions.
To demonstrate the feasibility of our proposed method, we have determined the masses of the three lowest-lying glueball states, corresponding to the 0 ++ , 2 ++ , and 0 −+ glueballs, in the continuum limit by using four lattice QCD simulations for the lattice spacings ranging from 0.026 to 0.068 fm. Our results of the 0 ++ , 2 ++ and 0 −+ glueball masses are consistent with the previous works [7,8]. Especially, it is worth emphasizing that the spatial gradient flow method can efficiently reproduce the recent high-precision results [8], which are slightly underestimated in comparison to the results given by the anisotropic lattice simulations [5,6], within the statistical uncertainties.
We have also showed numerical equivalence between the spatial gradient flow and the stout smearing in the glueball calculations at the relatively fine lattice spacing of 0.0513 (3) fm. This observation can be understood through the analytical derivation that demonstrates the equivalence between the gradient flow and stout smearing methods in the vicinity of the continuum limit as described in Appendix A. This fact can help to reflect actual efficiency for the glueball spectroscopy.
As discussed in Appendix B, although the spatial gradient flow requires several times fewer statistics to achieve the same statistical accuracy than the conventional method, its computational cost is roughly a factor of O(10) higher than the conventional one. Thus, by replacing the spatial gradient flow method with the stout smearing method, which is almost as computationally inexpensive as the conventional method, the gradient flow approach becomes really an efficient scheme for the glueball spectroscopy.
We plan to extend our research to calculate the glueball three-point function with the renormalized energy-momentum tensor formulated in the gradient flow method [28] in order to investigate the origin of the glueball masses. Such study is now in progress [29].
expressed with respect to a basis T a as where the operators ∂ a µ,x are defined by When the link derivative ∂ x,µ acts on the action, we may simply focus on the term that explicitly depends on U µ (x) in the action as where X µ (x) is the sum of all the path ordered products of the link variables, called the "staple". The staple X µ (x) is given by If we set Ω µ (x) = X µ (x)U † µ (x), each basis component is given as where Ω † µ (x) denotes the sum of all plaquettes that include U µ (x). Therefore, we finally get which becomes a Lie algebra valued quantity.
In the stout smearing, the link smearing is defined as the following recursive procedure [3].
Here, for simplicity, the stout smearing parameters ρ µν are taken as ρ µν = ρ. The link variables U (k) µ (x) at step k are mapped into the link variables U (k+1) µ (x) using where Q (k) µ (x) is given by Eq. (A8) with the stout link U (k) µ (x) [3]. By taking the logarithm of both sides of Eq. (A9), one can get where ∆ k represents a forward difference with respect to k as ∆ k f (k) ≡ f (k + 1) − f (k).
Next, let us introduce a continuous variable s = kρ, and then reexpress the link variable It is clear that Eq. (A12) is equivalent to Eq. (2) at the leading order. Since the variable s = kρ directly corresponds to the Wilson flow time τ , the perturbative matching relation of τ , ρ and the number of stout smearing steps n st as τ = ρn st found in Ref. [10] is also rigorously proved. Therefore, the gradient flow equation is certainly regarded as a continuous version of the recursive update procedure in the stout smearing at the smaller lattice spacing.
When S W is replaced by S sW in the gradient flow equation, and ρ µν is set as ρ ij = ρ and ρ 4µ = ρ µ4 = 0 in the stout smearing, above consideration fully supports our finding that there is the numerical equivalence between the spatial gradient flow and the spatial stout smearing in the glueball spectroscopy at the relatively fine lattice spacing of 0.0513(3) fm. 4 When U (s) is a matrix Lie group, ∂U ∂s U −1 is given in terms of ∂ ∂s log U as ∂U ∂s U −1 = ∂ ∂s log U + 1 2! log U, ∂ ∂s log U + 1 3! log U, log U, ∂ ∂s log U + · · ·.
In the case when U is the link variable, the power series of log U in the left-hand side can be neglected for the small lattice spacing a.