Revisiting the pseudoscalar meson and glueball mixing and key issues in the search for pseudscalar glueball state

We revisit the mixing mechanism for pesudscalar mesons and glueball which is introduced by the axial vector anomaly. We demonstrate that the physical mass of the pseudoscalar glueball does not favor to be lower than 1.8 GeV if all the parameters are reasonably constrained. This conclusion, on the one hand, can accommodate the pseudoscalar glueball mass calculated by Lattice QCD, and on the other hand, is consistent with the high-statistics analyses at BESIII that all the available measurements do not support the presence of two closely overlapping pseudoscalar states in any exclusive channel. Such a result is in agreement with the recent claim that the slightly shifted peak positions for two possible states $\eta(1405)$ and $\eta(1475)$ observed in different channels are actually originated from one single state with the triangle singularity interferences. By resolving this long-standing paradox, one should pay more attention to higher mass region for the purpose of searching for the pseudoscalar glueball candidate.


I. INTRODUCTION
The non-Abelian property of Quantum Chromo-Dynamics (QCD) predicts the possible existence of glueball states as a peculiar manifestation of the strong interaction in the non-perturbative regime. However, until now indisputable experimental evidence for the glueball states is still lacking. In the pseudoscalar sector, the flux tube model supports a low-lying pseudoscalar glueball with a mass around 1.4 GeV [1]. This was the mass region accessible by several experiments in 1980's and 1990's, for instance, Mark-III [2,3], DM-2 [4,5], OBELIX [6][7][8], and BES-II [9]. Reviews on the early experimental observations can be found in Refs. [10,11]. With the strong motivation of looking for glueball candidates in experiment, the observation of three possible pseudoscalar states with isospin 0 around 1.3∼ 1.5 GeV, i.e. η(1295), η(1405), η(1475), was regarded as the clues for the presence of a pseudoscalar glueball in association with the isospin singlets in the qq scenario. Note that there have been well-established states, i.e. π(1300) and K(1460), in the same mass region with which the first radial excitation of the qq pseudoscalar meson nonet with J P = 0 − can be formed [11][12][13]. For a long time following the rather vague experimental results, there have been tremendous efforts trying to understand the property of these three states among which the η(1405) has been assigned as the most-likely pseudoscalar glueball candidate. Other explanations for the out-numbering of isoscalar pseudoscalar states around 1.3 ∼ 1.5 GeV include dynamically generated states [14] and tetraquarks [15]. However, any explanation for the out-numbering problem should first confirm whether indeed an additional state is present.
The phenomenological studies of the pseudoscalar glueball candidate η(1405) have been focussed on the following main issues: • Whether there are mixings among the ground state pseudoscalar mesons η and η ′ , and the pseudoscalar glueball?
Following the questions from item one, most studies assume certain mixing mechanisms among η, η ′ and η(1405) and investigate the properties of η(1405) in gluon-rich processes such as J/ψ radiative decays. Also, the gluon contents inside η and η ′ can provide some hints of glueball states due to the mixing mechanism. As a consequence of such a mixing, one expects that observable effects can be measured in experiment which can make the glueball state different from the qq mesons. However, it is still difficult to conclude that the pseudoscalar glueball state has been observed in experiment taking into account the high precision measurements from the BESIII experiment and LQCD calculations. This is related to the questions raised above in items two and three.
During the past decade the progress of LQCD has brought many novel insights into the light hadron spectroscopy via numerical simulations of the non-perturbative strong interactions. Interesting and surprisingly, it shows that the lightest pseudoscalar glueball should have a mass around 2.4∼2. 6 GeV in a quenched calculation [22][23][24][25], while the later dynamical calculations [26,27] suggest that the mass of the lightest pseudoscalar glueball does not change much compared with the quenched result. This is obviously in contradiction with the data if η(1405) is assigned as a glueball candidate.
In parallel with the LQCD studies, great efforts have been made in experiment in order to establish η(1405) as an additional state apart from η(1295) and η(1475). A natural expectation is that since both η(1405) and η(1475) have the same quantum number and can couple to the same hadronic final states there should be channels that they both can have observable couplings, hence, nontrivial structures caused by two closely overlapping and interfering states should appear in the mass spectra. However, with the high-statistics measurements in various channels, e.g. in J/ψ radiative and hadronic decays at BESIII [32][33][34], there is no any evidence indicating that two nearby η(1405) and η(1475) have been produced together in the same channel. All the data so far only show one peak structure around 1.42 GeV and no need to introduce interfering states from two nearby states. These new measurements actually have brought serious questions on the need for an additional η(1405) apart from the radial excitation of the qq isoscalars η(1295) and η(1475).
The new data also raise new features for the radial excitation spectrum of the isoscalar states η and η ′ . One notices that the single peak positions for η(1405/1475) are slightly shifted in different channels. In particular, the observation of the significantly large isospin breaking effects in J/ψ → γη(1405/1475) → γ + 3π can be regarded as an indication of a special mechanism that causes the mysterious phenomena around 1.4∼1.5 GeV for the isoscalar pseudoscalar meson spectrum [32]. It was proposed by Refs. [29,30] that the presence of the so-called "triangle singularity (TS)" mechanism can enhance the isospin breaking effects and shift the peak positions of a single state by the interferences in exclusive decay channels. Similar analysis of Ref. [31] also confirms that the TS contribution is needed in order to understand the strong isospin breaking effects.
The TS mechanism was first investigated by Landau in 1950's [35] and followed up by many detailed studies later [36][37][38][39][40][41][42]. It states that for an initial state with energies near an intermediate open threshold, if the rescattering between these two intermediate states by exchanging another state (i.e. via a triangle diagram) into three-body final states would allow such kinematics that all the three intermediate states can approach their on-shell condition simultaneously, to be located within the physical region, then the triangle loop amplitude will be enhanced by the threebody singularity as the leading contribution. As a consequence, its interference with the tree-level transition amplitude of the initial state can shift the its peak position and even change the lineshape [29]. In the case of η(1405/1475) → 3π the mass of the initial state η(1405/1475) is within the TS kinematic region and has strong couplings to KK * + c.c. Thus, the intermediate KK * + c.c. and the exchanged kaon in the triangle loop can approach the on-shell condition simultaneously and results in the strong enhancement of the isospin breaking on top of the a 0 (980) and f 0 (980) mixing. The recognition of the TS mechanism here provides an alternative explanation for understanding the η(1405)-η(1475) puzzle and can resolve the contradiction between the LQCD results and experimental observations for the pseudoscalar glueball. Recent detailed analyses and discussions on the TS mechanism can be found in Refs. [43]. More recognitions of this special kinematic effects in various processes can be found in the literature [44][45][46][47][48][49][50][51][52][53][54][55][56][57][58] and recent reviews [28,59].
The above progress suggests that the pseudoscalar meson and glueball mixing mechanism should be re-investigated. Moreover, given that the pseudoscalar glueball mass in the quenched approximation is around 2.4∼2.6 GeV, its mixing with the cc(0 −+ ) should also be considered. An earlier study of the mixing mechanism has implemented the anomalous Ward identities with the corresponding equations of motion which connect the transition matrix elements of vacuum to η, η ′ and glueball to the pseudoscalar densities and the U(1) anomaly [60]. There, the physical glueball state was assigned to η(1405) and then the mixing effects on η, η ′ and η c were studied. Due to a large number of parameters in the mixing scheme of Ref. [60], it shows that a re-investigation of the parameter space is necessary. In particular, a detailed analysis of the sensitivity of the glueball mass range to the mixing parameters is necessary. This will help further clarify the puzzling situation around 1.4 GeV in the I = 0 pseudoscalar spectrum. One also notices that the recent analysis of Ref. [61] in a chiral Lagrangian approach with an axial anomaly coupling also leads to a much heavier glueball mass than the range of around 1.4 GeV.
As follows, we first introduce the formulation of the mixing scheme via the axial vector anomaly as studied in Refs. [62,63] in Sec. II. We then inspect the parameter space and impose constraints on these parameters in order to investigate the mass range of the pseudoscalar glueball. In particular, the sensitivities of the physical glueball mass to the parameters will be scrutinized. In Sec. III, the numerical results are presented and discussed. A brief summary is given in Sec. IV.
This work is organized as follows: In Sec. II we first introduce the formulation of the mixing scheme via the axial vector anomaly as studied in Refs. [62,63]. We inspect the parameter space and try to investigate the sensitivities of the physical glueball mass to the parameters. In Sec. III, the numerical results are presented and discussed. A brief summary is given in Sec. IV.

II. THE MIXING FORMALISM
A. η − η ′ − G − ηc mixing scheme As stated in Ref. [62,63], the well-known axial vector anomaly is, where j denotes the q, s, c quark respectively, and m j denotes the quark masses, G andG denote the strength tensor and the dual of the gluon field. The physical states are mixture of the pure states via a unitary matrix U as, where |η q , |η s , |g , |η Q denote |qq ≡ |(uū + dd)/ √ 2 , |ss , the unmixed glueball state, and the unmixed heavy quark state |cc .
Assuming that the decay constants in the flavor basis follow the same mixing pattern of the particle states [62], we have where all the OZI-suppressed off-diagonal elements are neglected. The pseudscalar meson decay constants are defined as follows, where M P is the diagonal mass matrix of the physical states that is explicitly written as, Noted that the meson state is with the dimension of mass −1 and the decay constant with the dimension of mass. Based on the above definitions and assumption, we can obtain the mass matrix on the flavor basis from two ways. On the one hand, the mass matrix on the flavor basis is related to the physical particle mass via a unitary transformation. On the other hand, according to the definition of the decay constants, the mass matrix is also related to the axial vector current divergences in a more dynamical and explicit way. Although some of the matrix elements cannot be well constrained and determined quantitatively, they are not going to affect our discussions here due to their small values that can be qualitatively determined. The mass matrix in terms of the physical masses can be written as, In order to obtain the mass matrix in terms of the divergences of the axial vector current, we firstly define the following abbreviations for pseudoscalar densities and the U(1) anomaly matrix elements as done in Ref. [60]: qq,qs,qg,qc ≡ √ 2 f q 0|m uū iγ 5 u + m dd iγ 5 d|η q , η s , g, η Q , m 2 sq,ss,sg,sc ≡ 2 f s 0|m ss iγ 5 s|η q , η s , g, η Q , Note that the definition for q (u, d) quark current is different from other quark flavors by a factor of √ 2 due to the definition of the |qq .
Then, the mass density matrix of the q, s, c dimension can be written explicitly in a dynamical way as in Ref. [60], The mass density matrix obtained from these two ways should be the same. Thus, the mixing information could be revealed.
To proceed, we first analyze the parameters involved in this mixing scheme by looking at the transformation matrix U . In Ref. [64], a general form for the unitary mixing matrix is presented with six independent rotation angles. It would not be realistic to determine all of them based on what we know about the pseudoscalar meson and glueball mixing. In order to implement constraints on the mixing matrix elements, we take a similar strategy of Ref. [60] to reduce the number of parameters.
Firstly, the mixing between the light favor octet state η 8 and glueball is neglected in the SU(3) flavor symmetry. Secondly, the heavy-flavor state mixing with the light-flavor state is also neglected, since they have a large mass difference and is OZI suppressed. These will reduce the number of undetermined parameters to only three mixing angles, i.e. the mixing angle φ Q between the heavy-flavor state and glueball, the mixing angle φ G for the glueball and light-flavor singlet state mixing, and the mixing angle θ between the octet and singlet light flavor states that mainly determines the structure of η and η ′ . So the mixing matrix between the flavor states and the physical states can be written as [60], where c and s are the short-handed notations for "cos" and "sin"; θ i are the ideal mixing angle between η q ≡ (uū + dd)/ √ 2 and η s ≡ ss. The mass density matrix element from the physical state mass through the U matrix can be obtained. The explicit expressions for each matrix element can be found in Ref. [60]. Here, we concentrate on the matrix elements that are relevant in the extraction of physical quantities of interest.

B. Constrain the parameters
The mixing mechanism discussed in Ref. [60] and summarized above allows us to express the mass matrix as follows: On the left hand side of the equation, there are four parameters, i.e. the physical glueball mass M G and three mixing angles. On the right hand side, more parameters emerge which are related to the mixing dynamics. Apart from f q , f s , f c , m cc , m qq , m ss that are more explicit and can be estimated phenomenologically by observable physical quantities, there are still nine pseudoscalar densities and four U(1) anomaly matrix elements to be determined. Note that parameters m qq and m ss are related to the relatively well defined η and η ′ mixing, they can be extracted from M 11 qsgc andM 22 qsgc in Ref. [60]. Besides, m qq is too small and in some cases the result even flips the sign [16]. Actually, since the masses of η and η ′ are rather far away from the glueball mass, the mixing effects due to the presence of glueball are expected to be small. In this sense, the constraint from the glueball contents of η and η ′ could be still marginal.
As discussed in Ref. [16], the OZI-violating light-flavor pseudoscalar density m qs , m sq scales as O(1/N c ) in the limit of the large color number N c , and m qg is of the order higher than m qq which is as small as m 2 π . Thus, we can drop these three parameters in this analysis and this is different from the treatment of Ref. [60]. We will show later that this is a reasonable assumption.
The flavor mixing angle θ for η and η ′ are constrained in a finite range −17 • < θ < −11 • . Although θ is not precisely fixed, its influence on the glueball property is rather small. A relatively small mixing angle θ < −10 • is favored by the unquenched LQCD calculation [65]. Thus, we adopt θ = −11 • which is the same as in Ref. [60].
The glueball component within the physical η c can be estimated by an empirical gluon power counting rule [66] to combine with the experimental data of the branching ratio of η c → γγ, as done in Ref. [60]. With the updated experimental data BR(η c → γγ) = (1.59 ± 0.13) × 10 −4 [12], φ Q = −2.7 • and 11.6 • are obtained with the central value. As discussed in Ref. [60], a negative φ Q is not favored by the radiative decays of J/ψ, ψ ′ → γη c . Therefore, we adopt the positive value of φ Q = 11.6 • .
The last mixing angle φ G is directly related to the physical glueball mass as shown by the third row of the mass density matrix (Eq. (10)). Thus, a reliable determination of this quantity is crucial for estimating the physical glueball mass range.
The above consideration has significantly reduced the parameter number but still there are more than 10 parameters to be determined in Eq. (10). In Refs. [16,60], a different treatment for the parameters was applied to estimate the glueball mass. By taking the ratio of elements, e.g.M 31 qsgc /M 32 qsgc in Eq. (10), and assuming the negligibly small values of m 2 qg and m 2 sg compared with √ 2G g /f q and G g /f s , the glueball mass will depend on the ratio of f s /f q , while its dependence on G g will be cancelled. A caveat of this treatment is that m 2 sg actually is not small enough to be neglected. This point will be discussed later. On the other hand, if m 2 qg and m 2 sg are neglected, it will lead to independence of the glueball mass on parameters G q,s,g,c . However, since G q,s,g,c describes the contributions from the pseudoscalar U(1) anomaly in Eq. (10), one would expect its direct connection with the physical glueball mass in the constraint relation. Our revisit to this issue is to examine how the glueball mass should depend on G q,s,g,c in an explicit way.
Still focussing onM 31 qsgc /M 32 qsgc in Eq. (10), we extend the discussions on the parameters slightly. These two elements have the following expressions: Note that in Eq. (11) it is safe to neglect m 2 qg and only keep term √ 2G g /f q since m 2 qg ≪ m 2 qq with m 2 qq about 36 times smaller than √ 2G g /f q . However, it is not obvious to neglect m 2 sg in Eq. (12) since so far we only know the relation of m sg ≪ m ss [16], but have no information about the values of m sg . Similar situation occurs with m 2 cg and m 2 cc when treating the elementsM 41 qsgc andM 42 qsgc , i.e. the only known information is m cg ≪ m cc . Apparently, if the value of m 2 sg is compatible with G g /f s , it will result in large uncertainties when taking the ratio ofM 31 qsgc /M 32 qsgc . Meanwhile, the sensitivities of the glueball mass to G g will be lost. This problem can be seen more clearly if one compares the following two equal ratios extracted from Eq. (10): where Eq. (13) leads toR 31/32 ≃ √ 2f s /f q after neglecting m 2 qg and m 2 sg . However, note that √ 2|G c /f q | ≃ 0.039 GeV 2 and G c /f s ≃ 0.023 GeV 2 both are much smaller than |m 2 qc | = 1.197 GeV 2 and |m 2 sc | = 0.092 GeV 2 in Eq. (14). The neglect of m 2 qg and m 2 sg in Eq. (13) and m 2 qc and m 2 sc in Eq. (14) cannot be justified. Thus, although the equivalenceR 31/32 =R 41/42 can be deduced rigorously from the left-hand side of Eq. (10), the relation of R 31/32 =R 41/42 ≃ √ 2f s /f q actually does not hold. To proceed, we take a slightly different strategy to determine the parameters and extract the pseudoscalar glueball mass. First, it should be noted that an explicit relation between G g and the glueball mass should be retained. Note that the value |G g | = (0.054 ± 0.008) GeV 3 , has been calculated by lattice QCD in the quenched approximation [22]. Although the sign of G g is not determined by LQCD, we will show that the positive value can be excluded since it will lead to negative values for the glueball mass. We also take the decay constant f q as an input. It is relatively well constrained to be (1 ∼ 1.1)f π while f s varies within a range of (1.3 ∼ 1.6)f π [62,63,67,68]. Another two parameters that we adopt are f c = 487.4 MeV [64] and m cc ≈ M 2 ηc [62]. These are reasonable approximations taking into account the success of potential model in the description of low-lying charmonium states. Note that the decay constant f J/ψ = 405 MeV [69] and the quenched mass m ηc = 3.024 GeV [70] are provided by LQCD. As the leading approximation we assume that the cc bare vector (J/ψ) and bare pseudoscalar (η c ) share the same wave function at the origin as expected in the heavy quark spin symmetry (HQSS) limit, although in reality the HQSS breaking effects cannot be neglected. Following the same reason, it is reasonable to adopt the physical η c mass for m cc in contrast with the LQCD quenched mass m ηc = 3.024 GeV.
With the above parameters fixed we are left with 12 equations with 13 undetermined parameters, i.e. M G , φ G , G q , G s , G c , m sg , m qc , m sc , m cq , m cs , m cg , m qq , m ss . Note that as mentioned earlier, m 2 sq , m 2 qs , m 2 qg are neglected since m 2 qs,sq ≪ m 2 qg ≪ m 2 qq with m 2 qq about 36 times smaller than √ 2G g /f q . Therefore, we make the approximation to Eq. (11) which leads to One notices that not all the parameters are explicitly correlated in a single relation in this mixing scheme. This allows us to investigate the relation between two unknown quantities in a single equation while the other parameters can be fixed with reasonable values. Following this consideration, Eq. (15) can be approximated by where all the cosine values of the small angles have been taken as unity, and the glueball mass sensitivity to φ G can be investigated. Note that φ G is not well constrained and its value varies in a wide range, depending on the parametrization of the mixing matrix, experimental inputs and fitting procedures [16]. For example, (12 ± 13) • is obtained in the radiative decay of V → P γ, P → V γ independent of f q and f s [67]. The variation range of φ G was found to be (32 +11 −22 ) • in the strong process J/ψ → V P in Ref. [71]. Similar scenario was also studied in Refs. [18,72]. In the next Section we provide relations for the other parameters in terms of φ G and the ratios f c /f q , f s /f q , and f s /f c . Although these three ratios are not independent, the idea is to investigate whether those undetermined parameters can have acceptable values located within a common regime of φ G and the ratios.
Note that for the determination of the pseudoscalar glueball mass via Eq. (11), the decay constant f s is not explicitly involved. However, its correlation with other parameters will affect the result to some extent, in particular, via the ratio of f s /f q . Note that f q ≃ (1 ∼ 1.1)f π is relatively well determined while f s ≃ (1.3 ∼ 1.6)f π is well estimated [62,63,67,68]. In order to determine all the parameters self-consistently, we will fix the ratio of f s /f q with commonly accepted values and solve the 12 equations with 12 parameters. The ratio of f s /f q contains the uncertainties of SU(3) symmetry breaking effect. We will show later that the uncertainties arising from the SU(3) symmetry breaking will not change the magnitude hierarchy of the correlated parameters. In particular, we will see that the calculated glueball mass should not be sensitive to the ratio f s /f q which is different from the result of Refs. [16,60] due to different ways of treating the parameters.

A. Pseudoscalar glueball mass and its correlations with other parameters
We first study the relation between M G and φ G in Eq. (15). The mixing angle between the flavor singlet and octet states is fixed as θ = −11 • , and the mixing angle between the pure glueball and the pure heavy quark state is fixed as φ Q = 11.6 • . One can check that the physical glueball mass M G keeps stable within the reasonable ranges of φ Q and θ. We adopt f q = f π = 131 MeV as an input. Parameter G g is fixed as |G g | = (0.054 ± 0.008) GeV 3 from the quenched LQCD calculation. As mentioned earlier, the positive G g is excluded in our model since it will result in negative values for M 2 G . This is consistent with analyses of Ref. [16,60]. Actually, within the favored space for all the other parameters the negative values for G g are always required. Therefore, we fix G g = −(0.054 ± 0.008) GeV 3 in this analysis and the φ G dependence of M G can be investigated.
As mentioned earlier, the light quark and glueball mixing angle φ G has a relatively large variation range, i.e. φ G ∈ (3, 25) • , we can then investigate the dependence of other quantities on the φ G within the range of ∈ (3, 25) • . Note that φ G = 0 corresponds to a vanishing mixing between light quarks and glueball. It simply means that as long as the mixing is introduced, the mixing angle will be constrained by other quantities (G g in this case) and deviate from zero. In Fig. 1 we present the physical glueball mass in terms of φ G with G g = −(0.054 ± 0.008) GeV 3 . It shows that M G is very sensitive to φ G . With a small value of φ G = 3 • , a large glueball mass of about 3.8 GeV can be extracted. This is even much larger than the pure gauge glueball mass and η c mass. So a small φ G like 3 • is certainly unphysical. As indicated by the central value of φ G = 12 • from a model analysis [67], the glueball mass is found to be M G ∈ (2.0, 2.2) GeV where the uncertainties are given by the uncertainties of G g . This is the range which is not far away from the pure gauge glueball mass by LQCD. It is worth quoting the unquenched LQCD calculations for the pseudoscalar glueball mass in the literature. For instance, the UKQCD Collaboration reported M G ≃ 2.5 ∼ 2.7 GeV at m π = 280 and 360 MeV, respectively [27], and M G = 2.56 ∼ 2.60 GeV at m π = 938 ∼ 650 MeV were also found by Ref. [26]. Although these results are extracted at relatively high pion mass region, it is very much unlikely that the physical state (a P -wave gluonic state) should have a mass lower than that for the scalar glueball (a S-wave gluonic state), i.e. around 1.5∼ 1.7 GeV. Our analysis also supports such a scenario. The results in Fig. 1 suggest that low glueball masses, e.g. lower than 1.8 GeV, cannot be accommodated by the mixing mechanism via the axial vector anomaly. As shown by Fig. 1, even for a much larger and unrealistic value of the mixing angle, the glueball mass will be still higher than 1.5 GeV. This eventually rules out the possibility of a light pseudoscalar glueball around 1.4 GeV.
By substitute Eq. (15) into the mass density matrix Eq. (6), we obtain the explicit expressions for the mass density matrix elements as follows, There are apparent features arising from the mixings described by the above equation array. In the light flavor sector the mixing between |qq and |ss is dominated by the flavor singlet and octet mixing as expected. The glueball mixing contributions are at order of G g /f q but it enters into the light quark submatrix with a suppression factor tan φ G .
Due to the small value of φ G , one would not expect significant contributions from the glueball mixings in η and η ′ as demonstrated by many studies. The mixing between the heavy and light flavors can be seen in M 14,24 qsgc , which is at order of tan φ Q and can be neglected. This is anticipated due to the large mass difference between η c and η (η ′ ). The mixing between glueball and heavy flavor cc can be seen from M 34 qsgc , where cancellations among the terms are present. Furthermore, this element is proportional to sin 2φ Q . Given the small value of φ Q , this factor also imposes a suppression to the mixing effects.
In the limit of small values for φ G and φ Q , the element M 33 qsgc can be approximated as where we keep the correction from η c to show the suppressed contributions from the heavy flavor part. The last line is obtained by substituting Eq. (16) into the equation with sin θ i = 2/3. It is interesting to compare the above expression with Eq. (16). It shows that the mixing effects on the mass of glueball from the quark states are indeed suppressed. Apart from the term of M 2 ηc sin 2 φ Q from the heavy flavor mixing, corrections from the light flavor singlet and octet mixings will introduce cancellations. The numerical calculation indeed suggests that such mixing effects cannot significantly change the pure glueball mass. Alternatively, it implies that the physical η c will have subleading mixing contributions from the glueball. This feature can be seen by the element M 44 qsgc , and is consistent with the experimental observations. A detailed investigation of this aspect has been presented in Ref. [60].
Combining the mass density matrix elements in Eq. (17) with those defined in Eq. (8), all the unknown parameters can be written in terms of φ G and the constrained parameters as follows, 2f c (sin 2θ cos 2θ i cos φ G + sin 2θ cos 2 θ i sin φ G tan φ G + cos 2θ sin 2θ i One notices that the elements between the light and heavy flavor mixings are suppressed explicitly either by the mixing angle or a term of (M 2 ηc cos 2 φ Q − m 2 cc ). This is understandable due to the large mass differences between η c and η (η ′ ). To see more clearly the dependence of the mixing matrix elements on the mixing angles and decay constants, we substitute the values of those fixed parameters, i.e. M η , M η ′ , M ηc , θ i , and G g , into the above equations. The mixing elements will be explicit functions of θ, φ G , φ Q and the decay constants. We can then investigate their relations by numerical calculations.
− f s f q (0.41 cos 2θ − 0.14 sin 2θ cos φ G (1 − tan 2 φ G )) , In the above equation array M G shows explicit dependence on G g and φ G . G q , G s and G c are proportional to f s , f q and f c , respectively. One notices that G c , m 2 qc , and m 2 sc contains the large cancellation term (M 2 ηc cos 2 φ Q − m 2 cc ). Thus, G c turns out to be sensitive to m 2 cc and φ Q . These quantities are also dependent on φ G due to the factor cot φ G there. There are also large cancellations in m 2 qq which is sensitive to both f s /f q and θ. In contrast, the cancellation in m 2 ss is relatively small, and it shows small dependence on f s /f q . m 2 sg shows sensitivities to f s /f q since it contains an SU(3) flavor symmetry breaking factor (f s /f q − 1). A cancellation also occurs in m 2 cg , and m 2 cg increases with the decreasing φ G . In contrast, it decreases when the φ Q gets smaller. In Eq. (20), we have also listed m 2 cq and m 2 cs in terms of φ Q , m 2 cc and φ G in comparison with m 2 qc and m 2 sc . However, they show little dependence on all these mixing angles.
By adopting G g = −0.054 GeV 3 , φ Q = 11.6 • , φ G = 12 • , f c = 487.4 MeV, f q = 131 MeV, and taking the ratio f s /f q = 1.2 and 1.3 in order to examine the sensitivity of the SU(3) flavor symmetry breaking, we can determine all the other quantities and they are listed in Table I for the two ratios of f s /f q , respectively. One notices that some of these parameters do not explicitly depend on f s as discussed above. Thus, they do not change values when taking different ratios for f s /f q .
The correlations among the parameters should be further discussed. In Table I it shows that m 2 qq is quite sensitive to the SU(3) flavor symmetry breaking ratio f s /f q . Namely, m 2 qq changes order with the ratio f s /f q varying from 1.2 to 1.3. In contrast, m 2 sg changes by about a factor of 1.5. Such a dependence can be seen from Eq. (19). Taking the where significant cancellations occur with the increasing ratio of f s /f q . With the dominance of the constant term the value of m 2 qq will decrease due to a cancellation caused by the increasing ratio of f s /f q . Similar phenomenon happens to m 2 sg . It is reasonable that f q and f s would affect the light quark mass term that are related to the light flavor states η q and η s . Since the other parameters keep stable with the varying ratio of f s /f q in a reasonable range, we will focus on the results with f s /f q = 1.2 in the following discussions.

B. Extracting topological susceptibilities for pseudoscalar mesons
It is noticeable that the anomaly matrix elements G q , G s and G c are of the same order of G g . G q and G s are quite large which is an indication of the important role played by the anomaly term in the U(1) Goldstone boson [73,74]. As also discussed in Ref. [16], the anomaly terms 0|α s GG/(4π)|η and 0|α s GG/(4π)|η ′ could be related to the topological susceptibility. In our scheme we obtain, where the central values of φ G = 12 • and m 2 cc = M 2 ηc are adopted. These results are consistent with the LQCD calculations, namely, 0|α s GG/(4π)|η ≈ 0.021 GeV 3 [75], 0|α s GG/(4π)|η ′ ≈ 0.035 GeV 3 [76], which have been determined in the chiral limit on LQCD, and G g = −(0.054 ± 0.008)GeV 3 calculated in the quenched approximation [22].
In order to estimate the uncertainties arising from the parameter ranges, we plot the topological susceptibility G P ≡ 0|α s GG/(4π)|P , with P stands for the physical states η, η ′ , G and η c , in terms of m 2 cc , φ Q and φ G in Fig. 2. We consider the dependence of G P on these three quantities in three cases, i.e. (a) G P dependence on m cc (with φ G = 12 • , φ Q = 11.6 • fixed); (b) G P dependence on φ G (with m 2 cc = M 2 ηc , φ Q = 11.6 • fixed); and (c) G P dependence on φ Q (with φ G = 12 • , m 2 cc = M 2 ηc fixed). It shows that G η and G η ′ are not sensitive to m cc , φ G and φ Q mainly because the mixing between η q,s and η Q are small. More significant sensitivities of G G and G ηc indicate the non-negligible effects arising from the mixing between |η Q and |g .
To be specific, in Fig. 2 (a) G ηc is the only one sensitive to the value of m 2 cc due to the presence of the dominant cancellation term (8.9 cos 2 φ Q − m 2 cc ) as shown in Eq. (20). By adopting G c = 0|α s GG/(4π)|η c = −0.079 GeV 3 , the corresponding value of m 2 cc ) becomes close to M 2 ηc as expected. In contrast, in Fig. 2 (b) when fixing m 2 cc = M 2 ηc and φ Q = 11.6 • the φ G dependence of both G ηc and G G are sensitive in the small value range of φ G 10 • . This suggests that φ G can be well constrained by the mixing angle in our scenario. In Fig. 2 (c), by fixing φ G = 12 • and m 2 cc = M 2 ηc , all the quantities appear to be stable in terms of φ Q in a relatively broad range. The favored value φ Q = 11.6 • corresponds to an overall reasonably good description of all the other quantities.

C. Constraints on the charmonium state
In Table I m 2 qq and m 2 ss are of the typical values as those extracted in Refs. [16,60]. But other mass terms are rather different. As emphasized earlier, our strategy of fixing these parameters is to retain the mass hierarchy |m 2 cc | ≫ |m 2 cg | ≫ |m 2 cq,cs |, and m 2 ss ≫ m 2 sg . However, m 2 qc and m 2 sc , labelled with " * " in Table I, seem to be where elements U 41 and U 42 are directly dropped as they are treated as small quantities. We now bring back these two elements by defining U 41 = x and U 42 = y, and investigate the effects on M 41 qsgc and M 42 qsgc due to the nonvanishing x and y. The mass matrix M qsgc with x and y as explicit parameters has the following expression, (24) If we assumed that x, y have the same order of magnitude but an opposite sign to the corresponding symmetry matrix elements in the U matrix in Eq. (23), then both x, y should be order of −0.01. This is truly very small compared to the other elements. Meanwhile, the two elements M 41 qsgc and M 42 qsgc can be expressed as, and The above two equations suggest that to keep both M 41 qsgc and M 42 qsgc small it requires intrinsic dynamic constraints on both m 2 qc and m 2 sc of which the effects will then show up via x and y in the U matrix. The dependence of m 2 qc and m 2 sc on m 2 cc , φ G and φ Q are encoded via Eq. (19). It means that more stringent constraints on m 2 cc , φ G and φ Q should be applied in order to keep both m 2 qc and m 2 sc small. Interestingly, by requiring |m qc , m sc | ≤ 0.25 GeV 2 , we find that m 2 cc should be restricted within (8.7 ∼ 8.8) GeV 2 as shown by Fig. 3. It corresponds to (M ηc − m cc ) = 13.5 ∼ 30.4 MeV which means that the glueball-cc mixing has resulted in a mass gap between the physical and pure states. Meanwhile, it shows that the glueball-cc mixing does not change the main character of η c as the ground state pseudoscalar charmonium. However, due to the mixing, some of the observables may have indicated effects arising from a small glueball component in the wavefunction of η c . This is consistent with the conclusion of Ref. [60].
To show the sensitivities of m 2 qc , m 2 sc and G c /f q to m cc , we set the value of m 2 cc to be 100 MeV 2 below the mass of η c squared, i.e. m 2 cc = M 2 ηc − 100 MeV 2 . As shown by Fig. 4, φ G will be restricted within (7,12) • and φ Q cannot be larger than 11.6 • . D. Pseudoscalar glueball production in J/ψ radiative decay In Fig. 5 the glueball mass in term of the favored range for φ G is plotted. The band indicates the boundary of G g with (−0.054 ± 0.008) GeV 3 . It should be noted that, if M G = 2.56 GeV from the LQCD calculation [26] is taken, φ G would be fixed as 7 • , and the corresponding U matrix is given as, The production rate of P in J/ψ → γP scales as (G P /M 2 ) 2 [16], whereM is a typical energy scale for 0|α s (M )GG/(4π)|(qq) 0 −+ . This energy scale also determines the strong coupling α s (M ). It is natural to expect that this energy scale is the same for the light pseudoscalar meson productions, i.e. η and η ′ , in the light flavor sector. However, it should be different for η c due to the much shorter range for the color force between c andc and larger momentum transfers to the gluons in the cc → gg transition. We can examine the fitted values in Eq. (22) for η and η ′ , and then estimate the production rate for the pseudoscalar glueball.
The branching ratio fraction for the production of two pseudoscalar mesons P 1 and P 2 in the J/ψ radiative decays can be expressed as where q 1 and q 2 are the three-vector momentum of the pseudoscalar meson P 1 and P 2 in the rest frame of J/ψ, respectively, whileM 1 andM 2 are the energy scales for the strong quark-gluon couplings for P 1 and P 2 , respectively, in the qq → gg transition. As mentioned earlier, it is a reasonable approximation to adopt the sameM value for η and η ′ . With the data for BR(J/ψ → γη ′ ) and BR(J/ψ → γη) from experiment [77] it allows us to extract G η ′ /G η = 2.39 +0.08 −0.15 , where the central values of the data give the ratio, and the boundaries are given by the upper and lower limit of the data uncertainties [77]. The theoretical value from Eq. (22) gives G η ′ /G η = 3.19 which is close to the data constraint taking into account the uncertainties arising from the parameters. The branching ratio fraction can be related to the η-η ′ mixing either with or without the glueball mixing which indicates the small glueball component in the η and η ′ wavefunction as found in the literature [18,62,63,67,68].
For the glueball production in J/ψ radiative decays, one can calibrate its production to the rate for J/ψ → γη via whereM g andM η are the energy scales for the glueball and η. Note that in case there is a large mass difference between the physical glueball mass and η it is unnecessary forM g =M η . Early studies of the scale relation can be found in Ref. [83]. As an approximation we assumeM g =M η to extract BR(J/ψ → γG) with a mass of M G = 2.1 GeV, i.e. BR(J/ψ → γG) ≃ 3.8 × 10 −3 . Taking into account thatM g is supposed to be larger thanM η , this rate sets up an upper limit to the branching ratio for the production of an M G = 2.1 GeV pseudoscalar glueball in the J/ψ radiative decays. This result is consistent with the analysis of Ref. [78] which pointed out the difficulty of reconciling the LQCD result with the experimental hint for the possible existence of the additional η(1405). One notices that in this mass region BESII reported pseudoscalar states X(2120) and X(2370) in the invariant mass spectrum of η ′ ππ in J/ψ → γX → γη ′ ππ [79,80], which was confirmed by BESIII later with high statistics [81]. The PDG also list η(2225) as an established state in J/ψ → γKKπ with BR(J/ψ → γη(2225) = (3.14 +0. 50 −0.19 ) × 10 −4 [77]. Whether these states are radial excitations of η and η ′ [13] or whether one of these states is the pseudoscalar glueball candidate should be further investigated in both experiment and theory.

IV. SUMMARY
In this work, we revisit the mechanism proposed and studied in Refs. [16,60] for the pseudoscalar meson and glueball mixings. On the one hand, we confirm many results from Refs. [16,60] on the correlations among the introduced parameters. On the other hand, we scrutinize the dynamical constraints on the glueball mass and clarify that the physically favored parameter space would lead to much higher glueball mass than that obtained before. In particular, we show that the approximation of neglecting both m 2 qg and m 2 sg in the extraction of the glueball mass was inappropriate. Although m 2 qg is indeed a negligible quantity in comparison with √ 2G g /f q , the value of m 2 sq is actually comparable with G g /f s and cannot be neglected. After properly treat these parameters and identify G g and φ G as the parameters that play a dominant role in the determination of the mixing pattern, we find that the glueball mass M G cannot be lower than 1.9 GeV which is much higher than 1.4 GeV determined by the approximation in Refs. [16,60]. We find that the mixing angle φ G ∈ (7, 12) • for the glueball and light flavor states is favored in our model. It allows the estimate of an upper limit branching ratio of the production of pseudoscalar glueball in J/ψ radiative decay.
This result is encouraging in such a sense that it resolves not only the apparent paradox between the LQCD results and some old experimental data for the pseudoscalar glueball mass, but also explains the single peak structure around 1.4∼1.5 GeV observed by the recent high-statistics measurements at BESIII in exclusive decay channels [34,80,82], although in some channels the peak positions are slightly shifted. The peak position shift is due to interferences from the triangle singularity mechanism [29,30,43] instead of by two-pole structures. Namely, we can conclude that there is no need for two light pseudoscalars η(1405) and η(1475) to be present as two individual states and with the η(1405) as the pseudoscalar glueball candidate. Consequently, as pointed out in Refs. [29,30,43], in order to search for the pseudoscalar glueball candidate one should look at the higher mass region at least above 1.8 GeV where some of the recently observed pseudoscalar states by BESIII [82] should be carefully examined.