Glueballs at Physical Pion Mass

We study glueballs on two $N_f=2+1$ RBC/UKQCD gauge ensembles with physical quark masses at two lattice spacings. The statistical uncertainties of the glueball correlation functions are considerably reduced through the cluster decomposition error reduction (CDER) method. The Bethe-Salpeter wave functions are obtained for the scalar, tensor and pseudoscalar glueballs by using spatially extended glueball operators defined through the gauge potential $A_\mu(x)$ in the Coulomb gauge. These wave functions show similar features of non-relativistic two-gluon systems, and they are used to optimize the signals of the related correlation functions at the early time regions. Consequently, the ground state masses can be extracted precisely. To the extent that the excited state contamination is not important, our calculation gives glueball masses at the physical pion mass for the first time.

Introduction In quantum chromodynamics (QCD), gluons carry color charges and have the direct strong interaction among themselves. Therefore it is expected that hadrons can be made up of gluons solely, namely glueballs. The property of glueballs has been one of the hottest topics in hadron physics for more than several decades and have aroused extensive and intensive experimental and theoretical studies. However, their existence has not been confirmed yet.
Glueballs are well defined objects in the quenched approximation where the quark-gluon transitions are switched off, and the quenched lattice QCD (QLQCD) studies have confirmed the existence of purge gauge glueballs and have predicted their spectrum [1,2]. In the presence of dynamical quarks, the situation is more complicated owing to the decays of glueballs and the possible mixing between glueballs and qq mesons. There have been a few preliminary full-QCD lattice studies [3][4][5] at pion masses much larger than the physical value, which have observed possible glueball states with masses close to the predictions from QLQCD. However, it is still challenging to verify the previous glueball studies at the physical pion mass due to the computational cost, since glueball relevant lattice studies usually require large statistics of thousands gauge configurations, which is computationally prohibitive for physical pion masses and large spatial volumes in the present stage. Fortunately, the cluster decomposition principle ensures that the correlation length between the glueball operators is insensitive to the volume, and can be implemented  [9]. to reduce the statistical errors of glueball correlation functions. In Ref. [6][7][8], a cluster decomposition error reduction (CDER) method was introduced to suppress the statistics requirement on lattices with large spatial volumes. In this work, we use the CDER method to trade the statistical uncertainty for negligible systematic one, and obtain clear mass spectrum and Bethe-Salpeter wave function with the pure glueball operators on only ∼ 300 configurations at physical pion mass. Numerical details We choose two RBC/UKQCD gauge ensembles (labeled as 48I and 64I) with 2 + 1 flavor domain wall fermion at physical pion and kaon masses and with spatial sizes around 5.5 fm. The parameters of the ensembles are shown in Table I of reference [9]. In order to extract glueball states, we adopt the same strategy in Ref. [1,2] to build the glueball operator set S(R P C ) = {O α , α = 1, 2, . . . , 24} for the scalar (R P C = A ++ 1 ), pseudoscalar (A −+ 1 ), and tensor (E ++ ⊕ T ++ 2 ) channels, where R = A 1 , E, T 2 are the irreducible representations of the lattice symmetry group O, P and C are the parity and the charge conjugate quantum numbers, respectively.
It is well known that the glueball relevant lattice study usually requires a large statistics, but we have only N conf ∼ 300 configurations available for the two ensembles. Fortunately, its lattice size is large enough to make the cluster decomposition error reduction method (CDER) proposed in Ref. [7] to be efficient on improving the signal-to-noise ratio in the calculation. The idea of CDER is the following: since C αβ (r, t) ≡ 0|O α (r, t)O β (0, 0)|0 behaves as ∼ ξ − 3 2 e −mξ for large ξ, where m is the lowest mass in this channel and ξ = (r 2 + t 2 ) 1/2 is the Euclidean separation between the sink and the source operator, the signals of C αβ will be undermined by statistical noises when |r| is beyond some length scale |r c | ∝ 1/m. Therefore the correlation matrix C αβ with the operator set S(R P C ) in a given R P C channel can be calculated as where the average over time slices is taken into account to improve the statistics and the kernel functions K αβ (r, t) is introduced as in terms of the Fourier transformed operatorŝ The cutting scale parameter r c is chosen empirically when the value of C αβ (r c , t) saturates. The efficacy of CDER is illustrated in Fig. 1, where a typical C αβ (r c , t) in A ++ 1 channel at t/a = 0, 1, 2, 3 is plotted versus r c /a, where we use the standard deviations instead of the statistical ones such that the varying of errors with respect to r c to be seen more clearly. One can see that the central values of C αβ (r c , t) at different t saturate uniformly beyond r c /a = 7 but the errors grow when increasing r c . In the practical calculation we choose r c /a = 7 for all the C αβ (R, t) in all the channels.
In each R P C channel, after the correlation matrix C αβ (t) is estimated through Eq. (1), we perform the standard variational method by solving the generalized α O α that is expected to couple most with the n-th lowest state |n contributing to C αβ (t). In practice, the eigenvectors v (n) α are normalized by C n (t = 0) = 0|O n (0)O n (0)|0 = 1, which implies 0|O n |m ≈ δ mn if the states are normalized as m|n = δ mn .
Although the operator basis in each channel is large, it is found that the correlation function of the operator O 1 optimized for the ground state has still a sizeable contribution from higher states, (the effective mass function does not show a plateau good enough in the early time range, see correlation function at first 4 times lices against the cutoff on 48 3 × 64 lattice. The rc is in lattice unit. As expected, the signal reaches maximum at r ∼ 7a and only noise is increasing afterwards. Here the errorbars indicate the standard deviations instead of statistic errors.

materials [10]
). In order to enhance the contribution of the ground state to the corresponding correlation functions, we construct another type of gluonic operators in terms of the gauge potential A µ (x) [11,12] in the Coulomb gauge. According to the definition of the gauge link oriented in the µ direction starting from the lattice site, i.e, U µ (x) = exp(−igA µ (x + aμ/2)), up to an irrelevant pre-factor we have where lnU µ (x) can be derived by lnU being the eigenvalues of U µ (x) and V (x) being the unitary matrix that diagonalizes U µ (x). On each time slice, the spatially extended operator O R A (r, t) can be constructed by A µ (x) at two different points.
where N r is the degeneracy of r with the same r and S R ij are the combinational coefficients related to the irreducible representation R of the spatial symmetry group O (octahedral group) on the lattice. For R = A 1 , the non-zero values of S R ij are S R 11 = S R 22 = S R 33 = 1. The E representation has two components, the first of which is given by S R 11 = −S R 22 = 1/ √ 2, and the other is given by wherer is the orientation vector of r.
In each channel, by the implementation of CDER similar to Eq. (1) and (2), we calculate the correlation function C An (r, t) of O R A and the n-th optimized operator O n , where the last parameterization is due to predominant coupling of O n to the n-th state with i 1 and Φ n (r) ∝ 0|O R A (0, t; r)|n is interpreted as the Bethe-Salpeter (BS) wave function [13] of the n-th state |n if it is a two-gluon glueball state. In principle, Φ n (r) can be approximated by C An (r, t) at any t according to Eq. (6). However, due to the rapid increase of the noises, we can only observe clear signals of Φ n (r) at the first few time slices. Figure 2 shows the profiles of Φ 1,2 (r) at t = 0 for A ++ 1 , A −+ 1 , E ++ and T ++ 2 channels on 48I and 64I lattices (normalized by Φ(r = 0) = 1 for P C = ++ channels and Φ(r = a) for the A −+ 1 channel), where the values of r are converted to physical units using the lattice spacings a in Table I. The data points are the calculated values and the curves are plotted using the fitted parameters through some phenomenological functional forms (see the Supplemental materials [10]). It is interesting that the r-behaviors of Φ 1,2 (r) in P C = ++ channels are very similar to the radial 1S and 2S Schrödinger wave function of a two-body system with a central potential, while those in the A −+ 1 channel have the feature of 1P and 2P wave functions. Note that by solving the Bethe-Salpeter equation of the 0 −+ glueball, the two-body Pwave feature of the 0 −+ Bethe-Salpeter amplitude was also observed in Ref. [14].
We would not like to put too much emphasis on the interpretation of the wave functions Φ n (r), but will use the feature of the r-behaviors of these functions to enhance the contribution of the ground states in the early time region. Since Φ 2 (r) in each channel has a radial node at r = r 1 , the correlation function C A1 (r 1 , t) in each channel is expected to be dominated by the ground state, since the second term in Eq. (6) is further suppressed by Φ 2 (r 1 ) ≈ 0. According to Fig. 2, the radial node of Φ 2 (r) appears at r ∼ 0.12, 0.19 and 0.22 fm in the scalar (A ++ 1 ), the tensor (E ++ and T ++ 2 ) and the pseudoscalar channel (A −+ 1 ), respectively. According to the lattice spacings in Table I, we choose r 1 (A ++ 1 )/a = 1, r 1 (E ++ , T ++ 2 )/a = √ 3 and r 1 (A −+ 1 )/a = 2 on the 48I lattice. On the 64I lattice, these r 1 /a's are chosen to be √ 3, At these values of r 1 , Φ 2 (r) are seen to be approximately zero. The effective mass of each C An (r 1 , t) for all the four channels on the two ensembles (48I and 64I) are illustrated in Fig. 3, where the vertical axis and the horizontal axis are plotted in physical units. It is seen that the effective masses show more or less plateaus at first several time slices, which indicate that C A1 (r 1 , t) can be described with a single exponential as expected. With the prescription mentioned above, the ground state masses in each channels are extracted through single exponential function forms in early time windows. In each panel of Fig. 3, the grey bands illustrate the fitted results and the time interval of the fits. In order to exhibit the efficacy of the CDER method, we also show the lattice data that are derived by calculating the correlation functions without the implementation of CDER in green (48I ensemble) and yellow bands (64I ensemble). Obviously, the data without CDER fluctuate drastically with respect to time and the errors are quite large. In contrast, the data points based on CDER method have much smaller errors and can be described by the single exponential function forms. It is impressive that the effective mass plateau of the A ++ 1 ground state lasts over 5 time slices. The fitted mass values are listed in Table II  and A −+ 1 channels, the ground state masses on 48I and 64I ensembles are consistent with each other within errors. Since the two ensembles have similar physical volumes but different lattice spacings, this consistence imply that the finite lattice spacing a effects are not large.
• The masses of ground states in E ++ and T ++ 2 channels are compatible with each other on the same ensemble. This also indicates that the discretization effects are not important since these states correspond to the 2 ++ state in the continuum limit.
• If the ground states observed are free from excited states, their masses are slightly higher than but not far from those of the quenched approximation predictions and preliminary full-QCD calculations [1-

5].
It should be noted that, the situation for glueballs is much complicated in the presence of dynamical quarks. Experimentally, most of hadrons are observed as resonances, therefore their decays should be taken into account. The energy levels obtained on the lattice are actually the eigen energies of the lattice Hamiltonian, and their connection with the resonances in the real world is highly nontrivial. In full QCD, glueballs can mix with qq and multi-hadron states. However, with the gluonic operators, theqq states are suppressed by O(1/N c ) and two meson states are suppressed by O(1/N 2 c ). This is perhaps the reason that we could see the effective mass palteaus in short time saparation. To have a complete picture, it is essential to include quark operators in the calculation, which is beyond the scope of the present work. There is an exploratory lattice study in this direction in the scalar channel [15], but no glueball states can be identifiable unambiguously yet.
Summary An exploratory study of glueballs is performed on two N f = 2 + 1 gauge ensembles with large lattice sizes and the dynamical quark masses being tuned at the physical point. The large lattice size enables us to use the CDER technique to reduce the errors of the correlation functions of glueball operators. By using the spatially extended glueball operators defined through the gauge potential A µ (x) in the Coulomb gauge, the Bethe-Salpeter wave functions of the scalar (A ++ 1 ), the tensor (E ++ 2 ⊕ T ++ 2 ) and the pseudoscalar (A −+ 1 ) states are derived, which are similar to the features of the non-relativistic wave functions of two-body systems. Even though the physical meaning of these BS wave functions may not be clear, we can make use of their radial behaviors to further improve the ground state domination of the related correlation functions at the early time region, where the ground state masses in each channel can be extracted through single-exponential function form. To the extent that the excited state contamination is small and the couplings toqq meson states are negligible (they could show up at much large time separation given enough statistics). Our calculation with physical dynamical quarks turns out to be compatible with the quenched lattice results. A scrutinized full-QCD study should be carried out by including theqq and multi-meson operators, as well as considering the resonance nature of most hadrons. As a noise reduction scheme, the same CDER technique can be applied to all the correlation functions involving the disconnected quark insertions. In this sense, our study indicates a breakthrough to improve the signal-to-noise ratios of correlation functions involving glueball operators, and pave the way to the final answer on the existence of the glueball states.  should be dominated by the contribution from the ground state. After introducing the effective mass function

ACKNOWLEDGMENTS
the ground state dominance can be monitored by the temporal behavior of m eff (t). Figure 4 shows m eff (t) of A ++ 1 , A −+ 1 , E ++ and T ++ 2 channels on ensemble 48I. It is seen that for each of these channels, m eff (t) has good signals but appreciable time dependence at the first several time slices. This indicates that the contamination from higher states has not been suppressed drastically even with the optimized operator O 1 .

S2. Parameterization of BS functions
The BS wave functions in Fig. 2 exhibit similar features to the radial wave functions of a two-body system with a central potential, so we try to parameterize them tentatively using phenomenological functions forms. For the scalar (A ++ 1 ) and the tensor (E ++ and T ++ 2 ) channels, the BS wave functions of the ground states and the first excited states, namely Φ 1 (r) and Φ 2 (r), are very similar to the 1S and 2S wave functions, therefore we describe them using the following function forms [12] Φ 1 (r) = Φ 1 (0)e −(r/r0) α + C Φ 2 (r) = Φ 2 (0)(1 + βr α )e −(r/r0) α + C where the constant terms C and C account for the constant tails in the large r region, which may be due to the unphysical modes of the domain-wall sea quarks [16]. For the pseudoscalar A −+ 1 channel, Φ 1,2 (r) show the typical feature Φ 1,2 (0) = 0 of P -wave radial wave functions. Therefore, the functions can be parameterized as Φ 1 (r) = A r e −(r/r0) α + D Φ 2 (r) = A r (1 − βr α )e −(r/r0) α + D We perform naive fits to Φ 1,2 (r) in each channel using the corresponding function forms mentioned above The fitted parameter values are shown in Table III. It should be noted that the lineshapes and the parameters of BS wave