Glueball spectrum from Nf = 2 lattice QCD study on anisotropic lattices

The lowest-lying glueballs are investigated in lattice QCD using Nf = 2 clover Wilson fermions on anisotropic lattices. We simulate at two different and relatively heavy quark masses, corresponding to physical pion masses of mπ ∼ 938 MeV and 650 MeV. The quark mass dependence of the glueball masses has 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 TrBiDjBk, 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 1 GeV) and compatible with the expected masses of the flavor singlet q q ¯ meson. This indicates that the operator εijk TrBiDjBk and the topological charge density couple rather differently to the glueball states and q 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.


Introduction
Due to the self-interactions among gluons, quantum chromodynamics (QCD) admits the existence of a new type of hadron made up of gluons, usually called the glueball. Glueballs are of great physical interest, 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 the reviews in Refs. [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). However, 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 finds that the tensor meson f 2 (2340) has large branching fractions in the processes J/ψ→γηη [17] and J/ψ→γφφ [18].
Even though quenched lattice QCD studies have provided some information on the existence of glueballs, full lattice QCD studies in the glueball sector are highly desirable. For the masses of the scalar and tensor glueballs, some preliminary unquenched lattice studies have given compatible results [19][20][21][22]. However, for the mass of the pseudoscalar glueball, a consensus has not been reached. For example, in Ref. [21] the authors observed a pseudoscalar glueball state with a mass close to the result in the quenched approximation, but this was 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 necessary 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 gauge configuration ensembles with two different bare quark mass parameters which correspond to the physical pion masses m π ∼ 650 and 938 MeV, respectively. The advantage of using an anisotropic lattice is two-fold: on the one hand, large statistics can be obtained with a relatively low cost in terms 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. First, we will focus on the lowest-lying glueball states, such as the scalar, the tensor and the pseudoscalar states. Secondly, we will then focus more on 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 the 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 2 contains a brief description for the generation of gauge field configurations. Section 3 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 4, 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 5.

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 fourth root of the average spatial plaquette value, incorporates the usual tadpole improvement and γ g designates the bare gauge aspect 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 093103-2 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 to be 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.5 with the bare quark mass parameters m 0 =0.05 and m 0 =0.06. The lattice spacings a s are set by calculating the static potential parameterized as V (r)=V 0 +α/r+σr. Using the Sommer scale parameter r −1 0 = 410 (20) MeV defined through r 2 dV (r) dr | r=r 0 = 1.65, we can determine the ratio where α and σa 2 s are derived from the fit of calculated potential V (r)=V (ra s ) withr being the spatial distance in lattice units. Finally, a s is inverted to the value in physical units by the Sommer's scale parameter r −1 0 =410(20) MeV. The ensemble parameters are listed in Table 1, where we also give the physical values of a −1 t for the two ensembles.  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 2. We use the conventional I = 1 vector and scalar quark bilinear operators as sink operators and the corresponding Gaus-093103-3 sian 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 the two-pion threshold and certainly below the η ′ π threshold. At m π ∼650 MeV, m η ′ is estimated to be m η ′ ∼ 890 MeV (see below in Section 4), thus the mass value of 1362(53) MeV is also below the η ′ π threshold and can be taken as the mass of the 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.

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 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 states 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 refer to Ref. [14] for further details.

Variational method
The continuum SO(3) spatially rotational symmetry is broken into a 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 3 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 , and tensor states with J =2 are reduced to the direct sum of E and T 2 , i.e. (J =2)↓O=E⊕T 2 .
As described in Ref. [14], we use Wilson loops (up to 8 gauge links) as 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 Ref. [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 links, such that different smearing schemes can provide more realizations of the gluonic operators. In practice, we have four different realization of each R by choosing different prototypes. For the smearing of gauge links, we adopt 6 smearing schemes by combining the single-link and double-link smearing procedures with different iteration sequences. Finally we have a set of 24 different gluonic operators, {φ (R) α ,α=1,2,...,24}, for each R P C . Based on these operator sets, we use the variational method to get the optimized operators Φ (R) which mostly project to specific glueball states. In each symmetry 093103-4 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 the A ++ 1 channel. Secondly, we solve the following generalized eigenvalue problem, where 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φ (R) α to build an optimal operator Φ (R) i which overlaps mostly to the i-th state,

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 a one-bin-eliminating jackknife analysis. For the 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) increase 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 functionC(t) with the shifted oneC(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. This subtraction may decrease the signal-to-noise ratio to some extent, since signals can be suppressed (especially for small δt) while the statistical errors may increase. For a too large δt, although signals may remain almost as they were, the subtracting termC(t+δt) can introduce more errors tō C(t). 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 channel) are plotted in Figs. 2, 3, 4 and 5, respectively. Even though the temporal extent of our lattices is T = 128a t , for all the four channels, the signal-to-error ratios of the correlation functions decrease rapidly as t increases and tend to be less than one beyond t>20a t . Therefore, in these figures we only show the effective masses in the time range t/a t ∈[0,20]. In each figure, the left-hand panel shows the result at m π ∼ 938 MeV, and the right-hand 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 full QCD spectrum 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 substantially to other states. 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 a 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 in large statistical errors. Since we focus on the ground states in the present study, in order to get more precise results for 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     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 number of data points). In the fitting procedure, the upper limits t max of the fit windows ofC (R) 1 (t) andC (R) 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 the t max ofC (R) 1 (t) can be larger than 10a t ). Actually, the fit results are insensitive to t max in these ranges since they are almost determined by the data points in the small t range where relative errors are much smaller. For each channel, we keep t max 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 Tables 4 and 5. Except for the t min =1 case in the T ++ 2 channel, all other fits are acceptable with reasonable χ 2 /dof . For all 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 increases 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 Figs .
We calculate these values at each t in the fit windows.
The widths of the bands show the errors estimated through standard error propagation using the covariance error matrix of the parameters, , the a i 's are the six parameters in Eq. (16), and the σ ij 's are elements of the 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. In the large t regions, however, 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 Tables 4 and 5, 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 derived accordingly. This averaging is illustrated in Fig. 6, where the data points are the fitted results 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 6. 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 masses 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 1). From Table 2 and Table 6 one can see that the masses of the ground state scalar meson and our scalar glueball are very close to each other. This may indicate there is mixing between qq and the scalar glueball, which needs further investigation.

Interpretation of the ground states
Generally speaking, the two-point function of an interpolating operator O(t) with definite quantum num-bers is usually parameterized as (18) where {|n ,n = 1,2,...} are eigenstates of a Hamiltonian with eigenvalue m n , which make up an orthogonal, normalized, and complete state set with n |n n|=1, n|n ′ =δ nn ′ .
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 only focus on the physical meaning of the fitted masses of the lowest states. We take the scalar channel as an example. 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 consider only a two-state system composed of the ground state scalar glueball |G and its adjacent state, which could be either |ππ or |f 0 . This then yields the fitting model in Eq. (16) that we introduced previously.
We compare the results in the present study with the previous quenched and unquenched results in Table 7. The tensor glueball masses are obtained by averaging the corresponding E ++ and T ++ 2 values. Despite the fact that glueball correlation functions in the unquenched QCD acquire more complicated spectrum decomposition than the quenched case, the masses 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 Table 7. Comparison of our results with previous results both from quenched lattice QCD studies [13,14] and a full-QCD study [22]. We average the masses of the E ++ and T ++ 2 states to obtain an estimate of the 2 ++ glueball mass. 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 no longer valid.

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 Φ (PS) .
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 ∼ (180MeV) 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 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 the 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 Ref. [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 the 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] to evaluate the topological charge density. Figure 7 shows C q (r) for m π ∼938 MeV and m π ∼ 650 MeV at flow times t = 0. 2   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. (25) to extract the parameter m PS . 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 8 lists the fit ranges, the fitted results of m PS 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. (25) with the fitted parameters. The m PS '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 Φ (PS) . 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 = (180MeV) 4 , f π ∼ 150 MeV for m π ∼ 650 MeV and f π ∼ 200 MeV for m π ∼ 938 MeV, m 2 0 is estimated to be approximately (610MeV) 2 and (460MeV) 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 PS 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 the 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 the O group, c A 1 R are the combinational coefficients corresponding to the A 1 irreducible representation, and W α is any of the four prototypes made up of 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 the 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 9 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 PS = 2.560(140) GeV [13] and 2.590(140) GeV [14]. This is exactly what it should be, since there are only pseudoscalar glueball propa-093103-12 Table 9. Masses of flavor singlet pseudoscalar mesons from 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. 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 flaover 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.

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 TrB 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 states 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 TrB i D j B k and the topological charge density (proportional to TrE · 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, but our current results are still helpful to clarify some aspects of unquenched effects of glueballs and serve as a starting point for further studies.