Non-Abelian axial anomaly, axial-vector duality, and the pseudoscalar glueball

We generalize the dispersive approach to axial anomaly by A.D. Dolgov and V.I. Zakharov to a non-Abelian case with arbitrary photon virtualites. We derive the anomaly sum rule for the singlet current and obtain the $\pi^0,\eta,\eta'\rightarrow\gamma\gamma^*$ transition form factors. Using them, we established the behavior of a nonperturbative gluon matrix element $\langle 0 |G\tilde{G} |\gamma\gamma^*(q^2) \rangle$ in both spacelike and timelike regions. We found a significant contribution of the non-Abelian axial anomaly to the processes with one virtual photon, comparable to that of the electromagnetic anomaly. The duality between the axial and the vector channels was observed: the values of duality intervals and mixing parameters in the axial channel were related to vector resonances' masses and residues. The possibility of a light pseudoscalar glueball-like state is conjectured.


Introduction
One of the fundamental features of QCD, the axial anomaly, has many theoretical and phenomenological applications. Being essentially a nonperturbative phenomenon of QCD, the axial anomaly provides valuable insights into low-energy properties of the theory, inaccessible to perturbative methods. Among phenomenological manifestations of special importance are the γγ decays of neutral pseudoscalar mesons π 0 , η, η ′ . It was the problem of the neutral pion decay, which, in fact, led to the discovery of quantum anomalies [1,2]. The processes involving virtual photons, like transition form factors (TFFs) of pseudoscalar mesons, are also deeply connected with the axial anomaly. The existence of the axial anomaly results in the nonconservation of the axial current. In particular, a singlet axial current J (0) µ5 = (1/ √ 3) iψ i γ µ γ 5 ψ i gets both electromagnetic and strong (non-Abelian) anomalies, while isovector (a = 3) and octet (a = 8) components of the octet of axial currents J (a) µ5 = (1/ √ 2) iψ i γ µ γ 5 λ a ψ i acquire only the electromagnetic anomaly: Here F and G are electromagnetic and gluon field strength tensors respectively, F µν = 1 2 ǫ µνρσ F ρσ and G µν,t = 1 2 ǫ µνρσ G t ρσ are their duals, N c = 3 is a number of colors, α s is a strong coupling constant, C The axial anomaly can be calculated as a result of ultraviolet regularization of the diverging triangle VVA diagram. Alternatively, it can be derived considering the imaginary part of the VVA diagram [3], where the anomaly arises as a sum rule for a structure function in the dispersion representation of the three point VVA correlation function [4,5,6] (for a review, see [7,8]). Such sum rules, combined with the global quark-hadron duality, then can be employed to study the γγ decays of the pseudoscalar mesons as well as their transition form factors γγ * → π 0 , η, η ′ . In particular, the anomaly sum rule (ASR) for the isovector axial current was related to the pion TFF [9], while the anomaly sum rule for the octet axial current was developed and used to study the η, η ′ TFFs, including the cases of spacelike [10,11,12,13] and timelike [14] photon virtualities. These anomaly sum rules for the isovector (a = 3) and octet (a = 8) currents when one of the photons is real (k 2 = 0) and another is real or virtual (Q 2 = −q 2 ≥ 0) read (in what follows we put m u = m d = m s = 0) [5,9,10], where the spectral density function is defined as A ≡ 1 2 Im(F 3 − F 6 ) and is determined from a decomposition of the vector-vector-axial (VVA) amplitude T αµν (k, q) = e 2 d 4 xd 4 ye (ikx+iqy) 0|T {J α5 (0)J µ (x)J ν (y)}|0 (4) as [15] (see also [16,17]) T αµν (k, q) = F 1 ε αµνρ k ρ + F 2 ε αµνρ q ρ + F 3 k ν ε αµρσ k ρ q σ + F 4 q ν ε αµρσ k ρ q σ + F 5 k µ ε ανρσ k ρ q σ + F 6 q µ ε ανρσ k ρ q σ , where the coefficients F j = F j (p 2 , q 2 ), j = 1, . . . , 6 are the corresponding Lorentz invariant amplitudes constrained by current conservation and Bose symmetry. The electromagnetic currents are defined as J µ = i e iψi γ µ ψ i with momenta k, q. The rhs of (3) is exactly the Abelian (electromagnetic) anomaly constant stemmed from the matrix element 0|FF |γγ * . Due to appearance of the strong anomaly term in the singlet current (1), the analogous anomaly sum rule has an additional part stemmed from the matrix element The corresponding non-Abelian contribution in the dispersive form requires a subtraction, so the singlet current ASR in the considered kinematics N (p 2 , k 2 = 0, q 2 = −Q 2 ) ≡ N (p 2 , Q 2 ) can be written [18] as where Therefore, in addition to the Abelian anomaly contribution, the rhs of Eq. (7) has also a non-Abelian anomaly contribution (6) given by the subtraction N (0, Q 2 ) and the spectral parts. Rigorous calculation of these form factors face difficulty due to their nonperturbative nature. However, they can be related to the physical observables by means of the quark-hadron duality hypothesis. To carry this out, we saturate the lhs of ASRs (3) and (7) with a full set of resonances and single out the lowest-lying contributing states in each channel in terms of the corresponding transition form factors and decay constants. The "continuum" contribution absorbs the rest (higher resonances) with the same function A 3 (s, Q 2 ). So, we apply the following quark-hadron model, Here the sum is over the hadron states P whose decay constants f (a) P and the form factors F P γ of the transitions γγ * → P are defined as Then the ASRs (3) and (7) read, The lowest hadron contributions to the ASRs (12) and (13) are given by the π 0 , η and η ′ mesons, while the rest of the hadrons are absorbed by the integrals with lower limits s 3 , s 8 , s 0 -the duality intervals in the respective channels. The numerical values of these parameters will be discussed in the next section.
It is important to note that the singlet current is renormalization-scale dependent (unlike the isovector and octet currents) resulting in ∼ 15% decrease of f (0) P at high Q 2 , which was recently taken into account in analyses of the TFFs [19,20,21]. Nevertheless, in the ASR (13), the decay constants are fixed on the respective meson mass scales as a result of the δ-like resonance model in (9). Figure 1: The Feynman diagrams of the processes: (a) e + e − → e + e − P scattering, (b) P → e + e − γ Dalitz decays, (c) e + e − → P γ annihilation.
Recently, the ASR for the singlet current was used to study the role of strong and electromagnetic anomalies in the processes with real photons η, η → γγ [18]. We found that the strong anomaly contribution to these processes is small in comparison with the electromagnetic one.
In this paper, using the same framework, we are going to generalize the obtained results and consider the case of one virtual photon. Combining ASRs for all three currents and incorporating the quark-hadron duality, we derive relations between the pseudoscalar TFFs and matrix element 0|GG|γγ * . Then, using the experimental data for the TFFs, we study the properties of this matrix element at different photon virtualities.
The paper is organized as follows: in Sec. 2, the hadrons contributions and duality intervals are studied. Section 3 is dedicated to analysis of the strong anomaly contribution in the spacelike region. The analytic continuation of the ASRs to the timelike region is considered and the investigation of the strong anomaly properties in this region is carried out in Sec. 4. The last Sec. 5 summarizes the results.

Hadron contributions and duality intervals
The lower integration limits s 3 , s 8 , s 0 in (12), (13), emerging as free parameters of the ASR approach, are in fact the duality intervals (continuum thresholds) of the corresponding currents.
In the isovector and octet currents, the duality intervals s 3 , s 8 can be determined from the large Q 2 limit of the ASRs (3) by making use of the pQCD asymptotes of the π 0 , η, η ′ TFFs and numerical values of the decay constants f (3,8) P [10]. In this way, the interval of duality of the isovector channel is evaluated as s 3 ≃ 0.67 GeV 2 . However, this parameter can have a weak Q 2 dependence [51,52] varying from 0.6 GeV 2 at Q 2 → 0 to 0.67 GeV 2 at Q 2 → ∞. In the case of the octet channel, the numerical value of the corresponding interval of duality was found to vary s 8 = 0.4 − 0.6 GeV 2 , depending on a particular mixing scheme [11,13]. However, in the chiral limit (m q = 0) it is natural to assume that the duality intervals of the isovector and octet currents cannot be much different from each other: s 8 ≃ s 3 within 20% uncertainty of the SU (3) symmetry breaking. For the purposes of numerical analysis, we take s 8 ≃ s 3 = 0.6 GeV 2 in this paper.
The duality interval of the singlet current s 0 is different from s 3 and s 8 . The main hadron contribution here is given by the η ′ meson, which is known to retain its mass in the chiral limit, contrary to the lighter π 0 and η mesons. So, the value of s 0 should be of order of m 2 η ′ and it is natural to assume that s 0 ≃ 1 GeV 2 . Although it is possible that the value of s 0 can have a weak Q 2 dependence, similar to that of s 3 .
The one-loop approximation for the spectral densities of the isovector and octet currents A (s, Q 2 ) (Fig. 2) in the chiral limit is given by [17] so the integration in the ASRs (12) leads to the following expressions for the hadron contributions, The case of the singlet current differs from the isovector and octet currents due to a new type of diagrams involving virtual gluons. In order to single out electromagnetic contribution, we split the spectral density into two parts, A The second part A QCD is the contribution of diagrams ( Fig. 3) with virtual gluons coupled to two photons through nonperturbative strong interactions (see also [18]). The first part A (0) QED represents the contribution of QED diagrams, whose lowest one-loop part ( Fig. 2) is given by a similar expression to Eq. (14) with an appropriate charge factor C (0) . Making use of it, we can rewrite the ASR (13) as The first and the last three terms in Eq. (17) represent the electromagnetic and the strong anomaly contributions to the ASR respectively. It is convenient to introduce a function that represents the ratio of contributions of strong and electromagnetic anomalies: As the integral of A QCD is suppressed as α 2 s at s 0 ≥ 1.0 GeV 2 , the function B(Q 2 , s 0 ) is predominantly determined by the first two terms. It reflects the properties of the nonperturbative matrix element 0|GG|γγ * . 1 Therefore, the study of the function B(Q 2 , s 0 ) gives us access to the non-Abelian anomaly contribution to the processes P → γ * γ at various photon virtualities, which is the main goal of the present paper.
As B(Q 2 , s 0 ) cannot be calculated perturbatively, we will extract it from the experimental data. This resembles the extraction of nonperturbative ingredients, like (transverse momentum dependent) parton-distribution functions in QCD and (nonlocal) vacuum condensates from data or lattice simulations [53]. In a sense, we are going to solve the inverse problem.
So, rewriting the ASR for the singlet current (17) in terms of the function B(Q 2 , s 0 ) (18), we get, Taking into account the lowest contributions, given by the π 0 , η, and η ′ mesons, the ASRs for the isovector, octet (15), and singlet (19) currents comprise a system of equations, whose solution leads to the expressions for the form factors, where P = π 0 , η, η ′ . The coefficients α P , β P , γ P are expressed in terms of the decay constants f (i) P , ∆ is the determinant of the decay constant matrix in (20): So, the expressions for the TFFs (21) are a consequence of the dispersive approach to axial anomaly and quark-hadron duality. In this way, it provides theoretical grounds for the Brodsky-Lepage interpolation formula for pion TFF [54] and to some of its generalizations to the η and η ′ TFFs [55].
The decay constants f (i) P are defined as projections of the currents onto the meson states (10). Since the meson states are not pure SU (3) states, the corresponding decay constant matrix has nonzero off-diagonal elements. In the pseudoscalar sector, η and η ′ mesons manifest the largest mixing, which results in substantial values of f (8) η ′ , f (0) η (see, e.g., [56,57,58,59,13,21]). Different mixing schemes used to determine these decay constants imply that the decay constants follow a particular symmetry (e.g., quark-flavor or octet-singlet) and effectively impose different restrictions on the decay constants [13], reducing the number of parameters which describe the mixing.
The pion dominates in the isovector current with f π 0 = f π = 0.1307 GeV, while its contribution to the octet and singlet currents as well as η and η ′ contributions to the isovector current are small: f [60,61,62,63]. Although in many cases the effects of π 0 mixing with the η − η ′ system can be neglected, they are found to be relevant in the case of timelike region [52]. In the present paper, we will confirm this observation.
Let us investigate the properties of the function B(Q 2 , s 0 ) (18), which represents the ratio strong to electromagnetic contributions, in the spacelike and the timelike regions using the experimental data for the TFFs. Our main goal is to establish the qualitative behavior of it, and as far as possible, to get its quantitative estimates too.

Spacelike region
The pion TFF in the spacelike region (21) is described predominantly by the first term, while the second and the third terms can be neglected as the η, η ′ admixtures to the isovector current are small (α π ≫ β π , γ π ; see Table 1). For this reason, the π 0 TFF does not depend on B(Q 2 , s 0 ) and coincides with the expression obtained earlier [9,13]. It is in fair agreement with the CELLO [23], CLEO [24], BaBar [25] and Belle [27] experimental data.
Contrary to the pion TFF, B(Q 2 , s 0 ) plays an important role in the η and η ′ TFFs. As one can see from Table 1, the η ′ TFF (21) is almost completely determined by the third term (γ η ′ ≫ β η ′ , α η ′ ), while the η TFF has a significant contribution from the second term (β η ∼ γ η ). Therefore, when evaluating B(Q 2 ), to avoid additional uncertainty from s 8 , we will rely on the η ′ TFF data and check the validity of the results for the η meson.

Evaluation of B(Q 2 , s 0 )
The case of two real photons k 2 = Q 2 = 0 was considered earlier [18], where the suppression of the strong anomaly contribution was found. In this limit the form factors are expressed in terms of the two-photon decay widths, F P γ (0) = 64πΓP →2γ e 4 m 3 P . The results for B(0) evaluated from F η ′ γ for different decay constant sets are given in Table 2. The obtained values B(0) ≪ 1 reaffirm our previous observation: for real photons, the strong anomaly contribution appears to be much smaller than the electromagnetic one, as B(Q 2 ) effectively represents their ratio.  Table 2: The gluon to electromagnetic anomalies ratio B(0, s 0 )| s0∼1GeV 2 in the case of two real photons for different mixing schemes (analyses). However, the assumption that the strong anomaly contribution is negligible in the case of a nonzero photon virtuality Q 2 > 0 does not hold. The plots of the η, η ′ TFFs (21) with B(0) from Table 2 (see Figs. 4,5) show a strong disagreement with the experimental data, pointing out that the matrix element 0|GG|γγ * should be significant. One can also observe that Q 2 |F η,η ′ (Q 2 )| tends to a constant value starting from Q 2 ∼ 5 − 10 GeV 2 which indicates that B(Q 2 ) is already close to its asymptotic value in this region.
Mix. sch.  Table 3: Results of the fit of η ′ TFF to the experimental data at different intervals of Q 2 .
As one can see from Table 3, the values of B(Q 2 ) at Q 2 > 0.6 GeV 2 considerably differ from B(0) ∼ 0 at Q 2 = 0 and are close to their asymptotical values B as in the region Q 2 > 0.6 GeV 2 . The consistency of the B values in different intervals supports our assumption that it is constant in this region. Note also, that the values B as − B(0) appear to be close in different mixing schemes.
We can conclude, that at Q 2 0.6 GeV 2 the strong anomaly comprises ∼ 25% of the electromagnetic anomaly contribution and is almost independent of the photon virtuality, while at smaller Q 2 it rapidly vanishes. Taking into account that the available data at Q 2 < 0.6 are scarce and have large uncertainties, it is reasonable to test a simple approximation of the function B(Q 2 ) by a step-function: B = B(0) at 0 Q 2 < 0.6 GeV 2 and B = B as at Q 2 0.6 GeV 2 . The η and η ′ TFFs for different mixing parameters is compared with the data in Figs. 6 and 7. The green shaded area shows the 20% uncertainty of s 8 (for the EGMS16 [21] mixing scheme).   The analysis in the previous section was carried out for s 0 = 1 GeV 2 . In order to examine the sensitivity of our approach to s 0 , we perform a combined fit of the η and η ′ TFFs with two parameters s 0 and B to the available data [22,23,24,26] using TMinuit [64]. The obtained values of s 0 , B, and the corresponding χ 2 η+η ′ /dof for different mixing schemes are listed in Table 4. As one can see, the obtained value s 0 ≃ 1 GeV 2 is consistent in all considered mixing schemes.   The contour plot of χ 2 η+η ′ /dof for the mixing parameters EGMS16 [21] is shown in Fig. 8. The results for other mixing schemes are similar. The black filled circle and the black filled square indicate the minima of χ 2 from the combined η+η ′ and separate η ′ data fits respectively; they show rather good agreement with each other. This plot demonstrates a correlation between B and s 0 . In particular, though the best-fit values are s 0 ≃ 1 GeV 2 , B ≃ −0.25, it allows B ∼ 0 provided s 0 ≃ s 3 ≃ 0.6 GeV 2 , which means a negligible role of the strong anomaly at larger Q 2 , not only at Q 2 = 0. We will call this scenario (s 0 ≃ 0.6 GeV 2 and B(Q 2 ) ≃ 0) a hidden strong anomaly case, as opposite to the open strong anomaly case (s 0 ≃ 1 GeV 2 and B(Q 2 ) ≃ −0.25). Further study in the timelike region will show that these regimes correspond to different kinematical regions.
The hidden anomaly scenario has a physical interpretation. The mass of the dominant contributor to the singlet ASR -the η ′ meson -originates from the strong anomaly (U (1) A problem) [65,66]. Neglecting the strong anomaly brings the value of s 0 close to the values of the octet and isovector intervals of duality s 8 and s 3 , dominated by the η and π 0 mesons. Also, the case of hidden anomaly corresponds, in some sense, to the hypothesis [58,26] (see also [13]) that an unphysical state (1/ √ 2)(|ūu + |dd ) differs only by a charge factor from the isovector state (1/ √ 2)(|ūu − |dd ). As we already mentioned, the singlet current is not renormalization invariant, which results in the renormalization-scale dependence of the corresponding singlet constants f (0) P . The model of quarkhadron duality that we apply (9) leads to the decay constants fixed on the mass scales of the respective mesons in the singlet ASR. However, more refined models could result in Q 2 -dependence of the scale of these constants in the ASR. Numerically, f (0) [19,20,21]. In this case, choosing the scale as large as 49 GeV 2 (exceeding the currently available experimental data range), results in B as | µ 2 =49GeV 2 = 0.86B as | µ 2 =1GeV 2 which is well within our expected theoretical uncertainty of 20 %.

B and s 0 from Padé approximants
It is instructive to test our results using available approximation formulas for the TFFs. As one of such recent descriptions, we will use the one that employs the Padé approximants [67,20,21]. In this approach, the TFFs are described by means of rational functions, The coefficients of the approximants -P 1 1 (Q 2 ) for π 0 [67], η ′ [21] and P 2 2 (Q 2 ) for η [20] -are listed in Appendix B. The values of f (0) π 0 ,η,η ′ were employed from the EGMS16 [21] mixing scheme. Substituting the TFFs (25) directly into the ASR for the singlet current (19), we evaluate the s 0 as a function of Q 2 . The results for s 0 (Q 2 ) for the hidden (B = 0) and open (B = B as = −0.26) anomaly cases are shown in Fig. 9 by dashed green and solid red curves respectively. One can observe that the function s 0 (Q 2 ) demonstrates close to constant behavior in both cases, confirming our evaluations: s 0 ≃ 0.6 is close to s 3,8 in the hidden anomaly scenario, while s 0 ≃ 1 GeV 2 in the open anomaly scenario.
So we can see that the comparison with the Padé approximations of the TFFs [67,20,21] confirm our conclusions made in the previous subsections.

Timelike region
The analytic continuation of the the ASRs for isovector and octet currents (3) to the timelike region was carried out earlier [14]. It was shown that for Eqs. (15), the analytic continuation to the timelike region leads to the following relations for the real parts of the TFFs, The corresponding imaginary parts of the TFFs are negligible (implying |F P (q 2 )| ≃ |ReF P (q 2 )|) everywhere except the vicinity of q 2 = s 3,8 [14]. The analytic continuation of the ASR for the singlet current (7) is not so straightforward. The rhs of Eq. (7) contains functions depending on Q 2 appeared due to the strong anomaly matrix element, also the spectral density A  (7) is beyond the scope of this paper, however, we can try to make it directly from the equation for the TFFs (19). The Eq. (19) differs from the analogous relations for the isovector and octet currents (15) by the unknown function B(Q 2 ) in the rhs In the previous chapter we established that for the spacelike photon virtualities (Q 2 = −q 2 > 0) the function B(Q 2 ) rapidly tends to a constant value. The asymptotic behavior of the TFFs (21) and B(Q 2 ) for the timelike and spacelike photon virtualities should coincide. Therefore, at high q 2 one can suppose that the analytic continuation for the singlet case is analogous to the isovector and octet ones, and so it results in a replacement of Q 2 with −q 2 , We will try to describe the TFFs in the form of Eq. (27) at all q 2 , not just at high q 2 . At the same time, at intermediate region, B(q 2 ) can be probably more complicated than in the spacelike region, e.g., it may not be monotonic. The particular behavior of B(q 2 ) will be extracted from the data. Solving Eqs. (26) and (27), we obtain the expressions similar to Eqs. (21), but for the real parts of the TFFs in the timelike region. Therefore, with the exception of the vicinities of q 2 = s 3 , s 8 , s 0 , we can write for the absolute values of the timelike TFFs as where P = π 0 , η, η ′ . This equation looks similar to that appearing in the vector meson dominance (VMD) model and provides a complementary justification for it based on the dispersive approach to axial anomaly (see Sec. 4.4).
Because of such structure of the TFFs (28), even small coefficients α P , β P , γ P become relevant in the timelike region, in contrast to the spacelike one. This appears to be of special importance in the pion TFF: the third term with γ π 0 ≪ α π 0 cannot be neglected in this case. Moreover, as we will see, γ π 0 = 0 is necessary for the correct description of the data, which also conforms recent analyses indicating a nonzero mixing of π 0 with the η − η ′ system [63]. Let us note, that a similar form of the transition form factor in timelike was also employed in Ref. [52].
For the purposes of evaluation of the TFFs (28), we will use the decay constants from Refs. [21,63] (denoted as EGMS16 in Appendix A).

Dalitz decays domain
Let us start our study of the η and η ′ mesons in the timelike region from the Dalitz decays domain (4m 2 l ≤ q 2 < m 2 P ). The η meson experimental data are available in two decay modes: a dielectron mode, measured by A2 collaboration [29], and a dimuon mode, measured by NA60 collaboration [30], covering the region of 0 < q 2 < 0.25 GeV 2 . Concerning the η ′ meson, the experimental data are available only in the dielectron mode from the BESIII collaboration [31] in the region 0 < q 2 < 0.56 GeV 2 .
In the previous section we discussed two different cases: the hidden anomaly (negligible strong anomaly contribution and s 0 ≈ s 3,8 ) and the open anomaly (significant strong anomaly contribution and s 0 ≃ 1 GeV 2 ). The TFFs (28) for both cases (solid red and dashed green curves respectively) are shown in Figs. 11 and 12; the shaded area indicates the estimated theoretical uncertainty (contributed by s 8 ) of 20%. The respective χ 2 values and the curve slopes at q 2 = 0, a P ≡ lim q 2 →0 ∂ ∂q 2 |FP (q 2 )| 2 |FP (0)| 2 , are given in Table 5.  (28) for different s 0 . The experimental data are by A2 [29] and NA60 [30] collaborations. The experimental data are by BESIII collaboration [31].  Table 5: χ 2 /dof and curve slopes a P of the timelike η, η ′ TFFs (28) for different s 0 .
One can see from Figs. 11, 12 and Table 5, that the data in the considered small-q 2 region support the hidden anomaly case. Let us note that this conclusion (B = 0 at low-q 2 ) remains valid also for the pion TFF, see the insert in Fig.13. The value of the slope a π 0 = 1.639 GeV −2 also agrees with the experimental value of the A2 collaboration 1.646 ± 0.549 GeV −2 [28].
The function B(q 2 ) at large q 2 in the timelike region should coincide with that in the spacelike region, B(q 2 → ∞) = B as , while at q 2 = 0, B(0) ≃ 0. Let us suppose that it has the same form as the estimate obtained earlier in the spacelike region, B(q 2 ) = B(0) = −0.024 at q 2 < 0.6 GeV 2 and B(q 2 ) = B as = −0.26 at at q 2 > 0.6 GeV 2 (results from Tables 2,3 for the EGMS16 mixing scheme). Employing the same values for the duality intervals as in the spacelike region, s 8 = s 3 = 0.6 GeV 2 , s 0 = 1 GeV 2 at q 2 > 0.6 GeV 2 , we plot the corresponding π 0 TFF in Fig. 13 (dashed green curve) and compare it to the data [32,36]. The insert shows the results for the Dalitz decay domain [28].
We see that the last term in (28) is necessary for reproducing the second peak in the data, implying that γ π 0 is nonzero and even small π 0 -η-η ′ mixing has impact here. For comparison, the dot-dashed blue curve shows indicates the limit γ π 0 = 0, i.e. when π 0 and η-η ′ mixing is neglected.
However, the vicinity of the second peak (at q 2 ≃ 1 GeV 2 ) is described incorrectly, so in order to improve the data description we need to adjust the function B(q 2 ). The interference pattern of the data indicates that the first two and the third pole terms in the TFF should have different signs in this region, which can be satisfied only if 1 + B(q 2 ) becomes negative at q 2 ≃ 1 GeV 2 . Therefore, B(q 2 ) has an extremum in this region B < −1. As a modification of B(q 2 ), we suggest a Gaussian function where B as = −0.26 is the asymptotic value, and b, µ, c are the parameters. Numerically, the parameters were found to vary, −1.25 b < −1.00, 0.1 c 0.2 GeV, µ ≃ 1 GeV (we took into account the η meson TFF constraints, which will be analyzed in the next section). The π 0 TFF for the parameters b = −1.1, c = 0.15 GeV, µ = 1.0 GeV is shown in Fig. 13 as a solid red curve. One can see that this modified B(q 2 ) significantly improves the data description. We also considered a different modification of B(q 2 ), in the form of a step function. It gives a worse description, especially in the case of the η TFF (dashed green curves in Figs. 13 and 14). The suggested Gaussian modification is also more physically motivated, providing a model for nonperturbative behavior of the non-local condensates [68]. The similar short-range behavior can be observed in the transverse momentum dependent parton distributions [69,70].  Figure 13: π 0 TFF (28) compared with the experimental data [28,32,36]. The solid red curve -TFF with a Gaussian B(q 2 ) (29), the dashed greed curve -TFF with a step B(q 2 ), the dot-dashed blue curve -mixing of π 0 and η − η ′ is neglected. The insert shows the Dalitz decay domain. Therefore, we can conclude that the pion TFF experimental data in the timelike region confirms the existence of π 0 − η − η ′ mixing. We found that the function B(q 2 ) should have a sharp minimum at q 2 ≃ 1 GeV 2 with B(q 2 ≃ 1) < −1.
We can conclude that the proposed B(q 2 ) with Gaussian minimum at q 2 ≃ 1 GeV 2 (29) gives a consistent description of the π 0 , η TFFs in a wide range of q 2 , providing correct reproduction of the data at small and large q 2 as well as the interference pattern near the poles.
The explored Gaussian-like minimum of B(q 2 ) can mean a resonance contribution in the gluon matrix element, ∞ 0 ImR(s, Q 2 )ds = X 0|GG|X X|γγ * . From the fact that the minimum is at q 2 ≃ 1 GeV 2 one can suppose that it is provided by pseudoscalar resonance(s) 2 with masses ∼ 1 − 1.3 GeV, so the virtual photon energy in the final state could be ∼1 GeV. A light glueball-like state is a possible candidate. It is also possible, that this minimum is a result of a sum of several resonances. Let us note, that some approaches [71,72,73,74,75] predict the lowest glueball state with the mass ∼ 1.4 GeV, which is close to this region. However, higher mass resonances are also often considered as candidates for the pseudoscalar glueball (see, e.g., [76,77]).
Finally, let us consider the high-q 2 region, where the η and η ′ TFFs were measured at the point q 2 = 112 GeV 2 by BaBar [37]. The results for the η and η ′ TFFs (28) with B = B as , s 0 = 1 GeV 2 , s 3 = s 8 = 0.6 GeV 2 at this point for different mixing schemes are listed in Table 6. We can see a good description of the η ′ data. For the case of η meson the results are consistent with the experiment at the order of 2σ.

Axial-vector duality
As we saw in the previous section, the values of the duality intervals s 3,8,0 evaluated in the spacelike region are directly related to the positions of the peaks in the timelike region of the TFFs. It is worth noting that the processes e + e − → π 0 (η)γ can be described within the framework of the VMD model, which implies the intermediate vector meson contributions to the photon propagator [78,79]. In such a model, a timelike pseudoscalar meson TFF is represented by a sum over the lightest vector meson contributions (ρ, ω, φ).
Comparing our results with the VMD model, we observe correlation between the quantities inherent to the axial channel (and pseudoscalar hadron resonances) and vector hadron resonances, observed earlier in the octet channel [13,14]: s 3 , s 8 , s 0 ←→ m 2 ρ , m 2 ω , m 2 φ . In this paper we additionally established that the mixing parameters are important for the data description in the vicinity of the peaks. So, the decay constants of pseudoscalar mesons (which determine α P , β P , γ P coefficients) correspond to the vector meson coupling constants (residues) of the VMD model. This means, in particular, that the large η-η ′ mixing in the pseudoscalar sector is correlated with the residues of the vector mesons.
These observations confirm that the axial anomaly in its dispersive form (i.e. respective ASRs) reveals the duality between axial and vector channels. This duality is also related [80] to the theorems [81] for longitudinal and transverse parts of two-point VA correlators in external electromagnetic fields. There are also relations between resonances in different channels in the holographic approach [82].

Conclusions
We considered the non-Abelian axial anomaly in the framework of dispersive approach to derive the anomaly sum rule for the singlet current (7). By means of this ASR and the ASRs for the isovector and octet currents we obtained the transition form factors of the π 0 , η, η ′ mesons in spacelike (21) and timelike (28) regions. In order to investigate the contribution of the strong (gluon) anomaly, stemming from the matrix element 0|GG|γγ * , to these processes, we introduced the function B(q 2 , s 0 ), which is a ratio of the strong and electromagnetic anomalies. This function provides the description of the nonperturbative vacuum structure and may be also studied in lattice simulations.
We confirmed our earlier observation [18] that at low photon virtuality |q 2 | the role of the strong anomaly in the corresponding processes is negligible (B(q 2 , s 0 ) is close to zero). At larger |q 2 | in both, spacelike and timelike regions, the relative strong anomaly contribution rapidly reaches its asymptotic value B ≃ −0.25. However, in contrast to the spacelike region, in the timelike region it should have a sharp extremum (minimum B ≃ −1.3) at q 2 ∼ 1 GeV 2 (Fig.15), so the strong anomaly contribution exceeds the electomagnetic one. This might be attributed to glueball-like pseudoscalar resonances.
We observed the correlation between the value of the duality interval for the singlet current s 0 . The analysis showed, that when B ≃ 0, the duality interval of the singlet current becomes close to the values of the duality intervals of the octet and isovector currents, s 0 ≃ s 3 ≃ s 8 ≃ 0.6 GeV 2 . This (hidden anomaly) scenario occurs in the region of small |q 2 |, while at larger |q 2 | the relative strong anomaly contribution B becomes significant while s 0 reaches 1 GeV 2 (open anomaly).
From our analysis we found that the mixing of the pion with the η and η ′ mesons is necessary for correct description of the pion timelike TFF. This is in agreement with recent estimations of mixing parameters from V P γ decays, where a nonzero isospin-symmetry breaking was evaluated [63].
The axial anomaly in its dispersive form reveals the duality between axial (pseudoscalar) and vector sectors: the intervals of duality of the axial channel were found to be numerically close to the vector meson masses, while the combinations of the mixing parameters and the intervals of duality correspond to the behavior near the poles and the coupling constants (residues) of the VMD model.