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

TL;DR

This study computes the masses of the lowest-lying glueballs in $N_f=2$ lattice QCD with anisotropic lattices, revealing operator-dependent differences and potential effects of dynamical quarks on pseudoscalar states.

Contribution

First lattice QCD calculation of glueball spectrum with dynamical quarks using anisotropic lattices and gluonic operators, highlighting operator dependence and quark mass effects.

Findings

Tensor glueball mass around 2.36-2.38 GeV

Pseudoscalar glueball mass around 2.57-2.59 GeV with gluonic operators

Light pseudoscalar states (~1 GeV) when using topological charge density

Abstract

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.

Tables

Table 1. Table 1: Parameters of configurations. The spatial lattice spacing a s subscript 𝑎 𝑠 a_{s} is set by the calculation of the static potentials and the Sommer’s scale parameter r 0 − 1 = 410 ​ ( 20 ) superscript subscript 𝑟 0 1 410 20 r_{0}^{-1}=410(20) MeV. We also give the value of a t − 1 superscript subscript 𝑎 𝑡 1 a_{t}^{-1} in the physical units.

$\beta$	$m_0$	$L^3\times T$	$\xi$	$a_s$(fm)	$a_t^{-1}$ (GeV)	$N_{conf}$
$2.5$	0.05	${12}^3\times 128$	5	0.114(1)	8.654(76)	4800
$2.5$	0.06	${12}^3\times 128$	5	0.118(1)	8.360(76)	10400

Table 2. Table 2: The masses of the ground state pseudoscalar, vector and scalar mesons (these are actually isovector mesons since we ignore the disconnected contributions). The measured values are also inverted to the values in physical units through a t − 1 superscript subscript 𝑎 𝑡 1 a_{t}^{-1} in Table  1

$m_{PS}{a}_t$	$m_{PS}\left(\mathrm{MeV}\right)$	$m_V{a}_t$	$m_V\left(\mathrm{MeV}\right)$	$m_S{a}_t$	$m_S\left(\mathrm{MeV}\right)$
$0.07508\left(50\right)$	$650\left(4\right)$	$0.1147\left(19\right)$	$993\left(16\right)$	$0.1574\left(61\right)$	$1362\left(53\right)$
$0.1119\left(4\right)$	$938\left(3\right)$	$0.1388\left(12\right)$	$1164\left(10\right)$	$0.1757\left(34\right)$	$1473\left(28\right)$

Table 3. Table 3: Reduction of subduced representation of S ​ U ​ ( 2 ) 𝑆 𝑈 2 SU(2) with respect to octahedral group O 𝑂 O up to J = 5 𝐽 5 J=5

	$0$		$1$		$2$		$3$		$4$		$5$
$A_1$	1		0		0		0		1		0
$A_2$	0		0		0		1		0		0
$E$	0		0		1		0		1		1
$T_1$	0		1		0		1		1		2
$T_2$	0		0		1		1		1		1

Table 4. Table 4: Fitted results using fit model Eq. ( III.2 ) with different t min subscript 𝑡 min t_{\rm min} at m π ∼ 938 ​ MeV similar-to subscript 𝑚 𝜋 938 MeV m_{\pi}\sim 938~{}\mathrm{MeV} .

$J^{PC}$	$t_{\min }$	$m_1{a}_t$	$m_2{a}_t$	$W_{11}$	$W_{12}$	$W_{21}$	$W_{22}$	${\chi }^2/d.o.f$
$A_1^{++}$	1	0.170(06)	0.556(09)	0.61(02)	0.36(02)	0.11(01)	0.85(01)	1.84
	2	0.168(07)	0.494(16)	0.59(03)	0.36(03)	0.06(02)	0.84(01)	0.44
	3	0.170(09)	0.495(26)	0.61(04)	0.34(05)	0.06(02)	0.84(02)	0.47
	4	0.169(12)	0.474(43)	0.59(06)	0.34(09)	0.05(03)	0.81(05)	0.51
	5	0.169(12)	0.546(67)	0.61(07)	0.43(14)	0.07(03)	1.02(19)	0.33
$A_1^{-+}$	1	0.307(09)	0.720(15)	0.66(02)	0.33(03)	0.15(2)	0.79(2)	2.51
	2	0.316(13)	0.665(28)	0.67(05)	0.28(05)	0.12(3)	0.79(2)	1.21
	3	0.306(19)	0.633(38)	0.62(08)	0.33(10)	0.09(04)	0.77(3)	1.39
	4	0.272(31)	0.530(57)	0.43(13)	0.53(14)	0.02(5)	0.67(4)	1.17
$E^{++}$	1	0.278(05)	0.691(09)	0.66(01)	0.32(01)	0.19(1)	0.77(11)	1.48
	2	0.278(07)	0.669(17)	0.66(02)	0.32(03)	0.18(2)	0.77(11)	1.92
	3	0.287(13)	0.568(33)	0.66(06)	0.26(07)	0.12(4)	0.72(2)	0.52
	4	0.280(20)	0.500(47)	0.60(10)	0.29(14)	0.05(6)	0.68(3)	0.43
	5	0.280(26)	0.499(78)	0.61(17)	0.29(20)	0.06(8)	0.67(6)	0.54
$T_2^{++}$	1	0.283(04)	0.657(06)	0.64(01)	0.33(01)	0.15(09)	0.80(1)	5.66
	2	0.285(06)	0.590(10)	0.63(02)	0.33(02)	0.10(1)	0.80(1)	1.63
	3	0.299(08)	0.564(18)	0.69(04)	0.22(05)	0.09(2)	0.78(1)	0.59
	4	0.293(13)	0.543(32)	0.64(07)	0.28(09)	0.06(3)	0.77(2)	0.63
	5	0.273(19)	0.517(55)	0.52(11)	0.45(13)	0.04(5)	0.74(7)	1.13

Table 5. Table 5: Fitted results using fit model Eq. ( III.2 ) with different t min subscript 𝑡 min t_{\rm min} at m π ∼ 650 ​ MeV similar-to subscript 𝑚 𝜋 650 MeV m_{\pi}\sim 650~{}\mathrm{MeV}

$J^{PC}$	$t_{\min }$	$m_1{a}_t$	$m_2{a}_t$	$W_{11}$	$W_{12}$	$W_{21}$	$W_{22}$	${\chi }^2/d.o.f$
$A_1^{++}$	1	0.176(09)	0.548(14)	0.62(03)	0.36(03)	0.09(02)	0.85(02)	1.71
	2	0.163(14)	0.481(23)	0.55(05)	0.42(06)	0.04(2)	0.85(2)	0.98
	3	0.173(17)	0.515(37)	0.59(07)	0.38(09)	0.06(3)	0.88(03)	0.99
	4	0.184(20)	0.549(87)	0.66(09)	0.27(13)	0.09(5)	0.92(15)	1.12
	5	0.150(25)	0.533(99)	0.49(11)	0.69(22)	0.06(4)	0.96(03)	0.84
$A_1^{-+}$	1	0.289(10)	0.710(16)	0.60(03)	0.38(03)	0.12(2)	0.82(2)	1.71
	2	0.320(15)	0.706(40)	0.70(06)	0.24(06)	0.15(4)	0.79(2)	0.71
	3	0.311(23)	0.673(60)	0.66(10)	0.29(13)	0.13(5)	0.77(4)	0.81
	4	0.286(41)	0.593(98)	0.52(19)	0.47(21)	0.07(8)	0.71(9)	0.93
$E^{++}$	1	0.268(07)	0.646(10)	0.62(02)	0.36(02)	0.14(2)	0.81(1)	1.42
	2	0.267(10)	0.601(20)	0.61(04)	0.35(04)	0.11(2)	0.80(1)	0.81
	3	0.255(13)	0.578(28)	0.55(05)	0.42(06)	0.09(3)	0.79(1)	0.71
	4	0.254(20)	0.563(49)	0.54(09)	0.42(11)	0.08(4)	0.77(6)	0.85
	5	0.272(36)	0.453(86)	0.62(24)	0.17(31)	0.00(11)	0.66(2)	0.46
$T_2^{++}$	1	0.283(06)	0.674(09)	0.66(02)	0.31(02)	0.18(1)	0.78(1)	2.67
	2	0.298(10)	0.627(21)	0.69(04)	0.25(04)	0.15(3)	0.77(2)	0.69
	3	0.289(15)	0.575(30)	0.66(06)	0.28(07)	0.11(4)	0.75(2)	0.40
	4	0.274(22)	0.542(46)	0.57(10)	0.39(13)	0.08(5)	0.74(4)	0.36
	5	0.266(30)	0.578(95)	0.54(15)	0.50(21)	0.09(6)	0.82(17)	0.33

Table 6. Table 6: Final results for the masses of the lowest state we obtain in the A 1 + + , A 1 − + , E + + superscript subscript 𝐴 1 absent superscript subscript 𝐴 1 absent superscript 𝐸 absent A_{1}^{++},A_{1}^{-+},E^{++} and T 2 + + superscript subscript 𝑇 2 absent T_{2}^{++} channels. These are the averaged values weighted by the inversed squared errors at each t min subscript 𝑡 min t_{\rm min} .

$m_{\pi }$(MeV)	$m_1\left(A_1^{++}\right)$(MeV)	$m_1\left(E^{++}\right)$(MeV)	$m_1\left(T_2^{++}\right)$(MeV)	$m_1\left(A_1^{-+}\right)$(MeV)
$938$	1417(30)	2332(31)	2392(26)	2573(55)
$650$	1498(58)	2294(43)	2474(39)	2585(65)

Table 7. Table 7: We compare our results with previous results both from the quenched lattice QCD studies  Morningstar and Peardon ( 1999 ); Chen et al. ( 2006 ) and the full-QCD study  Gregory et al. ( 2012 ) . We average the masses of E + + superscript 𝐸 absent E^{++} and T 2 + + superscript subscript 𝑇 2 absent T_{2}^{++} states to obtain the estimate of the 2 + + superscript 2 absent 2^{++} glueball mass.

	$m_{\pi }$ (MeV)	$m_{0^{++}}$ (MeV)	$m_{2^{++}}$ (MeV)	$m_{0^{-+}}$ (MeV)
$N_f=2$	$938$	1417(30)	2363(39)	2573(55)
	$650$	1498(58)	2384(67)	2585(65)
$N_f=2+1$ Gregory et al. (2012)	$360$	1795(60)	2620(50)	—
quenched Morningstar and Peardon (1999)	—	1710(50)(80)	2390(30)(120)	2560(35)(120)
quenched Chen et al. (2006)	—	1730(50)(80)	2400(25)(120)	2590(40)(130)

Table 8. Table 8: The fitting details for η ′ superscript 𝜂 ′ \eta^{\prime} meson mass from topological charge density correlator at m π ∼ 938 ​ MeV similar-to subscript 𝑚 𝜋 938 MeV m_{\pi}\sim 938~{}\mathrm{MeV} and m π ∼ 650 ​ MeV similar-to subscript 𝑚 𝜋 650 MeV m_{\pi}\sim 650~{}\mathrm{MeV} respectively.

$m_{\pi }$	fit range($a_s$)	$m_{{\eta }^{\prime }}{a}_s$	$m_{{\eta }^{\prime }}\left(\mathrm{MeV}\right)$	${\chi }^2/dof$
$938\mathrm{MeV}$	3.74-5.92	$0.856\left(21\right)$	$1481\left(36\right)$	$1.01$
$650\mathrm{MeV}$	3.87-5.48	$0.514\left(22\right)$	$890\left(38\right)$	$1.43$

Table 9. Table 9: The table collects the masses of flavor singlet pseudoscalar mesons from the quenched and unquenched lattice QCD studies. P ​ ( x ) 𝑃 𝑥 P(x) , q ​ ( x ) 𝑞 𝑥 q(x) and Φ ( PS ) superscript Φ PS \Phi^{({\rm PS})} stand for the quark bilinear pseuscalar operator, the topological charge density, and the pseudoscalar glueball operator, respectively.

	$P\left(x\right)$	$q\left(x\right)$	${\Phi }^{\left(\mathrm{PS}\right)}$
$N_f=0$	——	2563(34) MeV Chowdhury et al. (2015)	2590(140) MeV Chen et al. (2006)
$N_f=2$	768(24) MeV Helmes et al. (2017)	890(38) MeV (this work)	2585(65) MeV (this work)
$N_f=2+1$	947(142) MeV Christ et al. (2010)	1019(119) MeV Fukaya et al. (2015)	——
$N_f=2+1+1$	1006(65)MeV  Michael et al. (2013)	——	——

Equations

$$S_g=$$

$$R_{ss}(x,y)$$

$$R_{st}(x,t)$$

$$L({\xi }_g)=\sum _{x,y}\frac{(R_{ss}(x,y)-R_{st}(x,{\xi }_{g}y){)}^2}{(\Delta R_s{)}^2+(\Delta R_t{)}^2}$$

$$S_f$$

$${\hat{W}}_{\mu }$$

$${\nabla }_{\mu }f(x)$$

$${\Delta }_{\mu }f(x)$$

$$E(p{)}^2{a}_t^2=m^2{a}_t^2+\frac{|\vec{p}{|}^2{a}_s^2}{{\xi }_f^2}.$$

$$\frac{r_0}{a_s}=\sqrt{\frac{1.65+\alpha }{\sigma a_s^2}}$$

$$C_{\alpha \beta }^{(R)}(t)=\sum _{\tau }⟨0|{\bar{\varphi }}_{\alpha }^{(R)}(t+\tau ){\bar{\varphi }}_{\beta }^{(R)†}(\tau )|0⟩,\alpha ,\beta =1,...,24$$

$${\bar{\varphi }}_{\alpha }^{(R)}(t)={\varphi }_{\alpha }^{(R)}(t)-⟨0|{\varphi }_{\alpha }^{(R)}(t)|0⟩.$$

$${\mathbf{C}}^{(R)}(t_0){\mathbf{v}}_i^{(R)}={\lambda }_i(t_0){\mathbf{C}}^{(R)}(0){\mathbf{v}}_i^{(R)},$$

$${\Phi }_i^{(R)}(t)=\sum _{\alpha =1}^{24}{v}_{i,\alpha }^{(R)}{\bar{\varphi }}_{\alpha }^{(R)}(t).$$

$${\tilde{C}}_i^{(R)}(t)=\sum _{\tau }⟨0|{\Phi }_i^{(R)}(t+\tau ){\Phi }_i^{(R)†}(\tau )|0⟩,$$

$$m_{i,\mathrm{eff}}^{(R)}(t)=\ln \left(\frac{{\tilde{C}}_i^{(R)}(t)}{{\tilde{C}}_i^{(R)}(t+1)}\right).$$

$${\bar{C}}_i^{A_1^{++}}(t)={\tilde{C}}_i^{A_1^{++}}(t)-{\tilde{C}}_i^{A_1^{++}}(t+\delta t),$$

$${\bar{C}}_i^{A_1^{++}}(t)=\sum _j{W}_{ij}^{A_1^{++}}(1-e^{-m_j\delta t})e^{-m_{j}t}\equiv \sum _j{W}_{ij}^{\prime A_1^{++}}{e}^{-m_{j}t},$$

$${\tilde{C}}_1^{(R)}(t)$$

$${\tilde{C}}_2^{(R)}(t)$$

$$m_{i,\mathrm{eff}}^{(R)}(t)$$

$${\sigma }_C^2\approx \sum _{i,j=1}^6\frac{\partial C}{\partial a_i}{\sigma }_{ij}\frac{\partial C}{\partial a_j},$$

$$C(t)=⟨0|\mathcal{O}(t){\mathcal{O}}^{†}(0)|0⟩=\sum _n⟨0|\mathcal{O}|n⟩⟨n|{\mathcal{O}}^{†}|0⟩e^{-m_{n}t},$$

$$\sum _n|n⟩⟨n|=1,⟨n|n^{\prime }⟩={\delta }_{n{n}^{\prime }}.$$

$${\partial }_{\mu }{J}_5^{\mu }(x)=2mP(x)-\frac{N_f{g}^2}{16{\pi }^2}{ϵ}_{\mu \nu \rho \sigma }Tr{F}^{\mu \nu }{F}^{\rho \sigma },$$

$$C_q(x-y)=⟨q(x)q(y)⟩,$$

$${\chi }_t=\frac{1}{V_4}\int d^{4}x{d}^{4}y{C}_q(x-y),$$

$$C_q(x-y)=A{\delta }^4(x-y)+{\bar{C}}_q(x-y),$$

$${\bar{C}}_q(r)=N\frac{m_{\mathrm{PS}}}{4{\pi }^{2}r}{K}_1(m_{\mathrm{PS}}r),$$

$$K_1(z)\sim \sqrt{\frac{\pi }{2z}}{e}^{-z}(1+\frac{3}{8z}),|argz|<\frac{3}{2}\pi .$$

$$\frac{1}{q^2-m_{{\eta }^{\prime }}^2}$$

$$m_{{\eta }^{\prime }}^2\approx m_{\pi }^2+m_0^2.$$

Full text

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 $q\bar{q}$ mesons described in the constituent quark model. Glueballs have been intensively studied by lattice QCD and other theoretical methods Jaffe and Johnson (1976); Cornwall and Soni (1983); Hou and Soni (1984); Brower et al. (2000); Szczepaniak and Swanson (2003); Narison (2006); Sanchis-Alepuz et al. (2015), for more details of this subject, see reviews in Klempt and Zaitsev (2007); Mathieu et al. (2009); Crede and Meyer (2009); Ochs (2013). 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 Morningstar and Peardon (1997, 1999); Chen et al. (2006).

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/\psi$ 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 Gui et al. (2013); Yang et al. (2013). 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/\psi\rightarrow\gamma\eta\eta$ Ablikim et al. (2013) and $J/\psi\rightarrow\gamma\phi\phi$ Ablikim et al. (2016).

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 Bali et al. (2000); Hart and Teper (2002); Richards et al. (2010); Gregory et al. (2012). However, for the mass of the pseudoscalar glueball, a consensus has not been reached. For example, in Ref. Richards et al. (2010) 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. Gregory et al. (2012). 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 ($\eta^{\prime}$ in the $N_{f}=2+1$ case) in the presence of dynamical quarks. Therefore, it is mandatory to identify the contribution of the $\eta^{\prime}$ 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_{\pi}\sim 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 $\eta^{\prime}$ could be probed by the topological charge density operator Fukaya et al. (2015). In contrast, a similar study in the quenched approximation found a pseudoscalar with a mass compatible with that in the pure gauge theory Chowdhury et al. (2015). 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 Morningstar and Peardon (1997):

[TABLE]

where $P_{ij}$ is the usual plaquette variable and $R_{ij}$ is the $2\times 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 $\gamma_{g}$ designates the gauge aspect bare ratio of the anisotropic lattice, denoted as $\xi_{0}$ in our former quenched studies Su et al. (2006). Although $\gamma_{g}$ suffers only small renormalization with the tadpole improvement Lepage and Mackenzie (1993), we have to tune it by determining the renormalized anisotropy ratio $\xi_{g}$. As for the tadpole improvement parameter $u_{t}$ for temporal gauge links, we take the approximation $u_{t}\approx 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 Umeda et al. (2003); Klassen (1998). We measure the ratios

[TABLE]

and expect the spatial and temporal behaviors being the same at the correct $\xi_{g}$.

Therefore we find $\xi_{g}$ by minimizing

[TABLE]

where $\Delta R_{s}$ and $\Delta R_{t}$ are the statistical errors of $R_{ss}$ and $R_{st}$. We interpolate $R_{st}(x,\xi_{g}y)$ and its error with a cubic spline interpolation at non-integer $\xi_{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\in\{2,3,4,5\}$.

We adopt the anisotropic clover fermion action in the fermion sector Edwards et al. (2008):

[TABLE]

where $\hat{F}_{\mu\nu}=\frac{1}{4}{\rm Im}(P_{\mu\nu}(x))$ and the dimensionless Wilson operator reads

[TABLE]

The bare fermion aspect ratio $\gamma_{f}$ is also tuned to make sure that the measured aspect ratio $\xi_{f}\approx\xi_{g}\approx\xi=5$. $\xi_{f}$ is measured from the dispersion relation of the pseudoscalar and vector mesons

[TABLE]

where $\vec{p}=2\pi\vec{k}/L_{s}$ is the momentum on the lattice with periodic spatial boundary conditions.

We generate two gauge ensembles on the $12^{3}\times 128$ anisotropic lattice at $\beta=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}+\alpha/r+\sigma r$. Using the Sommer scale parameter $r_{0}^{-1}=410(20)$ MeV defined through $r^{2}\frac{dV(r)}{dr}|_{r=r_{0}}=1.65$, we can determine the ratio

[TABLE]

where $\alpha$ and $\sigma a_{s}^{2}$ are derived from the fit to calculated potential $V(r)=V(\hat{r}a_{s})$ with $\hat{r}$ 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_{0}^{-1}=410(20)$ MeV. The ensemble parameters are listed in Table 1, where we also give the physical values of $a_{t}^{-1}$ 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_{\pi}$’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 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 $\eta^{\prime}\pi$ (there is only one $I=0$ pseudoscalar meson for $N_{f}=2$, which is taken as the counterpart of the (approximately) flavor-singlet $\eta^{\prime}$ in the $N_{f}=3$ case). At $m_{\pi}\sim 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 $\eta^{\prime}\pi$ threshold. At $m_{\pi}\sim 650$ MeV, $m_{\eta^{\prime}}$ is estimated to be $m_{\eta^{\prime}}\sim 890$ MeV (see below in Sec. 4), thus the mass value of $1362(53)$ MeV is also below the $\eta^{\prime}\pi$ 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 well-described 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 Chen et al. (2006) for further details.

III.1 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^{PC}$ quantum number should be denoted by $R^{PC}$ with $R$ the irreducible representation of $O$ which may include the components of $J$ in the continuum limit. The parity $P=\pm$ and the charge conjugation $C=\pm$ 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\geq 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}$, tensor states with $J=2$ are reduced to direct sum of $E$ and $T_{2}$, i.e. $(J=2)\downarrow O=E\bigoplus T_{2}$.

As described in Chen et al. (2006), 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 Chen et al. (2006). 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, $\{\phi^{(R)}_{\alpha},\alpha=1,2,\ldots,24\}$, for each $R^{PC}$.

Based on these operator sets, we use the variational method to get the optimized operators $\mathcal{O}^{(R)}$ which mostly project to specific glueball states. In each symmetry channel $R$, we first calculate the $24\times 24$ correlation matrix $C^{(R)}(t)$,

[TABLE]

where $\bar{\phi}^{(R)}_{\alpha}$ is the vacuum-subtracted operator of $\phi^{(R)}_{\alpha}$,

[TABLE]

In practice, we only apply the vacuum subtraction to the operators in $A_{1}^{++}$ channel. Secondly, we solve the following generalized eigenvalue problem,

[TABLE]

where $\mathbf{v}_{i}^{(R)}$ is the $i$-th eigenvector, and $\lambda_{i}\equiv e^{-\bar{m}_{i}(t_{0})t_{0}}$ is the $i$-th eigenvalue where $\bar{m}_{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 $\mathbf{v}^{(R)}_{i}$ gives the linearly combinational coefficients of operators $\bar{\phi}^{(R)}_{\alpha}$ to build an optimal operator $\Phi^{(R)}_{i}$ which overlaps mostly to the $i$-th state,

[TABLE]

III.2 Data analysis

In this work, the correlation function of the optimal operator $\Phi^{(R)}_{i}$ for the $i$-th state is calculated as

[TABLE]

where we do the summation over the temporal direction to increase the statistics. Accordingly, the effective mass is defined as

[TABLE]

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,\,{\rm 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+\delta t)$,

[TABLE]

whose spectral expression is

[TABLE]

where $W_{ij}^{A_{1}^{++}}$ is the spectral weight of the $j$-th state in $\tilde{C}_{i}^{A_{1}^{++}}(t)$. Obviously, the possible constant term cancels with the spectrum unchanged. In practice, we take $\delta t=5a_{t}$.

We focus on the $R^{PC}=A_{1}^{++},A_{1}^{-+},E^{++}$, and $T_{2}^{++}$ channels in this work. For all these channels, the effective masses of $\tilde{C}_{1}^{R}(t)$ (red points) and $\tilde{C}_{2}^{R}(t)$ (blue points) ($\bar{C}_{i}^{A_{1}^{++}}$ for $A_{1}^{++}$ channel) are plotted in Fig. 2, 3, 4 and 5, respectively. In each figure, the left panel shows the result at $m_{\pi}\sim 938~{}\mathrm{MeV}$, and the right panel is for $m_{\pi}\sim 650~{}\mathrm{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 $q\bar{q}$ 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 $\Phi^{(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_{1,{\rm eff}}^{(R)}(t)$ tends to reach a plateau as $t$ increases, $m_{2,{\rm eff}}^{(R)}(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}\approx 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 to $\tilde{C}_{1}^{(R)}(t)$ and $\tilde{C}_{2}^{(R)}(t)$ simultaneously through the following function forms,

[TABLE]

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_{\rm max}$’s of the fit windows of $\tilde{C}_{1}^{(R)}(t)$ and $\tilde{C}_{2}^{(R)}(t)$ are chosen properly to include only the data points with good signal-to-noise ratios (The $t_{\rm max}$ of $\tilde{C}_{2}^{(R)}(t)$ are set to be from $7a_{t}$ to $9a_{t}$, while $t_{\rm max}$ of $\tilde{C}_{1}^{(R)}(t)$ can be larger than $10a_{t}$). Actually, the fit results are insensitive to $t_{\rm 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_{\rm max}$’s fixed and vary $t_{\rm 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 4 and  5. Except for $t_{\rm min}=1$ case in $T_{2}^{++}$ channel, all other fits are acceptable with reasonable $\chi^{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_{\rm min}$ increasing gradually. This signals that our fitting model in Eq. (III.2) 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_{11}^{(R)}$, $W_{12}^{(R)}$,$W_{21}^{(R)}$, and $W_{22}^{(R)}$ at different $t_{\rm min}$, which are listed in Table 4 and 5. The red and blue bands are obtained through the function

[TABLE]

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,

[TABLE]

where $C$ denotes $m_{i,{\rm eff}}^{(R)}(t)$, $a_{i}$’s are the six parameters in Eq. (III.2) and $\sigma_{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 $\chi^{2}$. It is expected that the fitted $m_{2}$ is generally (much) higher than the mass of the second state.

As shown in Table 4 and  5, most of the fits using different $t_{\rm{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_{\rm{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_{\rm{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_{\pi}\sim 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_{\pi}\sim 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_{\pi}\sim 650$ MeV, which is roughly one-half as large as that at $m_{\pi}\sim 938$ MeV (see Table  1). From Table 2 and Table 6 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 $q\bar{q}$ and the scalar glueball, which needs further investigation.

III.3 Interpretation of the ground states

Generally speaking, the two-point function of an interpolating operator $\mathcal{O}(t)$ with definite quantum numbers is usually parameterized as

[TABLE]

where $\{|n\rangle,n=1,2,\ldots\}$ are eigenstates of Hamiltonian with eigenvalue $m_{n}$, which make up an orthogonal, normalized, and complete state set with

[TABLE]

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^{PC}=0^{++}$ can be a bare scalar glueball $|G_{0^{++}}\rangle$, a bare $q\bar{q}$ scalar meson $|f_{0}\rangle$, or even $\pi\pi$ scattering states $|\pi\pi\rangle$. 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\rangle$ and its adjacent state, which could be of nature $|\pi\pi\rangle$ or $|f_{0}\rangle$. This then yields the fitting model in Eq. (III.2) 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 mass of the bare glueball states $|G\rangle$ can still be obtained by assuming the corresponding operators $\mathcal{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}^{-+}}\sim 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 $q\bar{q}$ pseudoscalar meson is expected to exist, but we do not observe this state from the correlation function of the glueball operator $\Phi^{(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),

[TABLE]

where $g$ is the strong coupling constant, $P(x)=\bar{\psi}(x)\gamma_{5}\psi(x)$ is the pseudoscalar density, and the anomalous gluonic operator $\epsilon_{\mu\nu\rho\sigma}F^{\mu\nu}F^{\rho\sigma}$ 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 $\eta^{\prime}$).

The correlation function of $q(x)$ is expressed as

[TABLE]

from which one can get the topological susceptibility

[TABLE]

where $V_{4}$ is the four-dimensional volume of the Euclidean spacetime. It is known that $\chi_{t}$ is positive and takes a value $\sim(180~{}MeV)^{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

[TABLE]

where $\bar{C}_{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 from $\bar{C}_{q}(x-y)$. It is expected that $\bar{C}_{q}(x-y)$ would be dominated by the contribution of the lowest pseudoscalar meson in the large $r$ range and can be parameterized as Shuryak and Verbaarschot (1995)

[TABLE]

where $N$ is an irrelevant normalization factor, $m_{\rm 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

[TABLE]

Therefore, one can obtain $m_{\rm 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 Chowdhury et al. (2015) and full QCD calculations Fukaya et al. (2015). In the quenched approximation, the extracted $m_{\rm PS}=2563(34)$ MeV is in good agreement with the pseudoscalar glueball mass $m_{\rm 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_{\pi}$, $m_{\rm PS}$ is obtained to be $1013(117)$ MeV, which is consistent with the mass of the physical $\eta^{\prime}$.

In this work, we adopt a similar strategy to that in Fukaya et al. (2015). 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 Moran and Leinweber (2008); Fukaya et al. (2015). The Wilson flow provides a reference energy scale $\frac{1}{\sqrt{8t}}$ Lüscher (2010). In practice, we use the code published by the BMW collaboration Borsanyi et al. (2012) to evaluate the topological charge density. Fig. 7 shows $C_{q}(r)$ for $m_{\pi}\sim 938~{}\mathrm{MeV}$ and $m_{\pi}\sim 650~{}\mathrm{MeV}$ at flow times $t=0.2,0.3,0.4,0.8$ respectively. On our lattices, these $t$ values correspond to $\sqrt{8t}\sim 0.15,0.18,0.21$ and $0.30~{}\mathrm{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_{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\sim 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}\in[3.8,5.4]$.

We carry out a correlated minimal-$\chi^{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 $\chi^{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_{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 $\Phi^{PS}$. Thus the light pseudoscalar state observed in $C_{q}(r)$ can be naturally assigned to be the flavor singlet $q\bar{q}$ state $\eta^{\prime}$. Theoretically, the mass of $\eta^{\prime}$ is acquired through the interaction of sea quark loops according to the Witten-Veneziano mechanism Witten (1979); Veneziano (1979). In this mechanism, the propagator of $\eta^{\prime}$ can be expressed as

[TABLE]

where the parameter $m_{0}^{2}$ is introduced to describe the gluonic coupling, such that

[TABLE]

On the other hand, $m_{0}^{2}$ is related to the topological susceptibility $\chi_{t}$ through

[TABLE]

where $f_{\pi}$ is the decay constant of $\pi$. For our case of $N_{f}=2$, if we take the values $\chi_{t}=(180\,{\rm MeV})^{4}$, $f_{\pi}\sim 150$ MeV for $m_{\pi}\sim 650$ MeV and $f_{\pi}\sim 200$ MeV for $m_{\pi}\sim 938$ MeV, $m_{0}^{2}$ is estimated to be approximately $(610\,{\rm MeV})^{2}$ and $(460\,{\rm MeV})^{2}$, respectively. Thus the $\eta^{\prime}$ mass can be derived as $m_{\eta^{\prime}}\sim 890$ MeV for $m_{\pi}\sim 650$ MeV, and $m_{\eta^{\prime}}\sim 1045$ MeV for $m_{\pi}\sim 938$ MeV. These values are not far from the $m_{\rm 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 $q\bar{q}$ pseudoscalar meson corresponding to the $\eta^{\prime}$ meson in the real world, which can be accessed by the topological charge density operator.

Now that the $\eta^{\prime}$ 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 $\Phi^{\rm PS}$). In order to clarify this, we check the continuum form of $\Phi^{({\rm 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,

[TABLE]

where $R$ stands for each rotation operation in $O$ group, $c_{R}^{A_{1}}$ is the combinational coefficients corresponding to the $A_{1}$ irreducible representation, $W_{\alpha}$ is any of the four prototypes made up of a specifically smeared gauge links. According to the non-Abelian Stokes theorem Simonov (1989), a rectangle Wilson loop $P_{\mu\nu}^{a\times b}(x)$ of size $a\times b$, with $a,b$ small, can be expanded as

[TABLE]

where $F_{\mu\nu}$ is the strength tensor of the gauge field. For simplicity, the factor $ig$ is absorbed into the quantity $F_{\mu\nu}$. The small $ab$ expansion of $P_{\pm\mu\pm\nu}(x)$ is similar to Eq. 30 by replacing $a$ and $b$ with $\pm a$ and $\pm 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,

[TABLE]

which is obviously different from the anomalous part of the PCAC relation, $\epsilon_{\mu\nu\rho\sigma}F^{\mu\nu}(x)F^{\rho\sigma}(x)\propto\mathbf{E}(x)\cdot\mathbf{B}(x)$. Actually, the operator $\Phi^{({\rm 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 $q\bar{q}$ 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. Chowdhury et al. (2015) use $q(x)$ as pseudoscalar operators and derive the ground state mass $m_{\rm 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 Morningstar and Peardon (1999) and $2.590(140)$ GeV Chen et al. (2006). 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 $\eta^{\prime}$ mass in the lattice simulation with dynamical quarks, and have given the results $m_{\eta^{\prime}}=768(24)$ MeV ($N_{f}=2$) Helmes et al. (2017), $m_{\eta^{\prime}}=947(142)$ MeV ($N_{f}=2+1$) Christ et al. (2010) and $m_{\eta^{\prime}}=1006(65)$ MeV ($N_{f}=2+1+1$) Michael et al. (2013), which almost reproduce the experimental result $m_{\eta^{\prime}}=958$ MeV. When the $q(x)$ operator is applied, $N_{f}=2+1$ lattice simulation gives the result $m_{\eta^{\prime}}=1019(119)$ MeV at the physical pion mass Fukaya et al. (2015), 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_{\rm PS}=890(38)$ MeV at $m_{\pi}=650$ MeV, which is compatible with the $m_{\eta^{\prime}}=768(24)$ MeV above (note that our $m_{\pi}$ is higher than that in Ref. Helmes et al. (2017)). The similar result for $m_{\eta^{\prime}}$ 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 $\eta^{\prime}$ meson. In contrast, the glueball operator $\Phi^{({\rm 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_{\pi}\sim 650$ MeV and $m_{\pi}\sim 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_{\pi}\sim 938$ MeV and $650$ MeV, respectively. In the pseudoscalar channel, using the gluonic operator whose continuum limit has the form of $\epsilon_{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 state derived in this way are much lighter (around 1 GeV) and compatible with the expected masses of the flavor singlet $q\bar{q}$ meson. This provides a strong hint that the operator $\epsilon_{ijk}TrB_{i}D_{j}B_{k}$ and the topological charge density (proportional to $Tr\mathbf{E}\cdot\mathbf{B}$) couple rather differently to the glueball states and $q\bar{q}$ 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.

Acknowledgements

The numerical calculations are carried out on Tianhe-1A at the National Supercomputer Center (NSCC) in Tianjin and the GPU cluster at Hunan Normal University. This work is supported in part by the National Science Foundation of China (NSFC) under Grants No. 11575196, No. 11575197, No. 11335001, No. 11405053, No. 11405178 and No. 11275169. Y. C., Z. L. and C. L. also acknowledge the support of NSFC under No. 11261130311 (CRC 110 by DFG and NSFC). Y. C. thanks the support by the CAS Center for Excellence in Particle Physics (CCEPP). C.L. is also funded in part by National Basic Research Program of China (973 Program) under code number 2015CB856700. M. G. thanks the support by the Youth Innovation Promotion Association of CAS (2015013).