Glueball spectrum from $N_f=2$ lattice QCD study on anisotropic lattices

The lowest-lying glueballs are investigated in lattice QCD using $N_f=2$ clover Wilson fermion on anisotropic lattices. We simulate at two different and relatively heavy quark masses, corresponding to physical pion mass of $m_\pi\sim 938$ MeV and $650$ MeV. The quark mass dependence of the glueball masses have not been investigated in the present study. Only the gluonic operators built from Wilson loops are utilized in calculating the corresponding correlation functions. In the tensor channel, we obtain the ground state mass to be 2.363(39) GeV and 2.384(67) GeV at $m_\pi\sim 938$ MeV and $650$ MeV, respectively. In the pseudoscalar channel, when using the gluonic operator whose continuum limit has the form of $\epsilon_{ijk}TrB_iD_jB_k$, we obtain the ground state mass to be 2.573(55) GeV and 2.585(65) GeV at the two pion masses. These results are compatible with the corresponding results in the quenched approximation. In contrast, if we use the topological charge density as field operators for the pseudoscalar, the masses of the lowest state are much lighter (around 1GeV) and compatible with the expected masses of the flavor singlet $q\bar{q}$ meson. This indicates that the operator $\epsilon_{ijk}TrB_iD_jB_k$ and the topological charge density couple rather differently to the glueball states and $q\bar{q}$ mesons. The observation of the light flavor singlet pseudoscalar meson can be viewed as the manifestation of effects of dynamical quarks. In the scalar channel, the ground state masses extracted from the correlation functions of gluonic operators are determined to be around 1.4-1.5 GeV, which is close to the ground state masses from the correlation functions of the quark bilinear operators. In all cases, the mixing between glueballs and conventional mesons remains to be further clarified in the future.

The lowest-lying glueballs are investigated in lattice QCD using N f = 2 clover Wilson fermion on anisotropic lattices. We simulate at two different and relatively heavy quark masses, corresponding to physical pion mass of mπ ∼ 938 MeV and 650 MeV. The quark mass dependence of the glueball masses have not been investigated in the present study. Only the gluonic operators built from Wilson loops are utilized in calculating the corresponding correlation functions. In the tensor channel, we obtain the ground state mass to be 2.363 (39) GeV and 2.384(67) GeV at mπ ∼ 938 MeV and 650 MeV, respectively. In the pseudoscalar channel, when using the gluonic operator whose continuum limit has the form of ijk T rBiDjB k , we obtain the ground state mass to be 2.573(55) GeV and 2.585(65) GeV at the two pion masses. These results are compatible with the corresponding results in the quenched approximation. In contrast, if we use the topological charge density as field operators for the pseudoscalar, the masses of the lowest state are much lighter (around 1GeV) and compatible with the expected masses of the flavor singlet qq meson. This indicates that the operator ijk T rBiDjB k and the topological charge density couple rather differently to the glueball states and qq mesons. The observation of the light flavor singlet pseudoscalar meson can be viewed as the manifestation of effects of dynamical quarks. In the scalar channel, the ground state masses extracted from the correlation functions of gluonic operators are determined to be around 1.4-1.5 GeV, which is close to the ground state masses from the correlation functions of the quark bilinear operators. In all cases, the mixing between glueballs and conventional mesons remains to be further clarified in the future.

I. INTRODUCTION
Due to the self-interactions among gluons, Quantum Chromodynamics (QCD) admits the existence of a new type of hadrons made up of gluons, usually called glueballs. Glueballs are of great physical interests since they are distinct from the conventional qq mesons described in the constituent quark model. Glueballs have been intensively studied by lattice QCD and other theoretical methods [1][2][3][4][5][6][7], for more details of this subject, see reviews in [8][9][10][11]. Early lattice QCD studies in the quenched approximation show that the lowest pure gauge glueballs are the scalar, the tensor, and the pseudoscalar glueballs, with masses of 1.5-1.7 GeV, 2.2-2.4 GeV, and 2.6 GeV, respectively [12][13][14].
Experimentally, there are several candidates for the scalar glueball, such as f 0 (1370), f 0 (1500), f 0 (1710), how- * Electronic address: sunw@ihep.ac.cn ever, none of them has been unambiguously identified as a glueball state. On the other hand, J/ψ radiative decays are usually regarded as an ideal hunting ground for glueballs. A few lattice studies have been devoted to the calculation of the radiative production rate of the pure scalar and tensor glueballs in the quenched approximation [15,16]. The predicted production rate of the scalar glueball is consistent with that of f 0 (1710), and supports f 0 (1710) to be either a good candidate for the scalar glueball or dominated by a glueball component. The predicted production rate of the tensor glueball is roughly 1%. It is interesting to note that the BESIII Collaboration find that the tensor meson f 2 (2340) has large branching fractions in the processes J/ψ → γηη [17] and J/ψ → γφφ [18].
Even though the quenched lattice QCD studies have provide some information on the existence of glueballs, it is highly desired that full lattice QCD studies can be performed in the glueball sector. For the masses of the scalar and tensor glueballs, some preliminary unquenched lattice studies have given compatible results [19][20][21][22]. How-arXiv:1702.08174v2 [hep-lat] 7 Sep 2018  ever, for the mass of the pseudoscalar glueball, a consensus has not been reached. For example, in Ref. [21] the authors observed a pseudoscalar glueballs state with a mass close to the result in the quenched approximation, but this is not confirmed by Ref. [22]. On the other hand, owing to the U A (1) anomaly, in the pseudoscalar channel, gluons can couple strongly to the flavor singlet pseudoscalar meson (η in the N f = 2 + 1 case) in the presence of dynamical quarks. Therefore, it is mandatory to identify the contribution of the η meson before one draws any conclusions on the pseudoscalar glueball.
In this work, we attempt to investigate the glueball spectrum using the N f = 2 clover Wilson fermion gauge field configurations that we generated on anisotropic lattices. In order to check the quark mass dependence, we have generated two guage configuration ensembles with two different bare quark mass parameters which correspond to the physical pion masses m π ∼ 600 and 938 MeV, respectively. The advantage of using anisotropic lattice is two-folds: on the one hand, a large statistics can be obtained by a relatively low cost of computational resources, on the other hand, the finer lattice spacing in the temporal direction can provide a better resolution for the signals of the desired physical states. As the first step, we will focus on the lowest-lying glueball states, such as the scalar, the tensor and the pseudoscalar states. Secondly, we will pay more attention to the pseudoscalar channel.
A recent N f = 2+1 lattice study showed that η could be probed by the topological charge density operator [23]. In contrast, a similar study in the quenched approximation found a pseudoscalar with a mass compatible with that in the pure gauge theory [24]. Motivated by this, we use conventional Wilson loop operators to study lowest pseudoscalar glueball state and check for the lowest flavor singlet meson state with topological charge density operator on the same gauge ensembles. This paper is organized as follows: Section II contains a brief description for the generation of gauge field configurations. Section III presents the calculation details and the results of the glueball spectrum. The study of the pseudoscalar channel using the topological charge density operator will be discussed in Section IV, where we will also analyze the difference of the topological charge density operator from the conventional gluonic operators for the pseudoscalar glueball in previous quenched studies. Finally, we will give a summary and an outlook in Section V.

II. LATTICE SETUP
The gauge action we used is the tadpole improved gluonic action on anisotropic lattices [12]: where P ij is the usual plaquette variable and R ij is the 2 × 1 Wilson loop on the lattice. The parameter u s , which we take to be the forth root of the average spatial plaquette value, incorporates the usual tadpole improvement and γ g designates the gauge aspect bare ratio of the anisotropic lattice, denoted as ξ 0 in our former quenched studies [25]. Although γ g suffers only small renormalization with the tadpole improvement [26], we have to tune it by determining the renormalized anisotropy ratio ξ g . As for the tadpole improvement parameter u t for temporal gauge links, we take the approximation u t ≈ 1 following the conventional treatment of the anisotropic lattice setup. We use the Wilson-loop ratios approach, with which the finite volume artifacts mostly cancel [27,28]. We measure the ratios and expect the spatial and temporal behaviors being the same at the correct ξ g . Therefore we find ξ g by minimizing where ∆R s and ∆R t are the statistical errors of R ss and R st . We interpolate R st (x, ξ g y) and its error with a cubic spline interpolation at non-integer ξ g y. Since small x, y may introduce short-range lattice effects and large ones contribute only fluctuations, we scan and test different ranges and finally choose x, y ∈ {2, 3, 4, 5}.
We adopt the anisotropic clover fermion action in the fermion sector [29]: whereF µν = 1 4 Im(P µν (x)) and the dimensionless Wilson operator readŝ The bare fermion aspect ratio γ f is also tuned to make sure that the measured aspect ratio ξ f ≈ ξ g ≈ ξ = 5. ξ f is measured from the dispersion relation of the pseudoscalar and vector mesons where p = 2π k/L s is the momentum on the lattice with periodic spatial boundary conditions. We generate two gauge ensembles on the 12 3 × 128 anisotropic lattice at β = 2.
where α and σa 2 s are derived from the fit to calculated potential V (r) = V (ra s ) withr being the spatial distance in the lattice units. Finally, a s is inverted to the values in the physical units by the Sommer's scale parameter r −1 0 = 410 (20) MeV. The ensemble parameters are listed in Table I, where we also give the physical values of a −1 t for the two ensemble.
The pion masses on the two ensembles are measured to be 938 MeV and 650 MeV respectively. In the following, we use these m π 's to label the gauge ensembles for convenience. Apart from the pion masses, we also calculate the masses of the vector meson and scalar meson for calibration, which are listed in Table II. We use the conventional I = 1 vector and scalar quark bilinear operators as sink operators and the corresponding Gaussian smeared wall source operators to calculate the correlation functions. There is no ambiguity for the vector meson masses m V 's since they are all below the two-pion threshold. For the scalar, we actually deal with a 0 whose two-body strong decay mode is mainly η π (there is only one I = 0 pseudoscalar meson for N f = 2, which is taken as the counterpart of the (approximately) flavor-singlet η in the N f = 3 case). At m π ∼ 938 MeV, the calculated mass in a 0 channel is 1473(28) MeV, which must be the mass of a 0 since it lies below two-pion threshold and certainly below the η π threshold. At m π ∼ 650 MeV, m η is estimated to be m η ∼ 890 MeV (see below in Sec. 4), thus the mass value of 1362(53) MeV is also below the η π threshold and can be taken as the mass of a 0 scalar at this pion mass. In order to calculate the I = 0 scalar meson mass, the disconnected diagrams (quark annihilation diagrams) should be considered. We have not done this yet, but as a rough estimate, we take the a 0 mass as an approximation to the mass of the isoscalar scalar meson.

III. NUMERICAL DETAILS
In this work, the spectrum of the lowest-lying glueballs in three specific channels, namely scalar, tensor and pseudoscalar will be explored. The interpolating operators for these states are pure gluonic operators which have been extensively adopted in the previous quenched lattice studies. In other words, in each specific channel, no operators involving quark fields are included. This of course is only an approximation, assuming that the gluon-dominated state that we are after can be welldescribed by gluonic operators. Needless to say, mixing with the quark operators should be considered later on, especially for cases where the mixing is severe. For completeness, we briefly recapitulate the major ingredients of glueball spectrum computation in the following. One can resort to [14] for further details.

A. Variational method
The continuum SO(3) spatially rotational symmetry is broken into the discrete symmetry described by the octahedral point group O on the lattice, whose irreducible representations R are labeled as A 1 , A 2 , E, T 1 , T 2 , and have dimensions 1, 1, 2, 3, 3 respectively. Therefore, the lattice interpolation fields for a glueball of J P C quantum number should be denoted by R P C with R the irreducible representation of O which may include the components of J in the continuum limit. The parity P = ± and the charge conjugation C = ± can be realized by considering the transformation properties under the spatial reflection and time reversal operations. Since the octahedral group O is a subgroup of SU (2), the subduced representation of SU (2) with respect to O is reducible in general (for integer spin, this occurs for J ≥ 2). Table III shows the reduction of the subduced representation of SU (2) up to J = 5. For instance, the scalar and pseudoscalar with J = 0 states are represented by A 1 , tensor states with J = 2 are reduced to direct sum of E and T 2 , i.e. (J = 2) ↓ O = E T 2 .
As described in [14], we use Wilson loops (up to 8 gauge links) shown in Fig. 1. Each irrep R of group O can be realized by the specific linear combination of its 24 copies of a prototype Wilson loop under the 24 rotation operations of O. The combination coefficients of each R can be found in [14]. So each prototype may provide a different realization of R. On the other hand, the Wilson loops mentioned above can be built from smeared gauge   α , α = 1, 2, . . . , 24}, for each R P C . Based on these operator sets, we use the variational method to get the optimized operators O (R) which mostly project to specific glueball states. In each symmetry channel R, we first calculate the 24 × 24 correlation matrix C (R) (t), In practice, we only apply the vacuum subtraction to the operators in A ++ 1 channel. Secondly, we solve the following generalized eigenvalue problem, where v is the i-th eigenvector, and λ i ≡ e −mi(t0)t0 is the i-th eigenvalue wherem i (t 0 ) is dependent on t 0 and is close to the energy of the i-th state. For all the R channels, we use t 0 = 1. It is expected that the eigenvector v (R) i gives the linearly combinational coefficients of operatorsφ

B. Data analysis
In this work, the correlation function of the optimal operator Φ (R) i for the i-th state is calculated as where we do the summation over the temporal direction to increase the statistics. Accordingly, the effective mass is defined as We divide the measurements into bins with each bin including 100 measurements. The statistical errors are obtained by the one-bin-eliminating jackknife analysis. For A ++ 1 channel, the subtraction of the vacuum is very subtle. Even though we have O(10 4 ) gauge configurations in each ensemble, when we perform the jackknife analysis above after subtracting the vacuum expectation values of the operator, we find there is still a residual (negative) constant term in the correlation function, which makes the effective mass m i, eff (t) going upward when t is large. This problem can be attributed to the large fluctuation of gauge configurations in the presence of sea quarks. To circumvent this difficulty, we adopt a vacuum-subtraction scheme by subtracting the correlation function C(t) with the shifted one C(t + δt), whose spectral expression is where W A ++ 1 ij is the spectral weight of the j-th state iñ Obviously, the possible constant term cancels with the spectrum unchanged. In practice, we take δt = 5a t .
We focus on the R P C = A ++ 1 , A −+ 1 , E ++ , and T ++ 2 channels in this work. For all these channels, the effective masses ofC R Fig. 2, 3, 4 and 5, respectively. In each figure, the left panel shows the result at m π ∼ 938 MeV, and the right panel is for m π ∼ 650 MeV. Even though we have a set of 24 operators for each channel, it is seen that the effective masses do not show plateaus from the very early time slices. This is very different from the case in the quenched approximation. One important reason for this is that, in each channel, the spectrum of the full QCD is much more complicated than in the quenched approximation due to the sea quarks. This is true in principle, since qq states and multi-hadron states with the same quantum number do contribute to the corresponding correlation function in the presence of sea quarks.
Given the limited number of independent operators, our optimal operator Φ (R) i is actually not optimized as expected, namely, it does not only overlap to the i-th state but also to other states substantially. As seen in the effective mass plots, when m (R) 1,eff (t) tends to reach a plateau as t increases, m (R) 2,eff (t) decreases gradually and finally merges into this plateau at large t (within errors). Even though one can carry out the single exponential fit to the mass of the ground state in the plateau range roughly beyond t/a t ≈ 6 or 7, the bad signal-to-noise ratio in this time range results with large statistical errors. Since we focus on the ground states in the present study, in order to get more precise result of the masses of the ground states, we adopt the following data-analysis strategy which also makes use of the measured data in the short time range. In each channel, we carry out a correlated fit toC 2 (t) simultaneously through the following function forms, where the second mass term is introduced to take into account the contribution of the second state and higher states (of course, one can add more mass terms, but more parameters will ruin the data fitting due to the limited data points). In the fitting procedure, the upper limit t max 's of the fit windows ofC  2 (t) are chosen properly to include only the data points with good signal-to-noise ratios (The t max ofC (R) 2 (t) are set to be from 7a t to 9a t , while t max ofC (R) 1 (t) can be larger than 10a t ). Actually, the fit results are insensitive to t max 's in these ranges since they are almost determined by the data points in small t range where relative errors are much smaller. For each channel, we keep t max 's fixed and vary t min to check the stability and the quality of the fit. The fit results for the scalar (A ++ 1 ), the pseudoscalar (A −+ 1 ) and the tensor channels ( E ++ and T ++ 2 ) at the two pion masses are listed in Table IV and V. Except for t min = 1 case in T ++ 2 channel, all other fits are acceptable with reasonable χ 2 /d.o.f . For all the four channels, the fitted parameters m 1 and W 11 are stable with respect to the various t min , while m 2 decreases as t min increasing gradually. This signals that our fitting model in Eq. (16) is not so good that we should include more mass terms to account for higher states, which, however, affect the second state more than the first state. Since we are interested only in the first states, we do not take m 2 seriously and treat it as an object accommodating the effect of higher states.
In Fig. 2, 3, 4 and 5, we also plot the shaded bands to illustrate the goodness of the fits. For each channel, after the correlated fit to the two correlations simultaneously, we get the six parameters m 1 , m 2 , W  Table IV and V. The red and blue bands are obtained through the function . (17) We calculate these values at each t in the fit windows. The widths of the bands show the errors estimated through the standard error propagation using the covariance error matrix of the parameters, i,eff (t), a i 's are the six parameters in Eq. (16) and σ ij 's are elements of covariance error matrix of the parameters, which are obtained directly from the fit. The extensions of the red and blue bands corresponds to the actual fit windows.
It is seen that the fit model describes the data of the ground state very well throughout the fit windows. For the second states, the fit model also fits the data more or less, especially in the small t region. While in the large t regions, the fitted results deviate somewhat from the data. This is understandable, since higher states, which do contribute, are missed in this model. This deviation actually contributes much to the χ 2 . It is expected that the fitted m 2 is generally (much) higher than the mass of the second state.     As shown in Table IV and V, most of the fits using different t min are statistically acceptable and the masses of the first states are relatively stable. Therefore, for the final result of m 1 in each channel, we take tentatively the average value of m 1 's at different t min weighted by their inversed squared errors. The statistical errors are accordingly derived. This averaging is illustrated in Fig. 6, where data points are the fitted result of m 1 at different t min and the shaded bands are the averaged values with averaged errors. The results are also listed in Table VI. At the heavy pion mass m π ∼ 938 MeV, m 1 (E ++ ) is very close to m 1 (T ++ 2 ), as expected by the rotational symmetry restoration in the continuum limit where they correspond to the mass of the same 2 ++ tensor state. However for the lighter m π ∼ 650 MeV, the two masses deviate from each other by 200 MeV. Since the lattice spacings at the two pion massed are very close, the extent of the rotational symmetry breaking should be similar. We tentatively attribute this large deviation to the relatively small statistics at m π ∼ 650 MeV, which is roughly one-half as large as that at m π ∼ 938 MeV (see Table I). From Ta-ble II and Table VI one can see that the masses of ground state scalar meson and our scalar glueball are very close to each other, this may indicate there are mixing between qq and the scalar glueball, which needs further investigation.

C. Interpretation of the ground states
Generally speaking, the two-point function of an interpolating operator O(t) with definite quantum numbers is usually parameterized as        For QCD on a Euclidean spacetime lattice, m n take discretized values and the connection of these discretized energy levels to the relevant S-matrix parameters should be established through other theoretical formalisms, such as Lüscher's. Here we would only focus on the physical meaning of the fitted masses of the lowest states.  We take the scalar channel for instance. A hadron system of the bare states with the scalar quantum number J P C = 0 ++ can be a bare scalar glueball |G 0 ++ , a bare qq scalar meson |f 0 , or even ππ scattering states |ππ . We simplify the matter further by assuming that the two adjacent states mix most, then we can only consider a two-state system composed of the ground state scalar glueball |G and its adjacent state, which could be of nature |ππ or |f 0 . This then yields the fitting model in Eq. (16) that we introduced previously. values. Despite the fact that glueball correlation functions in the unquenched QCD acquire more complicated spectrum decomposition than the quenched case, the mass of the bare glueball states |G can still be obtained by assuming the corresponding operators O couple weakly to other states. Therefore, it is naturally understood that the glueball spectrum in our full-QCD lattice studies is similar to that in the quenched approximation. The difference is still visible, however, and it is most evident in the scalar channel where one would expect that this weak coupling assumption is not valid anymore.

IV. FURTHER STUDY ON THE PSEUDOSCALAR CHANNEL
As presented in the last section, in the A −+ 1 channel, we obtain the mass of the ground state to be m A −+ 1 ∼ 2.6 GeV at the two pion masses, which is compatible with the pure gauge glueball mass. Theoretically, in the presence of sea quarks, the flavor singlet qq pseudoscalar meson is expected to exist, but we do not observe this state from the correlation function of the glueball operator Φ (P S) .
In order to check the existence of the flavor singlet pseudoscalar meson in the spectrum, we would like to study the correlation function of topological charge density operator q(x). This is motivated by the partially conserved axial current (PCAC), where g is the strong coupling constant, P (x) = ψ(x)γ 5 ψ(x) is the pseudoscalar density, and the anomalous gluonic operator µνρσ F µν F ρσ is the so-called topological charge density (up to a constant factor), which is usually denoted by q(x). Thus q(x) may have substantial overlap with the flavor singlet pseudoscalar meson (denoted by η ). The correlation function of q(x) is expressed as from which one can get the topological susceptibility where V 4 is the four-dimensional volume of the Euclidean spacetime. It is known that χ t is positive and takes a value ∼ (180 M eV ) 4 . On the other hand, q(x) is a pseudoscalar operator and requires C q (x − y) < 0 for r = ||x − y|| > 0. So C q (x − y) can be intuitively expressed as whereC q (x − y) is negative for r > 0. On the Euclidean spacetime lattice with a finite lattice spacing, the delta function will show up a positive kernel with a width of a few lattice spacings, and C q (x − y) has a negative tail contributed fromC q (x − y). It is expected thatC q (x − y) would be dominated by the contribution of the lowest pseudoscalar meson in the large r range and can be parameterized as [30] C q (r) = N m PS 4π 2 r K 1 (m PS r), where N is an irrelevant normalization factor, m PS is the mass of the lowest pseudoscalar, and K(z) is the modified Bessel function of second kind, whose asymptotic form at large |z| is Therefore, one can obtain m PS by fitting the negative tail of C q (x − y) in the large r range using the above functional form. This has been actually done by several lattice studies in both the quenched approximation [24] and full QCD calculations [23]. In the quenched approximation, the extracted m PS = 2563(34) MeV is in good agreement with the pseudoscalar glueball mass m PS = 2560 (35) MeV. This is as it should be, since the hadronic excitations of a pure gauge theory are only glueballs. In the full-QCD study with N f = 2 + 1 and pion masses close to the physical m π , m PS is obtained to be 1013(117) MeV, which is consistent with the mass of the physical η . In this work, we adopt a similar strategy to that in [23]. The topological charge density q(x) is defined by the spatial and temporal Wilson loops (plaquettes) as conventionally done. We use the Wilson gradient flow method as a smearing scheme to optimize the behavior of topological charge density correlator [23,31]. The Wilson flow provides a reference energy scale 1 √ 8t [32]. In practice, we use the code published by the BMW collaboration [33] [13,14] and the full-QCD study [22]. We average the masses of E ++ and T ++ 2 states to obtain the estimate of the 2 ++ glueball mass.  to evaluate the topological charge density. Fig. 7 shows C q (r) for m π ∼ 938 MeV and m π ∼ 650 MeV at flow times t = 0.2, 0.3, 0.4, 0.8 respectively. On our lattices, these t values correspond to √ 8t ∼ 0.15, 0.18, 0.21 and 0.30 fm. As shown in the figures, at large flow time, C q (r) is mostly positive, which implies that the gauge fields are over smeared.
In order to compare the large r behaviors of C q (r) at different flow times, we plot them in Fig. 8 in logarithmic scale, where one can see that their behaviors are similar in the large r region, but the C q (r) at t = 0.4 looks the smoothest and has the smaller errors. Therefore, we fit the C q (r) at t = 0.4 directly through the function form of Eq. 24 to extract the parameter m P S . In determining the fit range, we take the following two factors into consideration. First, the spatial extension of our lattices is L s = 12a s . In order to avoid large finite volume effects, the upper limit of the fit range should be smaller than 6a s , due to the periodic spatial boundary condition. Secondly, as shown in Fig. 7, the negative tail of C q (r) starts beyond r ∼ 3a s , which requires the lower limit of the fit range to be larger than 3a s . In the practical fitting procedure of C q (r) at t = 0.4, we choose the fit range to be r/a s ∈ [3.8, 5.4].
We carry out a correlated minimal-χ 2 fit to C q (r) at t = 0.4 in the r interval described above. Table VIII lists the fit ranges, the fitted results of m P S and the χ 2 /dof 's at the two pion masses. In order to illustrate the fit quality, we also plot C q (r) in Fig. 8 in red bands using the function form in Eq. 24 with the fitted parameters. The m P S 's we get are around 1 GeV and show explicit dependence on the pion mass. However, they are much smaller than the values around 2.6 GeV from the correlation functions of the pseudoscalar glueball operator Φ P S . Thus the light pseudoscalar state observed in C q (r) can be naturally assigned to be the flavor singlet qq state η . Theoretically, the mass of η is acquired through the interaction of sea quark loops according to the Witten-Veneziano mechanism [34,35]. In this mechanism, the propagator of η can be expressed as (26) where the parameter m 2 0 is introduced to describe the gluonic coupling, such that On the other hand, m 2 0 is related to the topological susceptibility χ t through   where f π is the decay constant of π. For our case of N f = 2, if we take the values χ t = (180 MeV) 4 , f π ∼ 150 MeV for m π ∼ 650 MeV and f π ∼ 200 MeV for m π ∼ 938 MeV, m 2 0 is estimated to be approximately (610 MeV) 2 and (460 MeV) 2 , respectively. Thus the η mass can be derived as m η ∼ 890 MeV for m π ∼ 650 MeV, and m η ∼ 1045 MeV for m π ∼ 938 MeV. These values are not far from the m PS 's we obtained.
Because these are very preliminary calculations and the systematic errors are not well under control, we do not want to overclaim the values of m P S we obtain. What we would like to emphasize is that there does exist in the spectrum a flavor singlet qq pseudoscalar meson corre-sponding to the η meson in the real world, which can be accessed by the topological charge density operator. Now that the η state exists in the spectrum, there comes a question of why it is missing in the correlation function of the conventional gluonic operator for the pseudoscalar glueball (denoted as Φ PS ). In order to clarify this, we check the continuum form of Φ (PS) involved in this work. Actually, in the construction of the gluonic pseudoscalar operators, only the spatially solid (instead of planar) Wilson loops (the last four prototypes in Fig. 1) are used, where R stands for each rotation operation in O group, c A1 R is the combinational coefficients corresponding to the A 1 irreducible representation, W α is any of the four prototypes made up of a specifically smeared gauge links. According to the non-Abelian Stokes theorem [36], a rectangle Wilson loop P a×b µν (x) of size a × b, with a, b small, can be expanded as where F µν is the strength tensor of the gauge field. For simplicity, the factor ig is absorbed into the quantity F µν .
The small ab expansion of P ±µ±ν (x) is similar to Eq. 30 by replacing a and b with ±a and ±b, respectively. Since the last four prototypes can be expressed as products of two rectangle Wilson loops, using the above relation one can obtain the leading term of the pseudoscalar operator, which is obviously different from the anomalous part of the PCAC relation, µνρσ F µν (x)F ρσ (x) ∝ E(x) · B(x). Actually, the operator Φ (PS) is a linear combination of these kinds of operators defined through differently smeared gauge fields. This may imply that the two operators couple differently to specific states. Along with the observation in the calculation of glueball spectrum, this proves to some extent that our operator for the pseudoscalar glueball couples very weakly to the qq meson state and almost exclusively to the glueball states.
We collect the existing lattice results of the masses of flavor singlet pseudoscalar mesons in Table IX for an overview. In the quenched approximation (N f = 0), the authors of Ref. [24] use q(x) as pseudoscalar operators and derive the ground state mass m PS = 2.563(34) GeV, which is almost the same as the mass of the pure gauge pseudoscalar glueball m P S = 2.560(140) GeV [13] and 2.590(140) GeV [14]. This is exactly what it should be, since there are only pseudoscalar glueball propogating along time if no valence quarks are involved. When dynamical quarks are included in the lattice simulation, the situation is totally different. There have been several works using P (x) to calculate the η mass in the lattice simulation with dynamical quarks, and have given the results m η = 768 (24) MeV (N f = 2) [37], m η = 947(142) MeV (N f = 2 + 1) [38] and m η = 1006(65) MeV (N f = 2 + 1 + 1) [39], which almost reproduce the experimental result m η = 958 MeV. When the q(x) operator is applied, N f = 2 + 1 lattice simulation gives the result m η = 1019(119) MeV at the physical pion mass [23], which is consistent with the result through the P (x) operator. We also calculate the ground state mass using the q(x) operator on our N f = 2 gauge configurations and obtain the result m PS = 890 (38) MeV at m π = 650 MeV, which is compatible with the m η = 768(24) MeV above (note that our m π is higher than that in Ref. [37]). The similar result for m η from the operators P (x) and q(x) can be understood as follows. Due to the U A (1) anomaly, q(x) is now related to P (x) through the PCAC relation. The relation implies that q(x) can couple substantially to the flavor singlet η meson. In contrast, the glueball operator Φ (PS) couples predominantly to the pseudoscalar glueball state either in the quenched approximation or in the presence of sea quarks.

V. SUMMARY AND CONCLUSIONS
The spectrum of the lowest-lying glueballs is investigated in lattice QCD with two flavors of degenerate Wilson clover-improved quarks. We generate ensembles of gauge configurations on anisotropic lattices at two pion masses, m π ∼ 650 MeV and m π ∼ 938 MeV. Focus has been put on the ground states of the scalar, pseudoscalar and tensor glueballs, which are measured using gluonic operators constructed from different prototypes of Wilson loops. The variational method is applied to obtain the optimal operators which couple dominantly to the ground state glueballs.
In the tensor channel, we obtain the ground state mass to be 2.363(39) GeV and 2.384(67) GeV at m π ∼ 938 MeV and 650 MeV, respectively. In the pseudoscalar channel, using the gluonic operator whose continuum limit has the form of ijk T rB i D j B k , the ground state mass is found to be 2.573(55) GeV and 2.585(65) GeV at the two pion masses. The masses of the tensor and pseudoscalar glueballs do not show strong sea quark mass dependence in our study. However, since our pion masses are still heavy, no decisive conclusions can be drawn on the quark mass dependence of glueball masses at present. In the scalar channel, the ground state masses extracted from the correlation functions of gluonic operators are determined to be around 1.4-1.5 GeV, which is close to the ground state masses from the correlation functions of the quark bilinear operators. One possible reason is the mixing between glueball states and conventional flavor singlet mesons, which requires further investigation in the future.
We also investigate the pseudoscalar channel using the topological charge density as the interpolation field operator, which is defined through Wilson loops and smeared by the Wilson flow technique. The masses of the lowest state derived in this way are much lighter (around 1 GeV) and compatible with the expected masses of the flavor singlet qq meson. This provides a strong hint that the operator ijk T rB i D j B k and the topological charge density (proportional to T rE · B) couple rather differently to the glueball states and qq mesons.
Admittedly the lattice volumes we used are relatively small and the continuum limit remains to be taken, our current results are still helpful to clarify some aspects of unquenched effects of glueballs and serves as a starting point for further studies.  IX: The table collects the masses of flavor singlet pseudoscalar mesons from the quenched and unquenched lattice QCD studies. P (x), q(x) and Φ (PS) stand for the quark bilinear pseuscalar operator, the topological charge density, and the pseudoscalar glueball operator, respectively.