$\eta$-glueball mixing from $N_f=2$ lattice QCD

We perform the first lattice study on the mixing of the isoscalar pseudoscalar meson $\eta$ and the pseudoscalar glueball $G$ in the $N_f=2$ QCD at the pion mass $m_\pi\approx 350$ MeV. The $\eta$ mass is determined to be $m_\eta=714(6)(16)$ MeV. Through the Witten-Veneziano relation, this value can be matched to a mass value of $\sim 981$ MeV for the $\mathrm{SU(3)}$ counterpart of $\eta$. Based on a large gauge ensemble, the $\eta-G$ mixing energy and the mixing angle are determined to be $|x|=107(15)(2)$ MeV and $|\theta|=3.46(46)^\circ$ from the $\eta-G$ correlators that are calculated using the distillation method. We conclude that the $\eta-G$ mixing is tiny and the topology induced interaction contributes most of $\eta$ mass owing to the QCD $\mathrm{U_A(1)}$ anomaly.


I. INTRODUCTION
Chiral symmetry is an intrinsic symmetry of quantum chromodynamics (QCD) in the massless limit of quarks and is spontaneously broken due to the nonzero quark condensate. The spontaneous symmetry breaking (SSB) of the SU L (3) × SU R (3) into the flavor SU(3) generates eight Goldstone particles sorted into the pseudoscalar octet composed of {π, K, η 8 } which are slightly massive due to the small u, d, s quark masses. If SSB also applies to the U L (1) × U R (1) sector of the chiral symmetry, it breaks into U V (1) that corresponds to the baryon number conservation expects the existence of an additional Goldstone particle, namely, a light flavor singlet pseudoscalar meson. The η ′ meson is predominantly a flavor singlet but is too massive to be taken as a candidate for this Goldstone particle. The η ′ mass puzzle has a direct connection with the QCD U A (1) anomaly that the anomalous gluonic term, the topological charge density, breaks the conservation of the flavor singlet axial current even in the chiral limit. Even though the anomalous axial vector relation can be written as the divergence of a gauge variant axial vector, which is zero and implies a U A (1) symmetry, Kogut and Susskind [1] pointed out that its spontaneous breaking generates a massless mode that violates the Gell-Mann-Okubo relation and then renders η ′ more massive. With respect to the nontrivial topology of QCD vacuum, Witten [2] and Veneziano [3] proposed a mechanism for the original of η ′ mass that the nonperturbative coupling of the topological charge density and the flavor singlet pseudoscalar induces a self energy correction m 2 0 , which is proportional to the topological susceptibility χ of gauge fields. Using the physical mass of η ′ , the value of χ is estimated to be around (180 MeV) 4 , which is supported by lattice QCD calculations [4]. Another interesting property of η ′ is its large production rate in the J/ψ radiative decays [5], which are abundant in gluons and are expected to favor the production of glueballs. This observation along with the original mechanism of η ′ mass manifests the strong coupling of η ′ to gluons and thereby prompts the conjecture that η ′ may mix substantially with pseudoscalar glueballs since they have the same quantum number. The KLOE Collaboration analyzed the processes ϕ → γη and ϕ → γη ′ and found that the η ′glueball mixing might be required and the mixing angle can be as large as (22 ± 3) • [6]. In contrast, another phenomenological analysis of the KLOE result gave the mixing angle (12 ± 13) • which is consistent with zero within the large error [7]. A phenomenological analysis of the processes J/ψ(ψ ′ ) decaying into a pseudoscalar and a vector final states obtained the η ′ -glueball mixing angle to be around 9 • by considering the η-η ′ -glueball mixing model [8]. Obviously, the determined mixing angle varies in a fairly large range. Based on the η ′ -glueball mixing picture, there have been theoretical discussions on the possibility of η(1405) as a pseudoscalar glueball candidate [8][9][10][11][12]. However, the quenched lattice QCD studies [13][14][15] predict that the mass of the pseudoscalar glueball is around 2.4-2.6 GeV, which is confirmed by lattice simulations with dynamical quarks [16][17][18][19]. This raised a question on η(1405) as a glueball candidate because of its much lighter mass. On the other hand, the strong hint for η(1405) to be a glueball candidate is the observation that there exist three isoscalar pseudoscalar mesons η(1295), η(1405) and η(1475) such that one of them is surplus according to the quark model. If η(1405) and η(1475) are the same state [20], there is no need for a pseudoscalar glueball state in this mass region. Some mixing model studies also favor the pseudoscalar glueball to have a mass heavier than 2 GeV [21,22].
In this work, we will investigate the possible mixing of isoscalar pseudoscalar meson and the pseudoscalar glueball in N f = 2 lattice QCD. The isoscalar pseudoscalar meson is named η throughout this work, which is the SU(2) counterpart of the SU(3) flavor singlet (approximately η ′ ) in the N f = 3 case. We have arXiv:2205.12541v2 [hep-lat] 25 May 2023 generated a large ensemble of gauge configurations with N f = 2 degenerated u, d quarks at a pion mass m π ≈ 350 MeV, so we can make a precise determination of η mass. The calculation of η ′ mass (and the η ′ − η mixing) has been performed in several N f = 2 + 1 lattice QCD studies, whose results are in agreement with the physical value after the chiral extrapolation [23][24][25][26]. A systematic and comprehensive lattice study on the properties of η and η ′ is presented in Ref. [27], where the masses, the decay constants and the gluonic matrix elements of η and η ′ have been derived to a high precision. There are also many studies on the η mass from N f = 2 lattice QCD [19,[28][29][30]. According to the Witten-Veneziano mechanism (WV), in the N f = 2 case, pion mass and η mass are related as m 2 η = m 2 π + m 2 0 , where m 2 0 is the correction from the topology induced interaction and is proportional to N f . As a check of WV, we would like to take a look at this relationship and use the obtained m 2 0 to predict the η ′ mass in the physical case (a pioneering work following this way in the quenched approximation can be found in Ref. [31]). After that, we calculate the correlation functions of the η operator and the glueball operator, from which the mixing angle can be extracted. The strategy of the study is similar to that used in the η c -glueball mixing [32]. As an exploratory study, we tentatively treat the pseudoscalar glueball as a stable particle and ignore its resonance nature in the presence of light sea quarks. Obviously, the numeric task involves the calculation of the annihilation diagrams of u, d quarks, so we adopt the distillation method [33] which enables the gauge covariant smearing of quark fields and the allto-all quark propagators (perambulators) simultaneously. Since we also have the perambulators of the valence charm quark, we also calculate the η-η c correlation functions and explore their properties. This paper is organized as follows: In Sec. II we describe the lattice setup, operator construction and formulation of correlation functions. Section III gives the theoretical formalism of the meson-glueball mixing, where the data analyses and the results can be found. The discussion and summary are given in Sec. IV, and the preliminary results from η-η c correlation functions are presented here.

A. Lattice setup
We generate gauge configurations with N f = 2 degenerate u, d on an L 3 × T = 16 3 × 128 anisotropic lattice. We use the tadpole-improved Symanzik's gauge action for anisotropic lattices [13,15] and the tadpoleimproved anisotropic clover fermion action [19,34]. The parameters in the action are tuned to give the anisotropy ξ = a s /a t ≈ 5.3, where a t and a s are the temporal and spatial lattice spacings, respectively. The anisotropy ξ ≈ 5.3 is checked by the dispersion relation of the pseudoscalar π (along with that of η calculated using the distillation method, see Sec. III), which takes the continuum form where X refers to a specific hadron state and ⃗ p is the spatial momentum ⃗ p = 2π Las ⃗ n on the lattice, and ⃗ n is the mode of spatial momentum. Figure 1 shows the energies obtained from the correlation functions of π, η and ρ at different spatial momentum modes up to ⃗ n = (1, 2, 2). The data points of π and η fall on straight lines perfectly and can be well described by Eq. (1) with ξ = 5.365(5) and 5.34(3), respectively (illustrated as shaded lines in the figure). The fit to energies in the ρ channel using Eq. (1) gives ξ = 5.58(1) which deviates from 5.3 drastically. This can be definitely attributed to the influence of nearby P -wave ππ scattering states that have the same center-of-mass momentum as ρ. So we do not use ρ meson to check the ξ value.
There are subtleties in the determination of the lattice spacing for our lattice setup. We intend to generate gauge configurations at a pion mass m π ∼ 300 MeV. As the first step, we calculate the static potential through Wilson loops to determine e c and the string tension σa 2 s (in lattice units). Then the spatial lattice spacing a s is determined in the unit of the Sommer's scale parameter r 0 according to the condition r 2 dV dr r=r0 =  [5]. Here n refers to the light u, d quarks. The right most column lists the m 2 V − m 2 P S (GeV 2 ). In the row of ss states, the mass of the ss pseudoscalar ηs is determined by the HPQCD collaboration from lattice QCD calculations [35].
We tentatively use the value r 0 = 0.491(9) fm, which is determined at the physical point (in the chiral, continuum and infinite volume limits) [36], to set a s of our lattice, from which the quark mass is tuned to give m π ≈ 300 MeV. However, given this value of a s , the mass of the vector meson ρ is determined to be around 750 MeV, which is obviously lower than expected (note that m ρ is 770 MeV at the physical pion mass m π ∼ 135 MeV). This can be understood as follows. The measured value of r 0 /a s actually also depends on the quark mass for a given lattice, such that an additional quark mass dependence is introduced to the values of m π and m ρ in terms of r 0 . Since our quark mass is substantially far from the chiral limit, it is more reasonable to use the value of r 0 that is extrapolated from the physical pion mass to the comparable one in our situation. An alternative scale setting scheme is to choose another quantity that is insensitive to quark masses. Experimentally, there is an interesting relation between pseudoscalar meson masses m P S and the vector meson masses m V of the quark configuration q lq , where q l stands for the u, d, s quark and q stands for u, d, s, c quarks. The PDG results of the masses of these vector and pseudoscalar mesons [5] are collected in Table I along with their mass squared differences. Even though the reason is still unknown for this relation, empirically these values are insensitive to quark masses. On the other hand, the mass of η s , the ss counterpart of π (not considering the ss annihilation effects in the calculation), is determined to be m ηs = 0.686(4) GeV from lattice QCD by the HPQCD collaboration [35]. Even though η s is not a physical state, the mass squared difference m 2 ϕ − m 2 ηs ≈ 0.570 GeV 2 also satisfies the empirical relation in Eq. (4). In this study, the dimensionless masses of π and ρ is determined to be m π a t = 0.05055 (13) and m ρ a t = 0.12046 (20). In this sense, we assume the relation of Eq. (4) is somewhat general for heavy-light mesons and use it to set the scale parameter a t . Of course, one should take caution to use this relation since ρ is experimentally a wide resonance and decays into P -wave ππ states by 99%. On our lattice the lowest P -wave ππ energy threshold in the rest frame of ππ is 2E π (p)a t ≈ 0.1795 with ξ ≈ 5.3, which is substantially higher than m ρ a t . This means ρ in its rest frame is a stable particle, such that the mass value m ρ is reliable. In practice, we make the least squares fitting to the mass squared differences over the nn, ns, nc and sc systems where n refers to the u, d quarks, and get the value ∆m 2 = 0.568(8) GeV 2 , which serves as an input to give the lattice scale parameter a −1 t = 6.894(51) GeV and the corresponding spatial lattice spacing a s ≈ 0.1517(11) fm. We emphasize that the errors quoted in the values of a t and a s include only the statistical errors of m π a t and m ρ a t , as well as the uncertainty of ∆m 2 . Accordingly, the u, d mass parameter in this study gives m π = 348.5(1.0)(2.6) MeV and m ρ = 830.5(6.3)(6.1) MeV, where the second error is due to the uncertainty of a −1 t . For most of this paper, we will not show the second error when listing our physical values. This means that most of the errors below are just statistical errors. We will bring the a −1 t error back to physical values in conclusion. Using the m π above, we have m π L s ≈ 3.9 for this lattice setup, which warrants the small finite volume effects. The number of configurations of our gauge ensemble is 6991, which is crucial for the glueballrelevant studies. The details of the gauge ensemble are given in Table II.
For the valence charm quark, we adopt the clover fermion action in Ref. [37], and the charm quark mass parameter is tuned to give (m ηc + 3m J/ψ )/4 = 3.069 GeV. With the tuned charm quark mass parameter, we generate the perambulators of charm quark on our ensemble, from which the masses of η c and J/ψ are derived precisely to be m ηc = 2.9750(3) GeV and m J/ψ = 3.0988(4) GeV. The according 1S hyperfine splitting is ∆ HFS = m J/ψ − m ηc = 123.8(5) MeV. We also check the dispersion relation in Eq. (1) for η c and J/ψ up to the momentum mode ⃗ n = (1, 2, 2). As shown in Fig. 2, the dispersion relation is almost perfectly satisfied with ξ = 5.341(2) and 5.307(5) for η c and J/ψ, respectively.
As a further check of our scale setting scheme, we also calculate the masses of D and D * and get m D  (7) GeV and m D * + − m D + = 0.14060 (7) GeV [5]. This manifests that our tuning of charm quark mass works well. The dispersion relation Eq. (1) is also checked to be correct for D and D * with ξ = 5.32(2) and 5.31 (3), respectively. The figure is similar to Fig. 1 and Fig. 2 and is omitted here to save space. Table III collects the results of η c , J/ψ, D and D * mesons. Along with the values of ξ derived from the dispersion relations of π and η, we can see that the values of ξ for different mesons are in agreement with the value ξ = 5.30 during the parameter tuning and are consistent with each other within 1%.

B. Operator construction and distillation method
The principal goal of this work is to investigate the possible mixing of the pseudoscalar glueball and the pseudoscalar qq meson, therefore the quark annihilation diagrams should be taken care of. For this to be done, we adopt the distillation method [33] which provides a smearing scheme for quark fields (Laplacian Heaviside smearing) and the calculation strategy of the all-toall propagators of the smeared quark fields that are distilled to be perambulators in the Laplacian Heaviside subspace of the spatial Laplacian operator −∇ 2 . Since we plan to investigate the η − η c correlation functions as well, we calculate the perambulators of u, d and c quarks on our large gauge ensemble in the Laplacian Heaviside space spanned by the N = 70 eigenvectors of −∇ 2 operator with the lowest eigenvalues. For the pseudoscalar glueball operator, we adopt the strategy in Refs. [14,15] to get the optimized Hermitian operator O G (t) = O † G (t) coupling mainly to the ground state glueball based on different prototypes of Wilson loops and gauge link smearing schemes (see Appendix for the details).
For the isoscalar η, the interpolation field can be defined as where Γ refers to γ 5 or γ 4 γ 5 , u (s) and d (s) are Laplacian Heaviside smeared u, d quark fields. Thus the correlation function of O Γ can be expressed as with C Γ (t) and D Γ (t) being the contributions from the connected and disconnected diagrams, respectively. We also consider the following correlation functions where the ∓ sign comes from the Hermiticity of O Γ , it takes minus sign for Γ = γ 5 (anti-Hermitian) and plus sign for is also defined in terms of the Laplacian Heaviside smeared charm quark field c (s) . Obviously, all of these correlation functions except for C GG (t) are contributed by quark annihilation diagrams and can be dealt with conveniently in the framework of the distillation method.

C. η mass as a further calibration
We calculate two types of correlation functions for η, namely C γ5γ5 (t) and C (γ4γ5)(γ4γ5) (t). We do observe the finite volume artifact that C γ5γ5 (t) approaches to a nonzero constant when t is large, as shown in Fig. 3. It has been argued that this constant term comes from the topology of QCD vacuum and can be approximately expressed as a 5 (χ top + Q 2 /V )/T where a is the lattice spacing (in the isotropic case), χ top is the topological susceptibility, Q is the topological charge, V is the spatial volume and T is the temporal extension of the lattice [30,38,39]. This can be understood from the anomalous axial vector current relation that the O γ5 has a direct connection with the topological charge density operator. It is interesting to observe that C (γ4γ5)(γ4γ5) (t) has normal large t behavior that it damps to zero for large t. As usual, the t-behavior of these correlation functions can be seen more clearly through their effective mass functions defined as where Γ refers to γ 5 or γ 4 γ 5 . Figure 4 shows these effective mass functions. Benefited from the large statistics of our gauge ensemble, the effective mass plateau starts from t ∼ 10, and the signal-to-noise ratio keeps good for t beyond 20 in the case Γ = γ 4 γ 5 . By a combined fit of C γ5γ5 (t) and C (γ4γ5)(γ4γ5) (t) using two mass terms, we obtain the best-fit mass of η to be where the error is obtained through the jackknife analysis. Considering the uncontrolled systematic uncertainties in our calculation (the continuum extrapolation and the chiral extrapolation have not been performed), this value is compatible with other determinations [28][29][30]. According to WV, the large mass of the flavor singlet pseudoscalar meson has a direct connection with the topology of QCD vacuum. In the N f = 2 case, to the leading order of chiral perturbation theory, m η is related to m π where the parameter m 2 0 is defined in terms of the topological susceptibility χ top and the pion decay constant f π , namely Usually, χ top is for the pure gauge case and is expected to be independent of flavor number N f . Using the values of m π = 348.5(1.0) MeV and m η in Eq. (9), m 2 0 can be derived to be m 2 0 = 0.3885(77) GeV 2 . According to the leading order chiral perturbation theory for N f = 2, the m π dependence of f π is [40] where m phys π = 134.98 MeV is the physical pion mass, ξ is approximated by m 2 π /(4πF ) 2 , the pion decay constant in the chiral limit F = 86(1) MeV and the low energy constant l 4 = 4.40 (28) are taken from FLAG 2019 [41].  [4]. In the physical N f = 3 case at physical pion mass, the GMOR relation implies the mass of the singlet counterpart of the pseudoscalar octet should be m 2 1 = (2m 2 K + m 2 π )/3 = 0.170 GeV 2 . Thus using the topological susceptibility we obtained, we can estimate the mass of the flavor singlet pseudoscalar meson to be which corresponds to m η1 ≈ 0.981(29) GeV and is not far from the experimental value m η ′ = 0.958 GeV. The GMOR relation also implies m 2 η8 = (4m 2 K − m 2 π )/3 ≈ 0.321 GeV 2 . This confirms again that the Witten-Veneziano mechanism for the flavor singlet pseudoscalar mass works fairly well both for the SU(2) and SU(3) flavor symmetry.

A. Theoretical consideration
In order to investigate the possible mixing between the pseudoscalar glueball and η, we must parametrize the correlation functions in Eq. (7). We adopt the following theoretical logic. As usual in the lattice study, the correlation function C XY (t) of operator O X and O Y can be parametrized as where the ± sign is for the same and opposite Hermiticities of O X and O Y , respectively, and |n⟩ are the eigenstates of the lattice HamiltonianĤ defined aŝ H|n⟩ = E n |n⟩ with E n being the corresponding eigen energies. For a given quantum number, |n⟩'s establish an orthogonal and complete set, namely, n |n⟩⟨n| = 1 with the normalization condition ⟨m|n⟩ = δ mn . In principle,Ĥ only exists heuristically, so we do not know the exact particle configurations of these eigenstates simply from the correlation function in Eq. (15). As far as the flavor singlet pseudoscalar channel is concerned in the N f = 2 QCD theory, each of the state |n⟩ should be a specific admixture of bare η states and bare glueballs if they exist theoretically (here we ignore the multihadron states temporarily) and can be taken as the states in the eigenstate set {|α n ⟩, n = 1, 2, . . .} of the free HamiltonianĤ 0 . Since we are working in a unitary lattice framework for u, d quarks, in principle this state set is orthogonal and complete with the normalization condition ⟨α m |α n ⟩ = δ mn . Now we introduce the interaction HamiltonianĤ I to account for the dynamics of the possible mixing, such that |n⟩ ofĤ =Ĥ 0 +Ĥ I can be expanded in terms of |α m ⟩ as with m |C nm | 2 = 1. In this sense, one can say that |n⟩ is an admixture of states |α m ⟩ whose fractions are |C nm | 2 , respectively. Furthermore, ifĤ I is small relative toĤ 0 , then to the lowest order of the perturbation theory, one has where E (0) n is the eigenenergy of |α n ⟩ and is ordered from low to high.
The experimentally observed isoscalar pseudoscalars are η, η ′ , η(1295), η(1405/1475), etc. [5], which are identified as I = 0 members of differentqq SU(3) nonets of different radial quantum numbers. As for the SU(2) case of this study, since there is only one isoscalar in each isospin quartet, the spectrum can be simplified largely. Because the smeared quark field suppresses the contribution of excited states in the correlation functions of η, as a simple approximation, we can truncate the spectrum of η state to be η and η * with η * taking account into all the excited η states. On the other hand, quenched lattice QCD predicted the mass of the lowest pseudoscalar glueball is around 2.4-2.6 GeV. This seems to be confirmed by the correlation function C GG (t) of the optimized operator O G , which is expected to couple most to the ground state. So we include the ground state pseudoscalar glueball |G⟩ and another state |G * ⟩ in the state basis {|α i ⟩, i = 1, 2, · · · } with |G * ⟩ standing for all the excited states of the pseudoscalar glueball. Finally, we have the following state basis |α i ⟩ = |η⟩, |η * ⟩, . . . , |G⟩, |G * ⟩, . . . , With this state basis, the free HamiltonianĤ 0 = diag{m η , m η * , m G , m G * } is diagonal with the matrix elements being the bare masses of the basis states, respectively, and being ordered from low to high. Theoretically, |η⟩ and |η * ⟩ are orthogonal, so do the states |G⟩ and |G * ⟩. Thus the interaction HamiltonianĤ I can be expressed as where x i , y i are called mixing energies sometimes. Accordingly, we have the following state expansion of |n⟩ The value of this correlator tends to zero at t = 0 within errors, and the error seems to be a constant independent of t. The correlator approaches a positive(negative) constant when t < T 2 (t > T 2 ) which might be due to the contribution from topology.

B. The Γ = γ5 case
Now we explore the possibility of glueball-η mixing through the correlation function C Gγ5 (t). Let us take a look at the t-dependence of C Gγ5 (t) shown in Fig. 5. We have the following observations: (1) C Gγ5 (t) is anti-symmetric with respect to t = T /2 and tends to 0 at t = 0. This is understood because O G is hermitian by construction (and even under the time reversal transformation T ) while O γ5 is anti-Hermitian and T -odd. At t = 0, since the product O G (0)O γ5 (0) is T -odd, its vacuum expectation value ⟨O G O γ5 (0)⟩ certainly vanishes.
(2) C Gγ5 (t) approaches a positive (negative) constant when t < T 2 (t > T 2 ). This may be due to the constant contribution from the topology similar to the case of C γ5γ5 (t) discussed in Sec. II C. Since C Gγ5 (t) is now T -odd, so does the topology contribution.
(3) Even though its central value is smooth, the error of C Gγ5 (t) is almost constant throughout the time range.
In order to find a function form to describe the time behavior of C Gγ5 (t), we take the following approximations by the assumptions ⟨0|O G |η i ⟩ ≈ 0 and ⟨0|O γ5 |G i ⟩ ≈ 0 similar to those adopted in the η − η ′ mixing studies [23,24].
Consequently, we have the following coupling matrix elements of the operators O G and O γ5 , and As an exploratory study and for the simplicity of the future data analysis, we temporarily neglect η * contributes to C Gγ5 (t). That is to say, we take the further approximation O † γ5 |0⟩ ≈ Z γ5,1 |η⟩. Thus after inserting Eqs. (20), (21), (22), (23) to Eq. (15) and ignore the terms relevant to |η * ⟩, we have the approximate expression of C Gγ5 (t) as The feature of this expression is C Gγ5 (t = 0) = 0 and is in accordance with the observation of item (1). In order to understand the almost constant error of C Gγ5 (t), we consider its variance [45] The first term on the right-hand side can be viewed as a correlation function of the operator O 2 γ5 and O 2 G , both of which have the vacuum quantum number 0 ++ (in the continuum limit) and are expected to have nonzero vacuum expectation values ⟨O 2 γ5 ⟩ ̸ = 0 and ⟨O 2 G ⟩ ̸ = 0. Thus we have The almost constant error of C Gγ5 (t) implies that the constant term ⟨O 2 γ5 ⟩⟨O 2 G ⟩ is large and dominate the variance. This is consistent with the argument in Ref. [45] that the variance of a correlation function is dominated by the possible lowest state, which corresponds to the vacuum state with E vac = The data points are lattice results and the shaded band illustrates the fitting results using Eq. (24). The signal to noise ratio of ∂tCGγ 5 is much better than that of CGγ 5 in Fig. 5. 0 in our case. This motivates us to consider the temporal derivative of C Gγ5 (t), namely, such that the constant term in C Gγ5 (t) and its constant variance can be canceled. This is surely the case. We plot ∂ t C Gγ5 (t) in Fig. 6, where one can see that ∂ t C Gγ5 (t) goes to zero when t is large and its relative error is much smaller.
In order for the mixing energies x 1 and y 1 to be extracted by using Eq. (24), one has to know the parameters m 1 , m 3 , m 4 , m η , m G , m G * , Z G , Z G * and Z γ5,1 , which, based on the assumptions of Eq. (21), are encoded in the correlation functions C γ5γ5 (t) and C GG (t) as where we take i = 1, 2 for C GG (t) and therefore G 1,2 refer to G and G * . It should be noted that the second state |G * ⟩ should be considered even though we have built the optimal operator O G based on a large operator set (seen in Appendix), it turns out that there is still a substantial contribution of higher states in C GG . To manifest this, we plot the effective mass function m eff (t) of C GG (t) in Fig. 7 by the definition in Eq. (8), where one can see that the ground state glueball G has not saturated C GG (t) before the signals are undermined by errors. With this observation, we add the second term relevant to G * to the first equation in Eq. (28) but do not consider its physical meaning. Since the mixing effect is expected to be small (our final results confirm this), we take the assumptions m 1 ≈ m η , m 3 ≈ m G and m 4 ≈ m G * . In the data analysis procedure, we first rearrange the N conf measurements into 139 bins with each bin including 50 measurements, and then we perform the one-eliminating jackknife analysis on these data bins. For each time of jackknife resampling, we first extract m η , m G , m G * , Z G , Z G * , Z γ4γ5,1 and Z γ5,1 from C GG (t), C γ5γ5 (t) and C (γ4γ5)(γ4γ5) (t) in the fixed time windows [t l , t h ] GG = [1,14], [t l , t h ] γ5 = [9,30] and [t l , t h ] γ4γ5 = [5,30], as shown in Table IV. Then we feed these parameters to ∂ t C Gγ5 (t) to determine the parameters x 1 and y 1 . The final results of m η , m G and |x 1 | with jackknife errors are obtained to be m η a t = 0.10358(84), m G a t = 0.3607(75), |x 1 |a t = 0.0155 (22). (29) Note that the definitions in Eqs. (21), (22), (23) are up to a plus or minus sign, we can only determine the absolute value |x 1 | of x 1 . The parameters of the fitting procedure and fitting results are collected in Table IV. The goodness of the fit of Eq. (24) to ∂ t C Gγ5 using the function is reflected by the χ 2 /dof = 0.96 in the fitting window [t l , t h ] Gγ5 = [3,30], and is also illustrated in Fig. 6 by the shaded curve. In mean time, making use of the Todd property of C Gγ5 (t), we average the t < T 2 part and t > T 2 part and find the errors can be reduced drastically around t = 0 except for C Gγ5 (t = 0), as shown in Fig. 8, where the function of Eq. (24) with fitted parameters is also plotted with shaded curves. It is seen that the function describes the t-dependence of C Gγ5 (t) very well up to an unknown constant term with opposite signs for t < T 2 and t > T 2 . Since |x 1 | is much smaller than m G − m η , to the lowest order of the perturbation theory, we can estimate the IV. Ground state mass and mixing angle fitted from operators with Γ = γ5 and Γ = γ4γ5 on the ensemble. Values of mη are the same in both cases because the value is derived from a combined fit to Cγ 5 γ 5 and C (γ 4 γ 5 )(γ 4 γ 5 ) . χ 2 are obtained from the fitting results of CGΓ.  C. The Γ = γ4γ5 case As a cross-check, we also carry a similar calculation by using the Γ = γ 4 γ 5 for the interpolation field operator of η. The corresponding correlation functions C (γ4γ5)(γ4γ5) (t) and C G(γ4γ5) (t) are calculated using Eq. (7). In contrast to the case of Γ = γ 5 , the correlation function C G(γ4γ5) (t) does not go to zero when t → 0. This is similar to the study of mixing of the pseudoscalar charmonium and the pseudoscalar glueball [32], which can be explained following the same logic that the QCD U A (1) anomaly may play an important role here. Obviously, the operator O γ4γ5 has the same operator structure as the temporal component of the isoscalar axial vector current j 5 µ (x) = 1 √ 2 ū(x)γ µ γ 5 u(x) +d(x)γ µ γ 5 d(x) , which satisfies the anomalous axial vector relation where j 5 (x) = 1 √ 2 ū(x)γ 5 u(x) +d(x)γ 5 d(x) is the pseudoscalar density and q(x) = g 2 32π 2 ϵ αβρσ G a αβ G a ρσ is the anomalous term stemming from the U A (1) anomaly with g being the strong coupling constant and G a αβ being the strength of color fields. Since j 5 (x) also has the same structure as O γ5 , based on the assumption in Eq. (21) we expect where f Gi is the decay constant of the pseudoscalar glueball G i . Accordingly we have Previous lattice studies show that pseudoscalar states can be accessed by the operator q(x) [15,16], thus the nonzero matrix element ⟨0|q(0)|G i ⟩ implies the coupling Z Gi |G i ⟩ still holds, then the correlation function C G(γ4γ5) (t) can be parametrized as which is similar to Eq. (24) apart from the first term due to nonzero coupling ⟨0|q(0)|G i ⟩ (here we ignore the contribution of excited η states). Alongside the correlation function and the first equation of Eq. (28), we carry out a similar jackknife data analysis to the case of Γ = γ 5 and get the mixing energy |x 1 |a t = 0.0112(55) and the corresponding mixing angle θ 1 = 2.5(1.2) • , which are consistent with the results of Γ = γ 5 case but have larger errors. The parameters of the fitting procedure and the fit results are also listed in Table IV for comparison. In Fig. 9, the colored curves with error bands are plotted using the fitted parameters. It is seen that the function forms mentioned above also describe the data very well. The fitted results are shown in Table IV and can be compared with the γ 5 -case directly. On this ensemble, it is clear that the results of the two cases are compatible with each other within errors.

IV. DISCUSSION AND SUMMARY
We generate a large gauge ensemble with N f = 2 degenerated u, d dynamical quarks on a 16 3 × 128 anisotropic lattice with the anisotropy ξ = a s /a t ≈ 5.3. Based on the experimental observation that the differences of mass squares of the heavy-light vector and pseudoscalar mesons are insensitive to quark masses, we propose a new scale setting scheme that the temporal lattice spacing a t can be estimated by this quantity. In practice, we use the value m 2 V − m 2 P S = 0.568(8) GeV 2 , which is a least squares fitting of nn, sn, ss, cn and cs mesons with n referring to the light u, d quarks. Our u, d quark mass parameter corresponds to a pion mass m π ≈ 350 MeV. We calculate the perambulators of u, d quarks and the valence charm quark on our ensemble using the distillation method where the quark annihilation diagrams can be calculated efficiently.
We calculate the mass of the isoscalar pseudoscalar meson η to be m η = 714 (6)(16) MeV. The error in the latter bracket is introduced by the estimation of lattice spacing in Sec. II A. This mass value can be matched, through the Witten-Veneziano relation, to the SU(3) flavor singlet pseudoscalar meson mass m η1 ≈ 981 MeV, which is in good agreement with the η ′ mass m η ′ = 958 MeV if the η − η ′ mixing is considered.
The mixing of η and the pseudoscalar glueball is investigated for the first time on the lattice. From the correlation function of the glueball operator and η operator, the mixing energy and the mixing angle are determined to be |x 1 | = 107(15)(2) MeV and |θ| = 3.46(46) • given the pseudoscalar glueball mass m G ≈ 2.5 GeV. This mixing angle is so tiny that the mixing effects can be neglected when the properties of η (and also η ′ in the physical N f = 3 case) are explored.
With perambulators of light and charm quarks at hand, we also checked correlation functions of η and η c , which are contributed only by annihilation diagrams. The Γ insertion are chosen to be γ 5 and γ 4 γ 5 . The signal-to-noise ratios of these correlation functions are fairly good, and the effective masses of them defined by Eq. (8) are presented in Fig 10, where the horizontal gray line illustrates the η mass in the lattice unit. Albeit the different behaviors at the small t range, all the effective masses approach m η when t is large. Because of the influence from excited states of η and η c , we cannot fit the curve by the simplification above, and optimized operators for these states are required to get better results. This kind of exploration can be performed in the future.
To summarize, we find that there is little mixing between the flavor singlet pseudoscalar (η is our case) and the pseudoscalar glueball. The topology of the QCD vacuum plays a crucial role in the origin of η mass due to the QCD U A (1) anomaly. This is in accordance with the common theoretical considerations [46].  [47] and QUDA library [48,49] are acknowledged. The computations were performed on the HPC clusters at Institute of High Energy Physics (Beijing) and China Spallation Neutron Source (Dongguan), and the CAS ORISE computing environment.

A1. Glueball operator construction
The interpolation field operators for the pseudoscalar glueball are built based on the four prototypes of Wilson loops shown in Fig. 11. The Wilson loops of each prototype are calculated using smeared gauge links. We adopt six different schemes to smear gauge links, which are different combinations of single link smearing and double link smearing [14,15]. The spatial symmetry group on the lattice is the octahedral group O whose irreducible representation A 1 corresponds to the spin J = 0, 4, . . . in the continuum limit. It is expected the glueball of J = 4 is higher than the J = 0 state in mass, so we can build operators in the A −+ 1 representation to couple with the ground state pseudoscalar glueball. Let W α (x, t) be one prototype of the Wilson loop under a specific smearing scheme, then the A −+ 1 operator in the rest frame of a glueball can be obtained by (36) where R • W α refers to differently oriented Wilson loops after one of the 24 elements of O (R) operated on W α , P is spatial reflection operation and c A1 R are the combinational coefficients for the A 1 representation. Thus we obtain a A −+ 1 operator set {ϕ α (t), α = 1, 2, . . . , 24} based on the four prototypes and six smearing schemes. In order to get an optimized operator O G that couples most to the ground state glueball, we first calculate the correlation matrix of the operator set, and then solve the generalized eigenvalue problem C αβ (t 1 )w β = λC αβ (t 0 )w β (38) to get the eigenvector w α of the largest eigenvalue λ, which serves as the combinational coefficients of O G , namely, In practice, we take t 0 = 1 and t 1 = 2 for the optimized operator coupling more to the ground state pseudoscalar glueball.