Internal particle width effects on the the triangle singularity mechanism in the study of the $\eta(1405)$ and $\eta(1475)$ puzzle

In this article, the analyticity of triangle loop integral with complex masses of internal particles is discussed in a new perspective, base on which we obtain the explicit width dependence of the absorptive part of the triangle amplitude. We reanalyze the decay pattern of $\eta(1405/1475)$ with the width effects included in the triangle singularity (TS) mechanism. Based on the present experimental information, we provide a self-consistent description of the $K\bar{K}\pi$, $\eta\pi\pi$, and $3\pi$ decay channels for $\eta(1405/1475)$. Our results confirm the claim that the TS mechanism plays a decisive role in the understanding of the $\eta(1405)$ and $\eta(1475)$ puzzle. Namely, the observed differences of $\eta$ resonances within the mass region of $1.40\sim 1.48$ GeV are originated from the same state. For the isospin violated process $J/\psi\to\gamma\eta(1405/1475)\to f_0(980)\pi\to 3\pi$, we identify an additional contribution to the $a_0(980)-f_0(980)$ mixing via the TS mechanism.


I. INTRODUCTION
It has been a long-standing question on the existence of glueball in hadron spectroscopy. This exotic object as the bound state of gluons predicted by QCD has been a crucial piece of information for our understanding of strong interaction theory in the non-perturbative regime. The corresponding theoretical study and experimental search for glueball states have been the topical subjects in hadron physics. However, although tremendous efforts have been made, the indisputable evidence for their existence is still lacking. In the glueball spectrum the low-lying states include scalar (J P C = 0 ++ ), tensor (2 ++ ) and pseudoscalar (0 −+ ). So far, the lattice QCD (LQCD) simulations [1][2][3][4][5][6] suggest that their typical masses are about 1.5 ∼ 1.7, 2.2 ∼ 2.3, and 2.4 ∼ 2.6 GeV, respectively. The mass hierarchy M 0 ++ < M 2 ++ < M 0 −+ has been a stable feature from LQCD. While there have been topical reviews on the glueball spectrum in the literature for the scalar and tensor, our focus in this work is on the pseudoscalar glueball. We will discuss the long-standing controversial issues involved in the identification of the pseudoscalar glueball candidate, and stress that a self-consistent picture can only be obtained with a special kinematic effect, i.e. the triangle singularity (TS) or Landau singularity mechanism [7][8][9][10], implemented.
In the literature the most promising candidate for the pseudoscalar glueball has been assigned to the η(1405) since it was introduced as an additional state to the nearby η(1475) and η(1295) around early 1990s [11][12][13][14][15][16][17][18]. These three states of isospin 0 and similar masses cannot fit the pattern arising from the SU(3) flavor symmetry of quark model in the light quark sector. A broadly accepted classification has been that the η(1295) and η(1475) belong to the isospin-0 radial excitation states in the SU(3) flavor multiplet due to the mixing between the flavor singlet and octet. The η(1405) as an out-numbered state was then proposed to be the ground state pseudoscalar glueball candidate. Such an assignment was based on phenomenological studies which predicted the mass of the pseudoscalar glueball around 1.4 GeV [19]. This proposal seemed to be able to accommodate the experimental observations with the theoretical prediction, and had attracted a lot of efforts to further explore the structure and production mechanism of the η(1405) as the pseudoscalar glueball candidate [20][21][22][23][24][25][26].
However, because of model-dependence it has been very difficult to make progress on establishing unambiguously the glueball nature of η(1405).
It is a challenge to bring down the pure gauge glueball mass from ∼ 2.5 GeV to ∼ 1.4 GeV. The very relevant issue is that whatever the mechanism could be it requires an abnormally strong coupling between the light quark states and pure gauge glueball. It has been quoted broadly that the QCD sum rules for pseudoscalar glueball led to relatively low masses. However, it should be noted that the pseudoscalar glueball sum rules are very sensitive to assumptions made in the calculations. As noted explicitly in Ref. [34] the pseudoscalar glueball mass in QCD sum rules has large uncertainties and is very sensitive to the gluon condensation. The question about the pseudoscalar glueball mass is inevitably correlated with the η-η ′ mixing because of the axial anomaly [35][36][37]. The deviation of the η-η ′ mixing angle from the ideal one between the flavor octet and singlet indicates the crucial role played by the anomaly. Meanwhile, the η-η ′ mixing angle does not determine the η ′ mass as a usual mixing scheme would suggest. On the contrary, its dependence of the topological charge density 0|G a µνG aµν |η ′ has to be taken into account. In Ref. [28] a dynamical approach for the η-η ′ -glueball mixing was explored by implementing the mixing into the equations of motion for the anomalous Ward identity and a low mass about 1.4 GeV for the physical pseudoscalar glueball. This approach was extended to accommodate the η c in Ref. [38] and a similar result was extracted. However, an analysis of Ref. [39] based on the same dynamics yields a lower bound of about 2 GeV for the pseudoscalar glueball mass. In Ref. [40] a revisit of the mixing scheme of Refs. [28,38] was carried out, and the numerical results of Refs. [28,38] were confirmed except that the approximation for extracting the pseudoscalar glueball mass appeared to be problematic [64]. After curing this problem, it shows that the physical glueball mass will favor to be higher than 2 GeV which is remarkably consistent with the conclusion of Ref. [39] and matches the LQCD simulations. It is interesting to note that the analysis of Ref. [41] also suggests a pseudoscalar glueball mass above 2 GeV, although the physical state should be lighter than the quenched state from the pure Yang-Mills gauge theory [41]. The mass difference between the quenched pure gauge state and the QCD state is of the order of 1/N c , which means that a low mass state around 1.4 GeV is unfavored. In fact, a lot of puzzling questions arised from not only conflicts between the experimental observation and LQCD results, but also between the early phenomenological studies and the first principle LQCD simulations.
The change of situation was triggered by the high-statistics experimental data from BESIII. There have been a number of the J/ψ exclusive decay channels measured with high precision, where contributions from J P C = 0 −+ can be clearly identified [42][43][44]. It shows that in the vicinity of 1.4 GeV there is only one Breit-Wigner peak structure in the invariant mass spectrum for η(1405/1475). Similar feature can be found in the radiative decays of ψ ′ . In the hadronic production channels, such as e + e − → ωηππ, φηππ, ωKKπ, etc., there is also only one Breit-Wigner peak present in the invariant spectrum. The interesting observation is that the peak positions somehow are slightly different in exclusive channels. Further insights into this puzzling problem was gained from the measurement of isospin breaking effects into J/ψ → γη(1405/1475) with η(1405/1475) → 3π at BESIII [42], where the isospin breaking effects are found to be unexpectedly large, i.e. Γ(η(1405) → f 0 π → π + π − π 0 )/Γ(η(1405) → a 0 π → ηπ 0 π 0 ) ∼ 17.9%. This value is nearly one order of magnitude larger than that produced by the a 0 (980) − f 0 (980) mixing. It was then discovered by the authors of Ref. [45] that the significantly enhanced isospin breaking effects are caused by the so-called "triangle singularity" (TS) mechanism. It was demonstrated in Ref. [45] and later a detailed analysis [46] that at the mass of η(1405/1475) the non-vanishing coupling of the initial η(1405/1475) to the intermediate K * K + c.c. and then their rescatterings into πf 0 (980) by the exchange of a Kaon or anti-Kaon allow a perfect satisfaction of the TS condition. While the detailed analytical properties of the triangle diagrams will be discussed later, a simple way to picture the TS mechanism is that it corresponds to such a kinematic condition that all the internal states of the triangle loop can approach their on-shell condition simultaneously. As a consequence of such a leading singularity within the loop function, it will provide significant interferences in exclusive decays of η(1405/1475) and produce the shift of peak positions of a single state in different channels and unexpectedly large isospin breaking effects in its decays into 3π [45,46]. The dominance of the TS mechanism in η(1405/1475) → 3π was later confirmed by Ref. [48] in a chiral unitary approach.
Although the TS mechanism seems to be promising for understanding so-far all the existing puzzles about the η(1405) and η(1475) signals, later according to Ref. [47], the non-zero width of K * in the triangle loop integral can lead to significant suppressions of the decay rate. Therefore, the dominance of triangle diagrams in the isospin violated channel may become questionable. Implications of such a possibility suggests that the width effects due to the internal states should not be neglected and may lead to a significant impact on the role played by the TS mechanism. In order to clarify this, we carry out a coherent and quantitative investigation of the decays of η(1405/1475) into KKπ, ηππ and 3π including the width effects in the TS mechanism.
As follows, in Sect. II we first introduce the TS conditions and discuss in detail the analytical properties of the triangle loop amplitude when non-zero K * width is considered. Then we will explore the decay patterns for η(1405/1475) → KKπ, ηππ and 3π, and clarify the role played by the TS mechanism in Sect. III. In particular, we will show that an additional transition process which enhances the a 0 (980) − f 0 (980) mixing via the TS mechanism should contribute to η(1405/1475) → 3π. The calculation results and discussions will be given in Sect. IV, and a conclusion will be given in Sect. V.

II. ANALYTIC PROPERTIES OF THE TRIANGLE LOOP INTEGRAL WITH NON-ZERO WIDTHS
To understand the impact of unstable internal K * meson on the triangle loop integral of η(1405/1475) (denoted by η ′′ in the following), we consider a typical triangle amplitude I shown in Fig. 1, where the m 1 , m 2 and m 3 are the masses for K * , K andK, respectively. Since K * meson is unstable, a finite width Γ has been introduced in its propagator. Due to the P -wave vertex and the polarization of K * meson, the amplitude I is actually a tensor integral. However, this tensor integral can be reduced to a sum of scalar 3-point integral and some 2-point integrals [47]. By studying the analytical property of the 3-point and 2-point integrals, we can learn the property of the physical amplitude I. The TS condition applies to the physical amplitude I where all the internal particles approach their on-shell condition simultaneously. In such a sense the reduction of I into 3-point and 2-point loop integrals means that the manifestation of the physical TS contributions is given by the sum of the reduced loop integrals although some contributions are from the 2-point loops. Such a clarification is essential for the reason that we actually deal with the physical process instead of single loop integrals which are only part of the dynamics of the physical process. The kinematic condition that all the internal particles are on-shell determines the kinematics for all the reduced loops which have to be considered simultaneously. To be more specific, within the TS kinematics the two-body on shell condition in the 2-point integrals has been contained in the 3-body on-shell condition. With the above clarification, the TS contribution in this work is referred to the overall contributions from the physical integral I instead of a reduced 3-point integral, and the influence of the finite width effects are also referred to its impact on the overall loop function. Some more features about the TS loops should be pointed out before we proceed to the detailed analysis: • The presence of the TS kinematics implies that the main contributions of the triangle loops come from the kinematic region near the on-shell condition for the internal particles. For physical processes where the internal motions of the internal particles can be treated non-relativistically the scalar triangle loop can be directly integrated out, and the leading logarithmic singularity can be explicitly extracted. In particular, for nonrelativistic heavy meson loop transitions where the TS mechanism is present, the loop amplitudes can be analyzed in the non-relativistic effective field theory (NREFT) framework and a power-counting scheme can be established.
• For loop transitions involving only light hadrons sometimes the non-relativistic approximation can hardly be justified. In such a case analysis of the triangle loop in the Mandelstam representation should be more appropriate. For most cases of the physical loop integrals, an empirical form factor has to be included to cut off ultraviolet divergence when the internal particles go off shell, which will inevitably introduce some modeldependence, although in general such an uncertainty can be under control. Even for convergent physical loops it is often checked that unphysical contributions from relatively large momentum transfers are reliably estimated and then removed [49]. For the physical triangle loop of Eq. (1), it converges with the choice of the vector propagator (−g µν + q µ q ν q 2 )/(q 2 − m 2 1 + im 1 Γ) for the K * (K * ). But in order to examine the sensitivities of the loop integrals to unphysical contributions from the ultraviolet region, we include a form factor F (q 2 ) and compare the numerical results for different cut-off parameters. The detailed discussions will be given later in Section III and IV.

A. Analytical expression
We first consider a typical scalar loop integral where the notation is the same as that in Fig. 1. Taking the Feynman parameterization, it can be expressed as where D is a homogeneous polynomial of x i , i.e.
In the case where the internal masses are real, the leading singularity of M is determined by the Landau equation This condition means that the extremum of D touches 0. Physically, it means that the internal particles are stable and simultaneously become on-shell. The Landau equation gives Y ij x j = 0. If this holds for 0 < x i < 1, then there exist solutions for the Landau equation within the physical region which satisfies det Y = 0. This leads to the kinematic bounds for the triangle singularity. In general, for the fixed external four-momentum squares s 2 and s 3 there are two solutions for s 1 which satisfy the Landau equation, i.e.
However, the conditions that s 1 = s + 1 with the fixed s 2 and s 3 or s 3 = s + 3 with the fixed s 2 and s 1 do not cause divergence because the corresponding solution x i is out of the region (0, 1). Only when s 1 (s 3 ) meets s − 1 (s − 3 ), the triangle singularity occurs within the physical region and can possibly produce detectable effects in experimental observables.
With the s 2 fixed (since it is an external particle), for any given (m 2 + m 3 ) 2 < s 3 < s 3c , there is a s − 1 lying within ((m 1 + m 3 ) 2 , s 1c ), where s 1c and s 3c are labeled as critical values for s 1 and s 3 , beyond which the triangle singularity no longer exists. Taking the same notation as Ref. [50], these critical values are given by A caveat arising from the above discussion is that the internal particles are stable ones, i.e. they do not have a width in the propagator. In reality the nonvanishing coupling for K * (m 1 ) to pion (s 1/2 2 ) and kaon (m 2 ) demonds that the propagator for K * must contain an imaginary part. Therefore, a detailed investigation of the width effects in the TS mechanism is necessary and useful for a better understanding of the underlying dynamics.
To accommodate the width effects in the triangle loops, we consider complex masses for the internal particles. With the help of Spence function: which is also called dilogarithm function Li 2 , as a special case for polylogarithm function Li n when n = 2, an analytic expression of the transition matrix M has been worked out by G.'t Hooft and M. Veltman [51]: where functions u i (y) have the following expressions: and a, b, c, d, e, f , and α are kinematic variables: This analytic expression is valid for both real and complex internal masses. In Eq. (13) N i (y) are functions of the integration variable y: and y = R(y 0 , y 1 ) + R(y 0 , y 2 ) where R(y 0 , y 1 ) with The η function arises from ln (z 1 z 2 ) = ln z 1 + ln z 2 + η(z 1 , z 2 ), with the argument in ln limited in (−π, π).
The following features with the TS kinematics will help understand better the analytical properties of the triangle loop amplitude: • In the vicinity of the TS kinematics the main contributions of the transition amplitude are given by the absorptive part. In particular, when the internal masses are real, the logarithmic divergence of the TS manifests itself in the absorptive part of scalar integral. So we will focus on the absorptive part in the analysis.
• For the case that the internal masses are all real, the absorptive part can be derived analytically according to the Cutkosky rule. This allows us to examine the width effects on the triangle loop amplitude by comparing them with the Cutkosky rule result. In particular, the width dependence can be highlighted in the absorptive part of the amplitude.
• For the physical case, namely the isospin-violating transition η(1405/1475) → f 0 (980)π, it should be noted that the dispersive part becomes negligible due to the cancellation between the charged and neutral triangle loop amplitudes. This actually leads to a rather model-independent behavior of the TS contributions to the isospin violations in η(1405/1475) → f 0 (980)π [65]. As stressed at the beginning, the only thing left behind is the width effect that should be quantified by explicit and self-consistent calculations.
In the following sections we focus on the derivation of absorptive part of the scalar integral under the influence of the finite width of the intermediate states. We will provide detail analysis of ℑm 1 < 0 and m 2 , m 3 ∈ R, for the reason that the physical widths of the K (K) by weak decays are much smaller than that of the K * meson. However, it is checked in the end that the same analytic expression is still valid when any of the internal state has a complex mass.

B. The motion of singularities
For convenience, we express the amplitude of Eq. (14) in a concise form with Our task in this subsection is to extract the imaginary part of the loop amplitude and investigate the movement of the singular kinematics manifested by the location of y (i) k in the complex plane. The general condition for the TS requires that the following kinematic constraints are satisfied, i.e. s 1 > (m 1 + m 3 ) 2 , s 3 > (m 2 + m 3 ) 2 , and √ s 2 < m 1 − m 2 . Moreover, the maximum allowed value for s 1 or s 3 is generally very close to the normal threshold. It allow us to make a substitution of s 3 by (m 2 + m 3 ) 2 in some steps as a reasonable approximation. Therefore, in the following discussion of the finite width effects on the imaginary part of the amplitude, we can apply this approximation to simplify the deduction without loss of accuracy. The numerical result of this approximation compared with the exact one will be discussed in the end of this section.
To proceed, we will start with W  1 is the larger one of the two roots. In the case of finite Γ, y (i) k may deviate slightly from real axis which will be our focus in this work.
According to the definition of S (i) (Eq. 24), the divergence of S (i) occurs when y In order to discuss how y (i) k moves when the K * meson has a small but finite width, we need to first write down the explicit expression of y (i) k as the roots of function u i (y) and then expand the y (i) k with respect to m 2 1 (the mass square of the K * ) to the first order of Γ.
Since u 1 (y) has nothing to do with the complex mass m 1 , y is not affected by Γ. Thus, their locations on the real axis are The only difference between y 1 as the larger root of u 2 (y) and u 3 (y), respectively, are Hence the variations of the locations of y 1 and y 1 in association with the presence of the finite width Γ for the K * can be examined by It should be noted that within the region s 2 < (m 1 − m 2 ) 2 , the value of E 2 + 1 is always negative, which means that y (2) 2 actually moves upward away from the real axis and y (2) 1 moves downward away from the real axis with the increasing Γ. These expansions give the sign for the imaginary parts, according to which we can conclude that both y For m 1 > m 2 > 0 and s 2 > 0, it can be easily verified that y j are quite clear, and we are to take a look at y but the derivatives, which can indicate their motion on the complex plane, are much simpler: where With the explicit expressions for y Each of such trajectories is illustrated by a set of solid dots (with dotted arrows in (a)) and empty circles (with solid arrows). The dashed vertical lines indicate the fulfilled TS conditions with the regime between the solid dots (in (a) the solid dots will move towards the dashed line and match the value of ℜy , these two points overlap horizontally which is marked by the thick (red) dashed vertical line.
• The thin (green) arrow lines (solid and dotted) indicate the situation with s 1 = s 1c . Then, y   only occurs when s 3 = s 3c , which is marked by the (red) cross. If s 1 is set to be a larger value, the trajectory of y are plotted in Fig. 2 (c). The condition that s 1 = (m 1 + m 3 ) 2 is also marked by a set of dots (solid and empty) with a thick (red) arrow line. The vertical dashed lines with crosses mark the TS condition when s 3 varies. In this case, as s 3 increases from the normal threshold, the empty (red) circle (y is coincident with the starting point of y (3) 0 , which is marked by the right (green) cross. If s 1 is smaller than the normal threshold, the y will rapidly run away from the real axis, so that there is no longer enhancement. Therefore, the valid region for the occurrence of the TS (s 1 = s − 1 , s 3 = s − 3 ) corresponds to the line segment from the right (green) cross to the left (red) one.
A crucial fact is that when s 1 = s − 1 , the overlaps of y occur simultaneously, which means that S (1) , S (2) and S (3) simultaneously become enhanced (or maximized). This exactly corresponds to the TS condition for the triangle loop.
There may be some other horizontal overlaps, which however, are not responsible for the TS. They happen either when s 1 = s + 1 or when s 3 = s + 3 , and these effects eventually cancel out when sum over all the S (i) in the final result of ℑM . We do not show these in Fig. 2 since they are not the focus of this analysis. However, the motions of y based on the Spence function, the occurrence of TS is recognized as the simultaneous (horizontally) overlaps of these singular points on the complex plane.
To extract the imaginary part from the Spence functions, we adopt the following formula where a is a real number and b is a positive real number. If b is much smaller than a, the second term is to the order of b and can be neglected. In this case, we obtain the following approximation ℑSp(a ± ib) = ±πΘ(a − 1) ln a.
We first discuss the sign of the real and imaginary part of z 12 | to the imaginary part of M in different kinematic regions. Energy

F. Reduction of ℑM
The imaginary part of M thus becomes ℑM = 1 16π 2 λ[s 1 , s 2 , s 3 ] 1/2 (ℑW tot + ℑT tot + ℑσ tot ). (69) In principle, Eq. (69) can be used to directly calculate the imaginary part of scalar loop integral in the case where the internal mass has a small negative imaginary part. In the case where m 2 or m 3 has imaginary part, one can reconsider the motions of poles in complex plane, redetermine the proper sign of each parameter to obtain new result, which will be identical to Eq. (69). That is to say, Eq. (69) can be used to do calculations when any of the internal mass has a small imaginary part. In the case where Γ = 0, it will be shown later that this expression is numerically identical to the ℑM calculated directly from the Cutkosky rule.
Although the logarithmic dependence is evident from Eq. (69), its form is still very complicated. For convenience, we denote the imaginary part derived by the Cutkosky rule by f c . Before doing any reduction, it is very interesting to notice that by replacing the internal mass m 2 1 (here m 1 should be a real mass) in f c with its complex form m 2 1 − im 1 Γ, we obtain the same result as what is calculated by Eq. (69). In Fig. 3, we show the imaginary part of the scalar loop corresponding to η ′′ → K * K + c.c. → KKπ → f 0 (980)π with and without the K * width. Different approaches give the same value for Γ = 0. Visible difference between these methods appear when Γ = 50 MeV. But the thick (red) solid line which is given by the replacement of m 2 1 with m 2 1 − im 1 Γ, and the thin (green) solid line which is given by Eq. (69), always match even when Γ is up to 200 MeV.
The disagreement between Eq. (14) (yellow dot-dashed line) and Eq. (69) (green dashed line) is due to several aspects. Firstly, when the width grows, the step function is no longer a good approximation, and a continuous function should be used to describe the parameters. Secondly, some substitutions using s 3 = (m 2 + m 3 ) 2 should be carefully dealt with. Thirdly, it is not a good approximation to omit the second term in Eq. (52). Instead, Eq (52) is suggested to be replaced by where D 2 is the Bloch-Wigner function [52], which a real analytic function in the cut plane C \ {0, 1}. Despite of this disagreement, the matching of thick (red) solid and (green) dashed lines encourages us to believe that the reduced result of Eq. (69) actually takes exactly the same form as f c with complex internal masses. This equivalence is non-trivial, for the fact that if f c is regarded as a function of internal masses, it is real so that it must be nonanalytic everywhere in the complex plane of internal masses, which makes it hard to be continued into the complex plane. A possible reason is that Eq. (69) and f c are the imaginary parts of two analytic functions of internal mass, respectively, then the uniqueness of analytic functions ensures their equivalence. However, this equivalence shows that the f c with complex internal masses can be regarded as the actual value of ℑM , given the width of internal particle is not too big to break the assumptions made in the previous discussion.
This conjecture has profound consequences. Since the analytic form of f c is much simpler than Eq. (69), if the equivalence holds, one can use f c with complex masses to isolate the Γ dependence, which may provide the answer to how sensitive the absorptive part of M is to Γ. Now we can concentrate on the process η ′′ → K * K + c.c. → KKπ → f 0 (980)π. If the equivalence is true, under the circumstance where internal masses all have a small imaginary part, the absorptive part of the corresponding scalar integral is where with m (ℜm = m 1 ) and m k being the masses of K * and K mesons, respectively. In the case where internal masses are real, as s 1 approaches the abnormal threshold s − 1 , both a 1 and a 2 approaches 1, which leads to a logarithmic  divergence. When internal masses are taken to be complex, as s 1 → s − 1 , a → 1 + iO(Γ), so that the divergence disappear and the ℑM is expected to be like ln β 2 Γ 2 . In the case where m 2 has an imaginary part −im 1 Γ, we can expand a 1 and a 2 to the first order of Γ, and take the limit s 1 → s − 1 to see how the maximum value of ℑM becomes divergent as Γ → 0. The expansion is not difficult so we just write down the result: with  where the kinematic variables are defined as follows: Since some of the Γ terms are contained in terms like λ[s 1 , m 2 , m 2 k ] with the complex masses, the closer to normal threshold the s 1 is, the smaller the convergence radius would be, which means this approximation becomes worse when s 1 is closer to (m + m k ) 2 . To see whether the amplitude is sensitive to Γ as it varies a small value in the vicinity around its physical value, we can calculate the logarithmic derivative of ℑM , i.e.
and the results for δΓ = Γ × 10% MeV with Γ = 10 MeV and 50 MeV, respectively, are illustrated by Fig. 4. From  Fig. 4, we can conclude that, if Γ varies 10%, the typical variation of max ℑM is around 4% for Γ = 10 MeV and 10% for Γ = 50 MeV, which are rather small effects. Naively, one may consider that the Γ dependence is described by Γ/ℜm. On the contrary, the quantity β in Eq. (74) should be identified as the energy scale to describe the strength of Γ dependence. Figure 5 shows that within the region (ℜm + m k ) 2 < s 1 < s 1c , the value of β is much smaller than m (∼ 0.9 GeV). Figure 6 shows that, when Γ = 10 MeV, the approximation holds well enough. When Γ = 50 MeV, the approximation breaks down for small values of s 1 . It shows that although Γ = 50 MeV is almost 1/18 of ℜm = m 1 , compared  to the scale β, it is actually a large value and an expansion to high order of Γ is required for small value of s 1 . This feature suggests that for relatively narrow widths the treatment of Ref. [53] can also provide a reliable evaluation of the TS contributions. Figure 6 also shows that when s 1 goes closer to s 1c , the amplitude becomes larger. Hence, it can be expected that the larger the allowed kinematic region for the TS, the stronger the enhancement would be. This is consistent with the analysis of Ref. [50].

A. Experimental aspects
In this section we will revisit the radiative decay J/ψ → γη(1405/1475) with the η(1405/1475) decays into KKπ, ηππ and 3π with the TS mechanism considered. The focus is to quantify the width effects based on the formulation developed in the previous section and present a coherent analysis of these exclusive channels on the same basis. For abbreviation, we denote η(1405/1475) by η ′′ as follows.
We are going to adopt the J/ψ radiative decay data from the BESIII Collaboration in this analysis. As pointed out in the Introduction, the present PDG [54] lists η(1405) and η(1475) as different states appearing in different channels. However, the high-statistic data analyses at BESIII so far have not observed two separated states in any exclusive channel. In order to keep self-consistent we only consider the J/ψ radiative decays and assume that there is only one state η ′′ contributing. We then examine all observables with the implementation of the TS mechanism as a self-consistent check of our one-state scenario. This strategy is the same as that in Ref. [46] except that we will quantify the width effects in the loop function here and point out an additional contribution to the isospin violation channel.
We adopt the PDG averaged combined branching ratios BR(J/ψ → γη ′′ → γKKπ) = (2.8 ± 0.6) × 10 −3 and BR(J/ψ → γη ′′ → γa 0 π → γηπ 0 π 0 ) = (8.40 ± 1.75) × 10 −5 [54] to constrain the parameters in this analysis. Note that the main weights for extracting these averaged branching ratios are from the BESIII measurements. One obvious feature is that in both processes there is only one state of (I = 0, J P C = 0 −+ ) has been observed around 1.4∼1. 48 GeV. Assuming that these two signals are from the same state, the much larger branching ratio for the γKKπ channel than γηππ should be one of the important observables for our understanding of the structure of η ′′ .
As studied in Refs. [45,46], the a 0 (980) − f 0 (980) mixing mechanism can only account for a few percent of isospin breaking effects. Thus, other mechanism is required to compensate such a significant deficit.
It is also interesting to note that the a 0 (980) − f 0 (980) mixing measured in J/ψ → φηπ and χ c1 → 3π at BESIII is at the order of 1% [55], which is consistent with the theoretical prediction of Refs. [56,57] and can be regarded as a direct evidence for additional mechanism that strongly enhances the isospin breakings in η ′′ → f 0 (980)π → 3π. Note that the TS condition is not fulfilled in J/ψ → φηπ and χ c1 → 3π. Taking into account this fact, the role played by the TS in η ′′ decays is rather natural and easy to test [45,46]. B. The amplitudes for η ′′ → f0π → π + π − π 0 Taking into account the TS mechanism, the isospin violating decay of η ′′ → 3π can be accounted for by Fig. 7, where Fig. 7 (a) and (b) describe the TS mechanism which can contribute via two processes, and Fig. 7 (c) describes the a 0 (980) − f 0 (980) mixing mechanism via the leading tree-level diagram. Namely, in Fig. 7 (a) the isospin violation arises from the mass difference between the neutral and charged triangle loops, and in Fig. 7 (b) the triangle loop provides an enhanced production mechanism for a 0 (980) which then mixes with f 0 (980) to contribute to the isospin breaking effects. It should be emphasized that although these two processes both are categorized as the TS mechanism, they represent different dynamics contributing to the significant isospin breakings.
It should be addressed that a coherent study including Fig. 7 (b) and (c) has not been carried out. In Refs. [45,46] process of Fig. 7 (a) is considered with the coupling vertices determined by experimental data. It can be regarded as a leading approximation for the vertex renormalization. In addition, since the dispersive parts from the charged and neutral loop amplitudes cancel each other, the dominance of the absorptive parts largely reduces the modeldependence in the estimate of the loop contributions. But the width effects were not estimated therein. In Ref. [47] the authors investigated the width effects arising from Fig. 7 (a), while the process of Fig. 7 (c) is estimated by the data for η ′′ → a 0 (980)π which is also a leading approximation for the coupling renormalization for the η ′′ a 0 (980)π vertex. However, the experimental data for η ′′ → a 0 (980)π still contain large uncertainties either in η ′′ → KKπ or ηππ. In Ref. [48], similar processes as Fig. 7 (a) and (c) are considered and the authors confirm the irreplaceable contributions from the TS mechanism. However, the width effects are not investigated. In order to obtain reliable results for the role played by the TS mechanism, we are going to coherently include Fig. 7 (a) (b) and (c), and the width effects within the triangle loops. Meanwhile, a combined analysis of η ′′ → KKπ and ηππ in association with the isospin violating decay of η ′′ → 3π is necessary.
FIG. 7: Isospin violating mechanisms for η ′′ → π + π − π 0 . (a) describes the direct isospin breaking via the TS mechanism, (b) describes the TS contribution to the production of the a0(980) which then mixes with the f0(980) via the KK loops, and (c) represents the tree-level contributions to the intermediate a0(980) production, which will then mixing with the f0(980) via the KK loops.
To proceed, the effective Lagrangian for vector-pseudoscalar-pseudoscalar (V P P ) vertex is defined as where V and P stand for the fields for the flavor SU(3) multiplets, respectively, and they have the following expressions: and The overall coupling strength g V P P in Eq. (78) for different channels can be determined by experimental data. For instance, with the partial decay width of K * → Kπ [54], the couplingg K * 0 K 0 π 0 ≡ g V P P / √ 2 = 3.20 can be determined, and then we extract g V P P = 4.53. The relative signs for g η ′′ K * K in different exclusive coupling channels are given by the charge conjugation symmetry, i.e.
The triangle amplitude for Fig. 7 (a) can then be expressed as where the function I 1,2 (s 1 , s 3 ) has the same expression as Eq. (1) and the subscriptions 1 and 2 denote the neutral and charged triangle loop integrals, respectively. Note that I 1,2 (s 1 , s 3 ) are complex functions of s 1 and s 3 . When the kinematic reflection is considered, s 3 should be replaced by the invariant mass squared of the corresponding 2-body system. Equation (1) converges with the choice of the gauge for the vector meson. However, since we are dealing with a hadronic loop which involves non-fundamental particles it is possible that unphysical contributions can be present in the dispersive part of the amplitude and should be subtracted. To examine such a possibility we include a form factor F (q 2 ) in the numerical calculations. By changing the different values for the cut-off parameter, we are able to examine uncertainties arising from unphysical contributions from the ultraviolet kinematic region: where p R with R = K * , K,K are the four-momenta of the corresponding intermediate particles and Λ R ≡ m R + αΛ QCD (with Λ QCD ≃ 200 ∼ 300 MeV) is the cut-off energy. This choice of form factor is not perfect but convenient and sufficiently efficient. It implies that the vertex couplings will be suppressed as long as the propagators become offshell towards the space-like region. One notices that the monopole form factor does not guarantee that the amplitude will be suppressed when the internal particles go off shell towards the space-like region. However, it can be easily checked that in the vicinity of the TS kinematics, in particular, when at least two particles are on shell, the third one will mostly stay in the space-like region. The kinematics when one state moves into the time-like region and the other two stay in the space-like region can still be efficiently cut-off when the kinematics deviate away from the TS condition. Such a situation can be recognized by the following decomposition: To be more specific, since Λ R is sufficiently larger than m R the artificial pole contribution from p 2 R = Λ 2 R is far away from the vicinity of p 2 R = m 2 R . Therefore, the artificial pole is not going to contribute to the physical value of the loop amplitude. Meanwhile, when p 2 R is sufficiently far away from the physical pole of m 2 R , the cancellation in Eq. (83) is obvious and traceable.
For the processes of Fig. 7 (b) and (c) both processes contribute to the production of a 0 (980) and then a 0 (980) − f 0 (980) mixing leads to the isospin violation. Note that in Fig. 7 (b) the neutral and charged triangle loop amplitudes will constructively add to each other. The contributions from the dispersive amplitude within the triangle loops will suffer from the uncertainties arising from the regularization of the loop divergence. This is very different from Fig. 7 (a) where the dispersive parts will cancel and the model-dependence is strongly suppressed. We will come back to this point in the numerical calculations in the next section.
For the a 0 (980) − f 0 (980) mixing we adopt the unitarized coupled-channel approach for their mixing via the intermediate KK meson loops and extract their propagators with the full loop corrections following the prescription of Ref. [58,59]. This will allow us to directly compare with the result of Ref. [47]. The general form of the propagators for f 0 (980) and a 0 (980) is written as the following: with where m is the physical mass defined at ℜΠ(m 2 ) = 0, k 2 ≡ s, and the subscription ab denotes the final state particles that a 0 (980) and f 0 (980) can decay into, i.e. ab = {ηπ, K 0K 0 , K + K − } for a 0 (980), and ab = {π 0 π 0 , π + π − , K 0K 0 , K + K − } for f 0 (980). Function Π(s) is the self-energy correction arising from the intermediate meson loops and iΠ ab is the two-point loop function. The g ab is the coupling of a 0 (980) → a + b or f 0 (980) → a + b, of which the adopted values are listed in Tab. VI. The interference among the amplitudes in Fig. 7 has an interesting implication of the partial decay width for η ′′ → ηππ. It is noticeable that the branching ratio of η ′′ → ηππ is much smaller than that of η ′′ → KKπ. So far, this question has not been seriously investigated because of the inconsistent treatment of the signals in these two channels, namely, they are treated as two irrelevant states η(1405) and η(1475). In our proposal their relative strength should contain dynamical information about the TS mechanism. One can see later that it can serve as a self-consistent examination of our solution.
Denoting the amplitudes of Fig. 7 (b) and (c) by M trimix and M treemix , respectively, we can write down their explicit expressions as follows: and where B 1 and B 2 are the loop functions for the neutral and charged kaon loops in the a 0 (980) − f 0 (980) mixing: In Eqs. (87) and (88) the relations g f0K 0K 0 = g f0K + K − , g a0K 0K 0 = −g a0K + K − and g η ′′ K * 0K 0 = g η ′′ K * + K − are implied for the vertex coupling constants.
C. The amplitudes for η ′′ → KKπ and ηππ A coherent study of the decays of η ′′ → KKπ and ηππ is required to quantify the contributions from the TS mechanism, and also to determine the parameters introduced in the formulation. For η ′′ → KKπ the decay can occur via the processes illustrated by Fig. 8, where (a) and (b) describe the tree-level transitions via the intermediates twobody decays η ′′ → K * K + c.c. and a 0 (980)π, respectively, and (c) describes the contribution from the TS mechanism. Note that the couplings of g η ′′ K * K and g η ′′ a0KK are bare couplings to be determined by experimental data with the dominant transition mechanisms included.
The amplitude for Fig. 8 (a) has the following form: The M tree1 is the sum of the amplitude of K * K channel and its charge conjugate: where the q and q ′ are the invariant mass of K * andK * , respectively. And for Fig. 8 (b) it reads .
The loop amplitude of Fig. 8 (c) reads where the triangle loop functions I 1,2 (s 1 , s 3 ) have been defined in the previous subsection. Similarly, the decay of η ′′ → ηππ is described by Fig. 9, where (a) is the tree-level transition via the intermediate a 0 (980)π and (b) describes the contribution from the TS mechanism. Again, we can see that the bare vertex couplings g η ′′ K * K and g η ′′ a0KK are the same as those appearing in Figs. 7 and 8. The corresponding amplitudes for Fig. 9 (a) and (b) are as follows: and

D. Coupling constants and relative phases
In Tab. VI all the vertex couplings are listed. Couplings g a0ηπ = 3.02 GeV and g a0K + K − = −g a0K 0K 0 = 2.24 GeV are determined by the KLOE Collaboration in φ → γa 0 (980)π → γπη [60], and g f0K + K − and g f0π + π − are based on the data of Ref. [61]. We note that the couplings from the KLOE Collaboration are chosen for consistency. Since the unitarized propagators for a 0 and f 0 are adopted, the couplings to KK should be extracted based on the same parametrization scheme. In fact, with the same couplings given in Tab. VI, significant difference between the Flatte propagator and the unitarized propagator can be identified. In our calculation the g η ′′ K * K and the g η ′′ a0π are the bare couplings determined by the combined analysis of the decays of η ′′ → KKπ, ηππ and 3π. Note that the corresponding VI: Coupling constants and phase angles determined in the combined analysis of η ′′ → KKπ, ηππ and 3π. The couplingg η ′′ a 0 π is defined with the TS contributions taken into account.
Parameters Values g f 0 K + K − 5.92 ± 0.13 GeV [61] g f 0 π + π − 2.96 ± 0.12 GeV [61] g a 0 K + K − 2.24 ± 0.11 GeV [60] ga 0 ηπ 3.02 ± 0.35 GeV [60] vertex appears in all these processes. The physical couplingg η ′′ a0π can be extracted by the combined contributions from the tree-level amplitude and TS mechanism in η ′′ → a 0 π. We also note that since experimentally the branching ratio of η ′′ → a 0 π is much smaller than that for η ′′ → KK * + c.c., we do not expect that the rescattering of η ′′ → a 0 π → KK * + c.c. has significant renormalization contributions to the coupling of g η ′′ K * K . Namely, we approximateg η ′′ K * K ≃ g η ′′ K * K . It should be cautioned that the determination of the couplingsg η ′′ K * K andg η ′′ a0π will strongly depend on the experimental accuracies of the corresponding branching ratios. For the coherent analysis done in this work it is sufficient for us to demonstrate the patterns arising from the TS mechanism. Therefore, we only list the central values for these two couplings in Tab. VI. Hence, we collect the amplitudes and express the total amplitudes for each channel as follows: The partial widths for η ′′ → f 0 π → π + π − π 0 , η ′′ → a 0 π → ηππ and η ′′ → KKπ channels are where Φ abc is the phase space of η ′′ → abc. For J/ψ radiative decay, the following expressions are adopted where Γ η ′′ ≃ Γ KKπ + Γ ηππ , m η ′′ = 1.42 GeV and m J/ψ = 3.1 GeV.
There are four types of couplings, corresponding to η ′′ → K * K , η ′′ → a 0 π, V → P P , and S → P P (S denotes a 0 or f 0 ). By taking the g V P P as a real coupling in the SU(3)-flavor symmetry, we introduce three energy independent phase factors e iφ1 , e iφ2 and e iφ3 to the couplings g η ′′ K * 0K 0 , g η ′′ a0π and g a0K 0K 0 , respectively. Namely, these phase factors are shared by the SU(3)-flavor multiplets because of flavor symmetry. For instance, in respect to the real coupling g V P P we express g a0K + K − = −g a0K 0K 0 = 2.24e iφ3 GeV and g f0K + K − = g f0K 0K 0 = 5.92e iφ3 GeV. By examining the transition processes of Figs. 7, 8 and 9, one notices that there are only two independent phase angles needed in the combined analysis. For convenience, we set φ 2 = 0, and redefine the two remained independent phase angles as the linear combination of φ 1 and φ 3 A feature arising from the experimental measurements of the η ′′ → KKπ and ηππ is that the branching ratio for KKπ is much larger than that for ηππ. In the combined analysis, such a difference is determined by the relative coupling strength between g η ′′ a0π and g η ′′ K * 0K 0 and the relative phase factors which are treated as parameters and to be constrained by experimental data. Our strategy of determining these couplings and relative phase angles is as follows: We first require that the relative branching ratio fractions between KKπ and ηππ satisfy the experimental The partial widths of η ′′ decays calculated for different amplitudes noted by the brackets. The second row of the experiment branching ratio is calculated by BR(J/ψ → γη ′′ → γa0π → γηπ 0 π 0 ) = (8.5 ± 1.75) × 10 −5 [42] and BR(J/ψ → γη ′′ → γKKπ) = (2.8 ± 0.6) × 10 −3 [54]. The first row of the experiment branching ratio is taken from Ref. [ data, and then the sum of the partial widths from these two channels can contribute to up to 50 MeV of the total width. Although both the total width and partial decay widths of η ′′ have large uncertainties, it is sufficient for understanding nearly all the existing puzzling questions about the η(1405/1475). With the cut-off parameter αΛ QCD = 500 MeV, i.e. Λ R = m R + 0.5 GeV, the couplings and phase angles are determined and listed in Tab. VI.

IV. RESULTS AND DISCUSSIONS
With the determined couplings in Tab. VI, the relative partial widths are listed in Tab. VII, where the dominance of the K * K (M tree1 ) channel in η ′′ decay is evident. For the KKπ final state the relatively small contributions from the intermediate a 0 π provide a small destructive cancellation which leads to Γ K * K /Γ KKπ = 1.06. This is also an indication that a partial wave analysis is need for the combined analysis of all these three decay channels into KKπ, ηππ and 3π. We will see later that the interference is more significant in the ηππ channel.
The KK and Kπ spectra at √ s = 1.42 GeV are shown in Fig. 10. In the left panel of Fig. 10, since K * K is dominant in the KKπ channel, the peak of the spectrum (red line) near 1 GeV is not generated by the pole of a 0 (980), but a projection of the peaks of K * andK * in the Dalitz plot into the KK spectrum. The K * peak is evident in the right panel of Fig. 10. It is worth mentioning that the intensity of triangle diagrams in KKπ channel is very small, although it becomes compatible with the tree diagram for η ′′ → a 0 π in η ′′ → ηππ.
For the η ′′ → ηππ channel the exclusive contributions from M tri2 and M tree3 are both larger than the experimental values. It requires a destructive interference between these two amplitudes which is crucial for understanding the rather small branching ratio for η ′′ → ηππ. To illustrate this destructive interference, the ηπ 0 mass spectrum at √ s = 1.42 GeV is plotted in Fig. 11 with the exclusive and coherent contributions from M tri2 and M tree3 . It is noticeable that the destructive interference leads to a nontrivial structure for the a 0 (980) lineshape. In particular, there is a double-peak structure arising from the K + K − and K 0K 0 thresholds of the triangle amplitude M tri2 . This prediction can be examined by future experiments with higher energy resolution.
It is helpful to understand why the contribution of triangle diagrams in ηππ channel is one order of magnitude larger than that in KKπ channel. On the one hand, this may be due to a larger phase space for η ′′ → ηππ. On the other hand, the most contributive part of a triangle diagram comes form the vicinity of the pole of a 0 (980), which is in the middle of the phase space for η ′′ → ηππ, but at the boundary for η ′′ → KKπ.
Due to the presence of the TS mechanism the effective couplingg ηa0π should include contributions from the tree diagram and the TS, which yieldsg η ′′ a0π = 0.5 GeV. It can be compared with the bare coupling g η ′′ a0π = 0.72 GeV extracted from the overall analysis. We also note that the overall analysis gives the bare coupling g η ′′ K * K = 3.97 which is slightly larger than that adopted in Refs. [45,46] since there is a weak destructive interference in the KKπ channel in our result. However, the dominance of the K * K in η ′′ → KKπ is evident.
Notice that partial wave analyses for these channels are still unavailable, and possible contributions from the ση channel to the ηππ final state are not excluded. Future improvement of the experimental data is strongly desired.
In the decay channel of η ′′ → 3π, the π + π − invariant mass spectrum at √ s = 1.42 GeV is depicted in Fig. 12. The typical narrow structure located between the charged and neutral KK thresholds is expected by the isospin violation mechanisms. In Fig. 12 the exclusive contributions from the isospin-violating TS diagram ( Fig. 7 (a)), and a 0 − f 0 mixing diagrams (Fig. 7 (b) and (c)) are plotted. Note that the a 0 − f 0 mixing after the triangle rescatterings appears to be more important than the a 0 − f 0 mixing through the tree diagram.
A unique feature with this channel is that the strength of the transition amplitude is quite insensitive to the cut-off energy. As discussed in Refs. [45,46], the dispersive parts of the charged and neutral kaon loop transition amplitudes largely cancel out, and the main contributions to the isospin violations are from the absorptive parts which are in the vicinity of the on-shell kinematic region. This explains that the width effects from the intermediate K * do not change the lineshape of the narrow structure of the ππ spectrum.
In Fig. 13 we show the spectra of η ′′ → KKπ, ηππ and 3π as a comparison of the different lineshapes for the same state in different channels. It confirms the main conclusion of Refs. [45,46]. Namely, due to the contributions from the TS mechanism the peak positions of the same state η ′′ are located in different masses in different channels. In particular, the spectrum of η ′′ in the KKπ channel appears to be asymmetric and apparently deviated from a Breit-Wigner form. This is the manifestation of the K * K dominance in η ′′ decays. As shown in Fig. 13 the shift of the peak positions due to the interferences from the TS mechanism can amount to almost 17 MeV.
It should be noted that the TS mechanism can also result in an energy-dependent description of the η ′′ K * K + c.c. coupling. This will be reported in a followed-up work. For the purpose of understanding the impact of the TS mechanism on the width effects from the intermediate K * and the relative strengths of exclusive decay channels for η ′′ , we actually fix the η ′′ K * K + c.c. coupling as a constant in the present calculations. dot-dot-dashed line is the spectrum with only triangle diagram of the f0π channel ( Fig. 7 (a)). The green dot-dashed line is the contribution from the a0 − f0 mixing with the bare g η ′′ a 0 π coupling (Fig. 7 (c)), and the orange dashed line is the contribution from the a0 − f0 mixing via the TS production of the a0(980) (Fig. 7 (b)). The black dotted line is the total contribution from the a0 − f0 mixing.

V. CONCLUSION
In this work we revisit the topical problem of η(1405/1475) in J/ψ radiative decays into KKπ, ηππ and 3π, where η(1405) and η(1475) are treated as being originated from the same state η ′′ but strongly affected by the triangle singularity mechanism. Different from the previous studies, we show that the non-zero width effects can be well-understood based on the analytical properties of triangle integral. We adopt the formulation by 't Hooft and Veltman [51] to introduce a complex mass for the internal state. The triangle loop integral can be expressed in terms of the Spence functions for which the analytical properties can be investigated. We then show that the triangle singularity conditions can be described by motions of a set of kinematic functions in the complex plane.
In the case where the internal states have a small width, the absorptive part of the scalar loop integral can be expressed in a very compact form, which is exactly the same as that obtained from the cutting rule. The simplification of the analytic absorptive part is difficult. But we make the conjecture that these two expressions are identical based on the numerical checks. To demonstrate this, we isolate the width dependence and prove that there is an typical energy scale β (determined by kinematic conditions) for the scalar loop integral, referring to which we can see the explicit influence of the finite width on the absorptive part.
Despite of the fact that the width of K * meson can be regarded as a small value compared with its mass, it is a big value compared with β. Thus, the large resulting suppression can be understood. We also suggest a possible way to get the accurate analytic absorptive part of scalar loop integral, which may be useful in the case where the width of internal masses is no longer a small value.
For the decays of η ′′ we obtain an overall description of η ′′ → KKπ, ηππ and 3π taking into account the K * width effects. It shows that the KKπ channel is the dominant one in its decays, and the tree diagram of η ′′ → KK * + c.c. is the main contribution. In this channel the TS contribution is relatively small in comparison with the tree process. For the decay channel of η ′′ → ηππ its small branching ratio can be understood by the destructive interference between tree and triangle diagrams. Note that the 3 P 0 quark model prediction for the branching ratios of η ′ (2S) → a 0 (980)π and KK * + c.c. are compatible [63]. Without the contributions from the TS mechanism it would be very difficult to understand the much weakened decay of η ′′ → a 0 (980)π as the first radial excited state of η ′ (958). An improved measurement of the branching ratio J/ψ → γη ′′ → γa 0 π → γηππ would be useful for further quantifying the TS contributions. For the decay of η ′′ → 3π we identify that the TS enhanced a 0 − f 0 mixing should also contribute to the isospin violating decay of η ′′ → 3π. Taking into account the width effects from the internal K * we obtain a consistent picture for understanding the role played by the TS concerning the η(1405)/η(1475) puzzle.
At this moment, there are still large experimental uncertainties with the measurement of the two-body decay widths. Therefore, partial wave analysis of these three channels combining the TS mechanism is strongly required. We would expect that experimental results from BESIII in the near future will be able to set all these controversial issues about η(1405)/η(1475).
The η(1405)/η(1475) puzzle has been a key issue for a better understanding of the pseudoscalar spectrum. For the purpose of searching for the pseudoscalar glueball state the settlement of this puzzle would provide a useful guidance for future efforts in both experiment and theory. Namely, our attention should be paid to the higher mass regime [45] and criteria for pseudoscalar glueball candidates should be further studied in detail.