Two-photon decays and transition form factors of $\pi^0$, $\eta$, and $\eta'$ in large-$N_c$ chiral perturbation theory

We present a calculation of $P\to \gamma^{(\ast)}\gamma^{(\ast)}$ processes, where $P=\pi^0,\ \eta,\ \eta'$, at the one-loop level up to and including next-to-next-to-leading order (NNLO) in large-$N_c$ chiral perturbation theory. The results are numerically evaluated successively at LO, NLO, and NNLO. The appearing low-energy constants are determined through fits to the available experimental data. We investigate the decay widths to real photons, the single-virtual transition form factors, and the widths of $P\to\gamma l^+l^-$, where $l=e,\ \mu$. Furthermore, we provide results for the slopes and curvatures of the transition form factors.


I. INTRODUCTION
In recent years, the two-photon interaction of the light pseudoscalar mesons has received considerable attention from both the experimental and theoretical sides [1]. To a large extent, this renewed interest was triggered by the muon anomalous magnetic moment discrepancy which states a 3.3 sigma deviation between experiment and theory (see Refs. [1-3] for reviews). On the theoretical side, the largest uncertainty in the anomalous magnetic moment a µ originates from the evaluation of hadronic contributions, namely, the hadronic vacuum polarization and the hadronic light-by-light (HLbL) scattering [3]. In this context, the two-photon decays of the light pseudoscalars enter the HLbL contribution in terms of pseudoscalar-exchange diagrams (see Fig. 35 of Ref. [2]).
Besides this more phenomenology-driven interest, the two-photon decays of the light pseudoscalars provide an ideal laboratory for investigating the symmetry-breaking mechanisms relevant in quantum chromodynamics (QCD). To be specific, the low-energy regime of QCD is characterized by an interplay between the dynamical (spontaneous) breaking of chiral symmetry, the explicit symmetry breaking by the quark masses, and the axial U(1) A anomaly. It is a generally accepted feature of QCD that the global U(3) L × U(3) R chiral symmetry of the QCD Lagrangian at the classical level for vanishing up-, down-, and strangequark masses is dynamically broken down to SU(3) V × U(1) V in the ground state (see, e.g., Ref. [5] for a discussion). Naively, one would then expect the appearance of nine massless pseudoscalar Goldstone bosons [6]. However, because of quantum effects, the singlet axialvector current is no longer conserved (U(1) A anomaly) and the corresponding alleged singlet Goldstone boson acquires a mass even in the chiral limit of massless quarks [7][8][9]. At this stage, the large-number-of-colors (LN c ) limit of QCD [10,11], i.e., N c → ∞ with g 2 N c fixed, provides another theoretical simplification aside from the assumption of massless quarks. Since the four divergence of the anomalous singlet axial-vector current is proportional to the square of the strong coupling constant g [12], it vanishes in the LN c limit. Therefore, the singlet pseudoscalar is also a Goldstone boson in the combined chiral and LN c limits, resulting in total in a pseudoscalar nonet (π, K, η 8 , η 1 ) as the Goldstone bosons [8,13]. Of course, massless LN c QCD is only an approximation to the real world. However, one may use it as a starting point for a perturbative framework, treating the symmetry breaking by the U(1) A anomaly and by the nonzero quark masses as corrections.
At leading order in the 1/N c and quark-mass expansion, the decays P → γγ (P = π 0 , η, η ) are driven by the chiral anomaly in terms of the Wess-Zumino-Witten effective action [14,15] (see, e.g., Ref. [5] for an introduction). Corrections to the WZW predictions originate from the axial U(1) A anomaly and the nonzero quark masses. Both mechanisms are also reponsible for generating the masses of the originally massless Goldstone bosons and for the η-η mixing. These modifications may be systematically calculated in the framework of large-N c chiral perturbation theory (LN c ChPT) [16][17][18], which can be viewed as an extension of conventional ChPT [19] by including, in addition to the pseudoscalar octet, the pseudoscalar singlet. In LN c ChPT, the most general effective Lagrangian is organized in a combined expansion in terms of momenta (derivatives), quark masses, and 1/N c . Observables are calculated perturbatively, according to a power counting with respect to a collective small expansion parameter δ [16].
In this article, we investigate the P → γ ( * ) γ ( * ) interaction at next-to-next-to-leading order (NNLO) in LN c ChPT. In Sec. II, we describe the effective field theory we will consider for our calculation by specifying the Lagrangian and the power counting. In Sec. III, we define the invariant amplitude and discuss its perturbative calculation including the η-η mixing. Section IV contains our numerical results for the decay rates and the transition form factors (TFF) at LO, NLO, and NNLO, respectively. In Sec. V, we discuss the decay rates for single Dalitz decays. Finally, in Sec. VI we conclude with a few remarks and an outlook on possible future work.

II. LAGRANGIANS AND POWER COUNTING
In the framework of LN c ChPT, one performs a simultaneous expansion of (renormalized) Feynman diagrams in terms of momenta p, quark masses m, and 1/N c . 1 Introducing a collective expansion parameter δ, the variables are counted as small quantities of order [16] The most general Lagrangian of LN c ChPT is organized as an infinite series in terms of derivatives, quark-mass terms, and, implicitly, powers of 1/N c , with the scaling behavior given in Eq. (1): where the superscripts (i) denote the order in δ.
The dynamical degrees of freedom are collected in the unitary 3 × 3 matrix where the Hermitian 3 × 3 matrix contains the pseudoscalar octet fields and the pseudoscalar singlet field η 1 , the λ a (a = 1, . . . , 8) are the Gell-Mann matrices, and λ 0 = 2/3 1. In Eq. (3), F denotes the piondecay constant in the three-flavor chiral limit and is counted as [8]. In addition to the dynamical degrees of freedom of Eq. (4), the effective Lagrangian also contains a set of external fields (s, p, l µ , r µ , θ). The fields s, p, l µ , and r µ are Hermitian, color-neutral 3 × 3 matrices coupling to the corresponding quark bilinears, and θ is a real field coupling to the winding number density [19]. The external scalar and pseudoscalar fields s and p are combined in the definition χ ≡ 2B(s + ip) [19]. The low-energy constant (LEC) B is related to the scalar singlet quark condensate qq 0 in the three-flavor chiral limit and is of O(N 0 c ) [16]. In general, applying the power counting of Eq. (1) to the construction of the effective Lagrangian in the LN c framework involves two ingredients. On the one hand, there is the momentum and quark-mass counting which proceeds as in conventional SU(3) ChPT [19]: (covariant) derivatives count as O(p), χ counts as O(p 2 ), etc. We denote the corresponding chiral order by D p . The LN c behavior can be determined by using the following rules (see Refs. [17,18] for a detailed account). In the LN c counting, the leading contribution to a quark correlation function is given by a single flavor trace and is of order N c [10,11,20]. In general, diagrams with r quark loops and thus r flavor traces are of order N 2−r c . Terms without traces correspond to the purely gluonic theory and count at leading order as N 2 c . This argument is transferred to the level of the effective Lagrangian, i.e., single-trace terms are of order N c , double-trace terms of order unity, etc. 3 Introducing ψ = √ 6η 1 /F [18], each power (ψ + θ) n is accompanied by a coefficient of order O(N −n c ). The reason for this assignment is the fact that, in QCD, the external field θ couples to the winding number density with strength 1/N c . In a similar fashion, D µ θ (as well as multiple derivatives) are related to expressions with O(N −1 c ). 4 Denoting the number of (ψ + θ) and D µ θ terms by N θ , the LN c order reads [17,18] The combined order of an operator is then given by A. Wess-Zumino-Witten effective action The two-photon decays arise from the odd-intrinsic-parity part of the effective field theory. At leading order, they are driven by the chiral anomaly, which is accounted for by the Wess-Zumino-Witten (WZW) action [14,15]. In U(3) ChPT, the WZW action (without external fields) reads where . . . denotes the (flavor) trace. For the construction of the WZW action, the domain of definition of U needs to be extended to a (hypothetical) fifth dimension, where Minkowski space is defined as the surface of the five-dimensional space for α = 1. The indices i, . . . , m in Eq. (7) run from 0 to 4, y 4 = y 4 = α, ijklm is the completely antisymmetric (five-dimensional) tensor with 01234 = − 01234 = 1, and U L i = U † ∂U/∂y i . In the presence of external fields, the anomalous action receives an additional term [21,22] 3 When applying these counting rules, one has to account for the so-called trace relations connecting singletrace terms with products of traces (see, e.g., Appendix A of Ref. [25]). 4 Note that we do not directly book the quantities (ψ + θ) or D µ θ as O(N −1 c ), but rather attribute this order to the coefficients coming with the terms.
given by with Z µνρσ (U, l, r) where U Lµ ≡ U † ∂ µ U and U Rµ ≡ U ∂ µ U † . The subtraction of the Z µνρσ (1, l, r) term is necessary to satisfy a boundary condition leading to an action that is consistent with the conservation of the vector current.

B. Normal-parity Lagrangians
In the NNLO calculation of the two-photon decays, the LO, NLO, and NNLO Lagrangians of even intrinsic parity enter as well. The leading-order Lagrangian is given by [16,18] where the covariant derivatives of U and U † are defined as The constant τ = O(N 0 c ) is the topological susceptibility of the purely gluonic theory [16]. Counting the quark mass as O(p 2 ), the first two terms of L (0) are of O(N c p 2 ), while the third term is of O(N 0 c ), i.e., all terms are of O(δ 0 ). The normal-parity part of the NLO Lagrangian L (1) was constructed in Refs. [16][17][18] and receives contributions of O(N c p 4 ), O(p 2 ), and O(N −1 c ). We only display the terms relevant for our calculation, in particular, here, we set v µ ≡ (r µ + l µ )/2 = 0 and a µ ≡ (r µ − l µ )/2 = 0 in the covariant derivatives: where and the ellipsis refers to the neglected terms. The first two terms of L (1) count as O(N c p 4 ) and are obtained from the standard SU(3) ChPT Lagrangian of O(p 4 ) [19] by retaining solely terms with a single trace and keeping only the constant terms of the so-called potentials which are functions of ψ + θ [18]. According to Eq. (15), the expression D µ ψD µ ψ implicitly involves two flavor traces (see footnote 7 of Ref. [18]), with the result that the corresponding term is O(N 0 c ). The SU(3) Lagrangian of O(p 6 ) was discussed in Refs. [23][24][25][26], and the generalization to the U(3) case has recently been obtained in Ref. [27]. For the present purposes, at NNLO, the relevant pieces of L (2) can be split into three different contributions of O(N −1 c p 2 ), O(p 4 ), and O(N c p 6 ), respectively: where The coupling v

C. Odd-intrinsic-parity Lagrangians
The WZW term accounts for the anomaly. Beyond the WZW action, the terms of the odd-intrinsic-parity sector are ordinary local Lagrangians expressible as closed expressions in U . However, also in this case, the U(3) unnatural-parity Lagrangian contains additional terms in comparison with its SU(3) counterpart. At O(p 4 ), there exist six independent invariants which obey charge conjugation invariance, and the effective Lagrangian at O(p 4 ) reads [18] where Due to parity, all potentials are odd functions of (ψ + θ), except forṼ 4 which is even.
In the combined LN c and chiral expansions, the WZW term starts contributing at O(N c p 4 ) = O(δ). Our aim is the calculation of the two-photon decays at the one-loop level, which corresponds to a NNLO calculation in the δ counting. Therefore, we need the odd-intrinsic-parity Lagrangians at NLO and NNLO. Up to and including NNLO, the effective odd-intrinsic-parity Lagrangian is denoted by where the superscripts (i) refer to the order in δ. The NLO Lagrangian L (2) receives contributions from O(p 4 ) and O(N c p 6 ). From Eq. (20) one can extract [18] The odd-intrinsic-parity Lagrangian at O(p 6 ) has been constructed in SU(3) ChPT in Refs. [23,24]. Reference [27] Table I. Since there is, at present, no satisfactory unified nomenclature for the coupling constants, for easier reference we choose the names according to the respective references from which the Lagrangians were taken.
The operators with the LECs L 6, i , which appear in SU(3) ChPT as well, are taken from Ref. [23]. They are given in terms of the building blocks where A refers to operators transforming under the chiral group G as A The other terms, genuinely related to the U(3) sector, are taken from Ref. [27]. Here, the corresponding building blocks are the same as in Eq. (19) with the additional structures Lagrangian LEC Operator SU(3)

D. Power counting
In the following, we provide the power-counting rules for a given Feynman diagram, which has been evaluated by using the interaction vertices derived from the effective Lagrangians of Eq. (2). Using the δ counting introduced in Eq. (1), we assign to any such diagram an order D which is obtained from the following ingredients: Meson propagators for both octet and singlet fields count as O(δ −1 ). Since meson fields are always divided by 2 ), a vertex with k meson fields derived from L (i) is O(δ i+k/2 ). The integration of a loop counts as δ 2 . The order D is obtained by adding up the contributions of the individual building blocks. The power-counting rules are summarized in Table II. External fields s and p 1  Recall that powers (ψ + θ) n come with expansion coefficients of O(N −n c ) even though we count (ψ + θ) as O(1).

III. CALCULATION OF THE INVARIANT AMPLITUDE
The invariant amplitude of the two-photon decay of a pseudoscalar meson P can be parameterized by where q µ 1 , q µ 2 denote the photon momenta, µ 1 , µ 2 the polarization vectors of the photons, and F P γ * γ * (q 2 1 , q 2 2 ) is the so-called transition form factor (TFF). In order to determine the invariant amplitude up to and including NNLO, we need to calculate the Feynman diagrams shown in Fig. 1. The vertices are derived from the Lagrangians given in Sec. II. The coupling of the electromagnetic field to the mesons is described by introducing an external field which couples to the electromagnetic current operator where Q is the quark-charge matrix. For N c = 3, the quark-charge matrix is given by However, as Bär and Wiese pointed out [28], in order for the Standard Model to be consistent for arbitrary N c , the ordinary quark-charge matrix should be replaced by The matrix element M is calculated using both versions, Q(3) and Q(N c ). In the Q(N c ) case, we first perform the δ expansion up to and including NNLO and then set N c = 3. The Feynman diagrams are evaluated using the Mathematica package FEYNCALC [29]. For the decays of η and η , we take into account the η-η mixing at NNLO 5 . A detailed derivation of the η-η mixing at NNLO can be found in Ref. [30]. First, we calculate the coupling of two photons to the octet and singlet fields φ b , collected in the doublet η A ≡ (η 8 , η 1 ) T , at the one-loop level up to and including NNLO in the δ counting. The result, which should be interpreted as a Feynman rule, is represented by the "matrix elements" F b = γ * γ * |b . In a next step, we transform the bare fields η A to the physical states using the transformation T in Eq. (51) in Ref. [30]: The resulting ("physical") matrix elements are then given by 5 No mixing with π 0 .
For the calculation of the loop diagrams, we employ the LO mixing.
Without the 1/N c expansion of Q, the results for the form factors of π 0 , η, η at LO and NLO read where θ [i] is the corresponding mixing angle at LO (NLO) obtained from Eq. (49) in Ref. [30]. The parameter λ 1 is a QCD-scale-invariant combination of parameters violating the Okubo-Zweig-Iizuka (OZI) rule, given by [18] Including the 1/N c expansion of Q, the results at LO and NLO now take the form At NNLO, the expressions for the form factors are quite long. Therefore, we do not display all terms explicitly. The loop contributions corresponding to the loop diagrams shown in Fig. 1 are provided in Appendix A. The expressions for the full NNLO form factors, with tree-level contributions, are available as Mathematica notebooks.

A. Observables
The decay amplitude for real photons is recovered by setting q 2 1 = q 2 2 = 0 in Eq. (26). The decay width is then given by [31] where λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2yz − 2xz is the Käellén function, λ 1 , λ 2 denote the polarizations of the photons, and dΩ is the solid angle of one of the photons. Using The single-virtual TFF F P γ * γ (q 2 ) := F P γ * γ * (q 2 , 0) can be measured in single Dalitz decays P → γl + l − . The slope of the TFF is defined as One can also define the dimensionless quantity b P = M 2 P × slope. The curvature is given by and the corresponding dimensionless quantity reads c P = M 4 P × curv. Experimental extractions of the slope parameter are often performed using a vectormeson-dominance model (VMD) [32] to fit the data. Introducing G P γ * γ (Q 2 ) = F P γ * γ (q 2 ) with Q 2 = −q 2 , in this case, the TFF is given by a normalized single-pole term with an associated mass Λ P [33]: Expanding this expression in Q 2 leads to Now, we can read off the slope and curvature VMD predictions, which are given by

IV. NUMERICAL ANALYSIS
We perform the numerical analysis of our results successively at LO, NLO, and NNLO. In the following, we distinguish between two cases: (a) using the normal quark-charge matrix for N c = 3 and (b) taking the 1/N c expansion of Q into account, denoted by Qexp. Performing the 1/N c expansion of Q shifts some of the LECs to higher orders. The LECs L 6, 8 and L 6, 19 , stemming from the NLO Lagrangian in Table I in Sec. II, appear only at NNLO in the expression for the two-photon decay of the π 0 . At LO, no unknown LECs show up and we can calculate the desired quantities directly.

Determinations of the parameters
First, we consider the Q(3) case. Since the decay width of the π 0 to two photons depends only on L 6, 8 , we start by fixing L 6, 8 to the experimental value of Γ π 0 →γγ . We then fit λ 1 simultaneously to the experimental results for Γ η→γγ and Γ η →γγ . The experimental values for the decay widths are taken from Ref.
[4] and are displayed in Table IV. Finally, the parameter L 6, 19 is determined through a simultaneous fit to the experimental values of the π 0 , η, and η slopes, given in Table IV. For the fits we employ the Mathematica routine NonlinearModelFit. The errors of the fit parameters and of the results are the ones obtained from the fit routine. The different steps are performed successively and we do not take the errors of LECs determined in a previous step into account. We also do not consider the errors due to neglecting higher-order terms. In principle, a systematic error of at least 10%, corresponding to δ 2 = 1/9, should be added to all quantities determined up to NLO. The results for the LECs are given in Table III and the results for the decay widths and slopes  in Table IV, labeled NLO 1.
Next, we examine the case where we fit L 6, 8 and λ 1 simultaneously to all three decay widths Γ π 0 →γγ , Γ η→γγ , and Γ η →γγ . The constant L 6, 19 is then again fixed to the slopes of π 0 , η, and η . The results are shown in Table III, labeled NLO 2. To consistently take the errors of L 6, 8 , L 6, 19 , and λ 1 into account, we consider another scenario where we determine these three LECs through a simultaneous fit to the decay widths of π 0 , η, η and the slope parameters of π 0 , η, η . The results are given in Table III, labeled NLO 3.
In the Q(N c ) case, the π 0 form factor at NLO is independent of L 6, 8 and L 6,19 , which are shifted to the NNLO expression. The width Γ π 0 →γγ takes the LO value of the Q(3) case, and the slope is equal to zero at NLO. Therefore, we determine L 6, 8 , L 6, 19 , and λ 1 via a simultaneous fit to Γ η→γγ , Γ η →γγ , b η , and b η . The results are displayed in Table III, labeled NLO, Qexp.
The parameters L 5 and θ [1] have been determined in the NLO analysis in Ref. [30] with a small error. Therefore, we have not taken these errors into account in the analysis of the two-photon decays. However, to obtain an estimate of the effect of the errors, we recalculated the quantities in the NLO 2 scenario varying L 5 and θ [1] within their errors. This led to only small variations in the last digit of the results for the decay widths and the slopes. We conclude that the influence of the errors of L 5 and θ 1 is very small, and we omit them in the following.
The LO values for the decay widths labeled LO agree within 20% with the experimental values. The slopes are equal to zero at that order. At LO, taking the 1/N c expansion of Q into account leads to results that are far from the experimental values. The NLO calculations improve the description of the decay widths. In the NLO 1 case, the π 0 decay width is equal to the experimental value, because L 6, 8 is fixed to it. In the NLO 2 case, where the parameters where fitted to all three decay widths, the description of Γ π 0 worsens, while Γ η and Γ η come closer to the experimental values. Our NLO 1-3 results for the slope of the η agree well with the experimental value. The description of the η slope, however, is very bad. Due to the small error of b η the fit favors this value, contributing to the poor description of b η . In the simultaneous fit to all decay widths and slope parameters (NLO 3), the results for the decay widths show larger deviations from the experimental values in comparison to NLO 2, marginally improving the values for the slopes. In the NLO, Qexp scenario, the π 0 decay width is given by the leading-order value of the Q(3) case. Since, then, the two parameters L 6, 8 and λ 1 need to be fixed by Γ η and Γ η alone, we reproduce the experimental values for these widths. The results for b η and b η are very poor in this case. In the NLO, Qexp case, the errors of the LECs and the results for the decay widths and slopes are very large. This further reflects the fact that the NLO, Qexp calculation is not appropriate to describe the data, and the LECs cannot be fixed in a sensible way. We thus conclude that omitting the 1/N c expansion of Q leads to a better description of the experimental data at LO and NLO. However, in general, the NLO calculation is not sufficient to adequately describe the decay widths and slopes of π 0 , η, and η , which motivates taking higher-order corrections into account.

Parametrization of the TFFs and determination of the parameters
At NNLO, a lot of new LECs appear both from the even-intrinsic-parity sector and the odd-intrinsic-parity sector. Moreover, our power counting demands taking terms of the O(p 8 ) Lagrangian into account, which has not been constructed yet. We therefore make the following ansatz for the q 2 dependence of the single-virtual TFFs up to and including NNLO: The A P and B P are combinations of LECs from the higher-order Lagrangians in Sec. II and, in principle, receive contributions from the O(p 8 ) Lagrangian as well. The C P stem solely from the O(p 8 ) Lagrangian. The expression loops P (q 2 ) denotes the q 2 -dependent part of the loop corrections, while the q 2 -independent parts are absorbed in the parameters A P . We determine the parameters A P , B P , C P through a simultaneous fit to the real-photon decay widths Γ P →γγ and to the experimental data for the TFFs. In the following, we perform several fits for the π 0 , the η, and the η TFF considering different NNLO contributions. We start by fitting the full NNLO expressions. Then, we consider the case without loops, which means switching off the q 2 -dependent pieces loops P (q 2 ). To study the influence of the C P terms, we also perform fits where we put C P = 0. Finally, we discuss the case where both C P and loops are neglected. In addition, we examine each of these four scenarios taking the 1/N c expansion of Q into account, denoted by Qexp. The fits are performed using the Mathematica routine NonlinearModelFit, and the errors of the fit parameters are the ones obtained by this routine. A systematic error of at least 4%, corresponding roughly to δ 3 = (1/3) 3 , should be added to all results determined up to NNLO. The π 0 TFF is fitted to both the time-like experimental data in Refs. [39][40][41][42] and the space-like data from Ref.
[1]. For each scenario, we fit the TFF to four different regions of the photon virtuality q 2 , which are given by −0.55 GeV 2 ≤ q 2 ≤ 0.55 GeV 2 , −0.5 GeV 2 ≤ q 2 ≤ 0.5 GeV 2 , −0.45 GeV 2 ≤ q 2 ≤ 0.45 GeV 2 , and −0.4 GeV 2 ≤ q 2 ≤ 0.4 GeV 2 . The results for the fit parameters obtained in the range −0.5 GeV 2 ≤ q 2 ≤ 0.5 GeV 2 are provided in Table V, while the results for the other fits are shown in Table XVI in Appendix B.
The TFF of the η is fitted to the time-like experimental data obtained in Refs. [42][43][44][45][46]. For each case, we perform fits up to three different values of the invariant mass of the lepton pair,   Table VI, and the results of the other fits can be found in Table XVII in Appendix B. For the η TFF there are also data points in the space-like low-energy region available. Therefore, we fit the TFF to the space-like data from Ref. [47] and to the time-like data from Ref. [48]. Here, we choose four fit regions for each scenario. The different fit ranges for the photon virtuality q 2 are −0.53 GeV 2 ≤ q 2 ≤ 0.43 GeV 2 (I), −0.53 GeV 2 ≤ q 2 ≤ 0.40 GeV 2 (II), −0.50 GeV 2 ≤ q 2 ≤ 0.43 GeV 2 (III), and −0.50 GeV 2 ≤ q 2 ≤ 0.40 GeV 2 (IV). Table  VII shows the results for the parameters fitted in the range −0.53 GeV 2 ≤ q 2 ≤ 0.43 GeV 2 , and the results of the other fits are displayed in Table XVIII in Appendix B.
The q 2 dependence of the π 0 TFF is displayed in Fig. 2 in the time-like region and in Fig. 3 in the space-like region. Here, the TFF is normalized to 1 at q 2 = 0 and plotted together with the experimental data. In this case, the TFF was fitted in the range −0.5 GeV 2 ≤ q 2 ≤ 0.5 GeV 2 . The bands show the 1σ error bands obtained by the Mathematica fit routine NonlinearModelFit.
In Figs. 4 and 5, the results for the different fit ranges are displayed in the time-like and in the space-like region, respectively.
The q 2 dependence of the normalized η TFF is shown in Fig. 6, where the TFF is plotted as a function of the invariant mass of the lepton pair m(l + l − ) together with the experimental data. In this case, the TFF was fitted up to 0.47 GeV. The bands show the 1σ error bands.   Figure 7 shows the results of the fits for the different fit ranges. As the fit range is extended to higher m(l + l − ) values, the curves become steeper. The q 2 dependence of the normalized η TFF, fitted between −0.53 GeV 2 and 0.43 GeV 2 , is shown in Fig. 8 together with the experimental data. The bands are the 1σ error bands due to the errors of the fit parameters. The results for the η TFF fitted to different ranges are displayed in Fig. 9.

Discussion of the results
In the following, we will interpret the results of the NNLO analysis. We start with the discussion of the LECs determined at NNLO, shown in Tables V-VII. Switching on or off the loop contributions corresponds to keeping the q 2 -dependent parts, loops P (q 2 ), or neglecting them, respectively. As a result, the parameters A P remain the same in both cases. The inclusion of the 1/N c expansion of Q has almost no visible effect on the shape of the TFFs. However, this expansion has an influence on the parameters A P , i.e., the absolute normalizations of the form factors, which change notably, since the LO expressions for the TFF (see Sec. III) are different with or without the 1/N c expansion of Q. The q 2 -dependent loop corrections loops P (q 2 ) give numerically quite similar contributions to the TFF with or without the 1/N c expansion of Q. Therefore, the parameters B P and C P do not vary very much in these two cases. As Figures 13 -16 in Appendix C show, the influence of the loop contributions on the shape of the TFFs is very small. However, the effects of the loops can be seen in the variation of B P and C P with and without loops, where we observe rather small changes. Neglecting the loops leads to an increase of the values of B P and C P in order to compensate for the missing contributions, which add positively to the TFFs.
Table V in Appendix B shows the variation of the fit parameters for the π 0 TFF with decreasing fit range. For C 0 π equal to zero, the parameter B 0 π increases only slightly if the fit range is decreased. In general, since the B 0 π term is not able to provide curvature, the curves are rather flat. If the C 0 π term is included, B 0 π and C 0 π decrease for the first three decreasing ranges, but increase for the last (smallest) range. When fewer data points are included, the curves become less steep, matching the decreasing parameter values. In the last scenario (IV), however, the fit is dominated by the single data point in the space-like region, leading to more curvature and therefore larger parameter values.
The results for η TFF fit parameters are displayed in Table XVII in Appendix B for decreasing fit ranges. If C η is put to zero, the parameter B η decreases as the fit range decreases. This behavior is in accordance with the fact that the curves become steeper as the fit range is extended to higher q 2 values. If we include the C η term in the fit, there is an interplay between C η and B η . For decreasing fit range, the C η values tend to decrease while B η increases. In addition, the errors of B η and C η become larger. This is to be expected, since less data points are included in the fit, there seems to be a correlation between B η and C η , and the C η term becomes more important at higher values of q 2 .
In the case of the η , the fit range is varied both in the time-like and the space-like region. The variation of the parameters is displayed in Table XVIII in Appendix B. Decreasing the time-like fit range yields smaller values for both B η and C η . This is to be expected, since the TFF curves show less curvature as the fit range gets smaller. If we exclude the last space-like data point, the values for B η and C η increase. In this case, the fit focuses more on the time-like region and the parameters adjust to the steep rise of the time-like TFF.

Slope and curvature
Employing the results for the fit parameters, we calculate the slopes and the curvatures of the TFFs as defined in Eqs. (47) and (48). The errors are due to the errors of the fit parameters. As a first estimate, we assume that the fit parameters are independent. Taking into account their correlations is beyond the scope of this work. The main results are given in Tables VIII-X. The values for the slopes with and without loop contributions agree within their uncertainties. This is the case, because the influence of the loops is already compensated by different values for the fit parameters B P . The 1/N c expansion of the quark-charge matrix plays a negligible role. For C π 0 = 0, the π 0 slope is smaller than in the full NNLO calculation. This is in accordance with the fact that the curves are less steep for C π 0 = 0, since the fit is dominated by the values at very low q 2 . If we neglect the C η term, B η compensates for the missing contribution, and, as a result, the η slope increases. This effect is diminished, if the fit range is restricted to lower q 2 values. In the case of the η , if we set C η = 0, the slope gets smaller. This behavior is different from the one of the η slope due to the inclusion of the space-like data. As a further check, we have investigated the case where the fit is performed only to the time-like data. Then, the η slope increases if we put C η = 0, which    is similar to the η case. Figures 10 -12 show the comparison of our results for the π 0 , η, and η slopes together with other experimental and theoretical determinations. Our values b η = 0. 55 fit up to 0.4 GeV, and b η = 1.47(31) from fit I agree within the errors with most of the other theoretical and experimental results. The result b π 0 = 0.025(3) from the fit up to 0.5 GeV 2 , however, is smaller than the other theoretical predictions and most of the experimental results. This is due to the inclusion of higher-q 2 data. Our result for the fit only up to 0.5 GeV, b π 0 = 0.028(5), e.g., is closer to the other predictions. In general, our result for b η is slightly lower than the other determinations, whereas our value for b η is slightly higher than the other results. From Table IX one can observe that a decreasing fit range leads to values for b η which come closer to the other theoretical and experimental determinations. The main results for the curvatures of π 0 , η, and η are displayed in Tables VIII -X. As the slope, the π 0 curvature decreases with decreasing fit range except for the smallest range, where the curvature takes its largest value. In this case, the TFF is adjusted to  the single data point in the space-like region leading to a large curvature. The η curvature is reduced if the fit range is restricted to smaller q 2 values. The η curvature decreases if the time-like range is decreased. As the space-like fit range becomes smaller, the fit is dominated by the steeply rising time-like data and the curvature is almost twice as large.
The main contributions to the curvature stem from the C η terms. If we put them to zero, the remaining curvature is given by the loop contributions, which is rather small. Our values for the curvatures can be compared with other theoretical determinations. A dispersive analysis finds c π 0 = 1.14(4) × 10 −3 [59], and works using Padé approximants obtain c π 0 = 1.06(26) × 10 −3 [56], c η = 0.339(15) stat (5) sys [38], and c η = 1.72(47) stat (34) sys [33]. If we use a simple VMD estimate as given in Eq. (52) with Λ P = 0.77 GeV [31], we obtain c 0 π = 0.9×10 −3 , c η = 0.26, and c η = 2.40. Our results for c π 0 are slightly smaller than the other predictions, being closest to the VMD one. For c η , our values are mostly larger than the other predictions. Only if the fit range is decreased, our results come to agreement  with Ref. [38], whereas the naive VMD prediction is even smaller. In the cases including the full space-like data, the η curvature is slightly smaller than the one from Ref. [33], but shows agreement within the errors. The VMD prediction for c η is larger and lies on the upper end of the error band in Ref. [33]. None of our values reaches the VMD value within the error range. Note that the errors are only the ones provided by the fit. The results for c η in the cases where the space-like fit range is restricted are much larger than the ones from the fit to all space-like data as well as the ones from the other references.

V. SINGLE DALITZ DECAYS
Having performed the numerical evaluation of the single-virtual TFFs of π 0 , η, and η , we are now able to calculate the widths of the decays to one photon and a lepton pair. To obtain the invariant amplitude for the decay P → γl + l − , where P = π 0 , η, η and l = e, µ, we use Eq. (26) and therein define q µ 1 , with q 2 1 = s, and µ 1 = (e/s)[ūγ µ v] as virtual-photon momentum and polarization, respectively. The momentum of the real photon is denoted by q 2 with q 2 2 = 0, and 2 is its polarization. The decay width can be written as [31] Γ Defining the leptonic tensor and employing the identity one obtains To evaluate this expression numerically, we make use of the LECs determined in Secs. IV A and IV B. At NNLO, we employ the LECs determined from the fit of the π 0 TFF up to 0.5 GeV 2 , the η TFF up to 0.47 GeV, and the η TFF in the range −0.53 GeV 2 ≤ q 2 ≤ 0.43 GeV 2 (I). The results for the decay widths to one photon and a lepton pair are shown in Tables XI and XII. The errors are calculated from the errors of the LECs which are assumed to be uncorrelated. They can be viewed as upper limits for the errors. Taking the correlations into account is beyond the scope of this work. The values for Γ P →γe + e − , P = π 0 , η, η , behave like the corresponding values for the decays to two real photons. The disagreement of the two-photon-decay widths in some scenarios and the experimental data is reflected in the values for Γ P →γe + e − as well. Therefore, we calculate the relative branching ratios (BR) using the values for Γ P →γγ obtained in the different scenarios. The results are shown in Tables XI and XII. Now, the values for the relative BRs do not vary very much within the different cases and orders. The π 0 relative BRs for the decay to an e + e − pair agree with the experimental value and also the η relative BR is very close to the experimental one, while the η relative BR is somewhat smaller than the experimental result, especially in most of the NLO cases. This is related to the value of the η slope. The slope is very large in the NLO 1 case, which leads to a large relative BR, and the negative values for b η in the other NLO cases are reflected in a reduced relative BR even when compared to the LO value. The decay width of P → γe + e − receives its main contribution at values where the virtual photon is in the vicinity of its mass shell. Therefore, the relative BRs are well described already at LO. The decay P → γµ + µ − provides a better probe of the virtual behavior of the TFF at larger photon virtualities. However, the values for Γ η→γµ + µ − are still related to the two-photon-decay widths, but the higher-order corrections in q 2 , parameterized by slope and curvature, become important. Here, we calculate the relative BRs as well. The LO relative BR of the η is now lower than the experimental value and increases at NLO and NNLO. Especially at NNLO, we obtain a very good agreement with the data for both Γ η→γµ + µ − and BR rel η→γµ + µ − . The LO relative BR for η → γµ + µ − is only 30% of the experimental value. In the NLO scenarios, it becomes even smaller except for NLO 1. This is related to the slope of the η , which is very large in the NLO 1 scenario, but poorly described in the other NLO cases, even cases with negative values. The full NNLO value is larger than the LO one and most of the NLO values. However, it is still smaller than the experimental result. If we neglect the C η term, the relative BR decreases again. This is connected to the description of the η TFF data. The time-like TFF is underestimated for higher values of q 2 and even more so if one does not take the (q 2 ) 2 term into account. The decay width of η → γµ + µ − receives contributions in q 2 ranges where vector-meson resonances become important [77], which are not included in our framework.

BR rel
η→γe + e − [10 −2 ] QED [72] 1.18 ± 0 Hidden gauge [73] 1.19 ± 0 Mod. VMD [73]  Our full NNLO results for the relative BRs are compared with other theoretical determinations. The comparison for the π 0 , η, and η relative BRs can be found in Tables XIII -XV. Our values for BR rel P →γe + e − , P = π 0 , η, η , agree well with the other determinations. In the case of BR rel η→γµ + µ − , as already stated, the simple QED prediction is too small. Here, our result agrees with the other works, except for Ref. [74] which gives a slightly smaller value.  Our result for BR rel η →γµ + µ − is smaller than the others. It agrees within errors with Ref. [77], and the determinations including vector mesons are larger.

VI. SUMMARY AND OUTLOOK
We have studied the P → γ ( * ) γ ( * ) interaction, where P = π 0 , η, η , at the one-loop level up to and including NNLO in LN c ChPT. Besides the loop corrections, all contact terms appearing at NNLO have been calculated, except for those of the O(p 8 ) Lagrangian, which has not been constructed yet. However, in the expressions for the form factors describing the decays, possible structures originating from the O(p 8 ) Lagrangian have been introduced phenomenologically, accompanied by free parameters. Furthermore, the η-η mixing at NNLO has been consistently included. The numerical analyses of the decays have been performed successively at LO, NLO, and NNLO. At NLO, we employed the values for the LECs and mixing angle determined in the NLO analysis of the η-η mixing in Ref. [30] labeled NLO I. The LECs from the odd-intrinsic-parity sector were fixed to the experimental data of the decay widths to real photons and the slope parameters of π 0 , η, η . We have found that the NLO results are not sufficient to describe all data simultaneously. If the 1/N c expansion of the quark-charge matrix is taken into account, the results worsen. At NNLO, the LECs have been determined through a fit to the experimental data for the π 0 , η, and η transition form factors. We have achieved a good description of the π 0 TFF between −0.45 GeV 2 and 0.45 GeV 2 , of the η TFF up to 0.45 GeV, and of the η TFF between −0.25 GeV 2 and 0.3 GeV 2 , which is mainly caused by the inclusion of (q 2 ) 2 terms, whereas loops do not play an important role. In addition, we have calculated the slopes and the curvatures of the TFFs and the decay widths of P → γl + l − , where l = e, µ, and compared them to other works. In general, our NNLO results for those quantities tend to agree with the other experimental and theoretical determinations.