Electromagnetic decays of the neutral pion

TL;DR

This paper calculates electromagnetic decay widths of the neutral pion using a non-perturbative QCD approach, providing results consistent with previous theories and experiments, and explores a numerical method for complex integrals.

Contribution

It introduces a symmetry-preserving Dyson-Schwinger equation approach to compute pion decay widths and employs a novel numerical contour deformation method for rare decay calculations.

Findings

Results agree with previous theoretical and experimental data.

The rare decay calculation aligns with dispersive methods but shows a 2 sigma discrepancy with experiment.

The numerical method can be generalized to other integrals with known analytic structure.

Abstract

We complement studies of the neutral pion transition form factor pi^0 --> gamma^(*) gamma^(*) with calculations for the electromagnetic decay widths of the processes pi^0 --> e^+ e^-, pi^0 --> e^+ e^- gamma and pi^0 --> e^+ e^- e^+ e^-. Their common feature is that the singly- or doubly-virtual transition form factor serves as a vital input that is tested in the non-perturbative low-momentum region of QCD. We determine this form factor from a well-established and symmetry-preserving truncation of the Dyson-Schwinger equations. Our results for the three- and four-body decays match results of previous theoretical calculations and experimental measurements. For the rare decay we employ a numerical method to calculate the process directly by deforming integration contours, which in principle can be generalized to arbitrary integrals as long as the analytic structure of the integrands are known. Our result for the rare decay is in agreement with dispersive calculations but still leaves a 2 sigma discrepancy between theory and experiment.

Tables

Table 1. Table 1: Result for the Dalitz decay.

Collaboration	${\Gamma }_{{\pi }^0\to e^{+}{e}^{-}\gamma }[{10}^{-11}$ GeV]
PDG Patrignani et al. (2016)	$9.06\left(18\right)$
Terschlüsen et al. Terschlüsen et al. (2013)	$9.26$
Hoferichter et al. Hoferichter et al. (2014)	$9.065$
Our result	$9.11\left(4\right)$

Table 2. Table 2: Result for the double Dalitz decay.

Collaboration	${\Gamma }_{{\pi }^0\to 2e^{+}2e^{-}}[{10}^{-13}$GeV]
PDG Patrignani et al. (2016)	$2.58\left(12\right)$
Terschlüsen et al. Terschlüsen et al. (2013)	$2.68$
Escribano et al. Escribano and Gonzàlez-Solís (2015)	$2.62$
Our result	$2.63\left(1\right)$

Table 3. Table 3: Our result for the rare decay, obtained either with the dispersion relation (DR) or directly from the contour deformation, compared to other theoretical calculations and experiment (after removing the final state radiative corrections).

Collaboration	B(${\pi }^0\to e^{+}{e}^{-}\left)\right[{10}^{-8}]$
Experiment Abouzaid et al. (2007); Vasko and Novotny (2011); Husek et al. (2014)	$6.87\left(36\right)$
Dorokhov et al. Dorokhov and Ivanov (2007); Dorokhov (2010)	$6.23\left(9\right)$
Husek et al. Husek and Leupold (2015)(THS)	$6.14\left(8\right)$
Masjuan et al. Masjuan and Sanchez-Puertas (2015)	$6.23\left(5\right)$
Our result (DR)	$6.21\left(3\right)$
Our result (direct)	$6.22\left(3\right)$

Equations

$${\Lambda }^{\mu \nu }(Q,Q^{\prime })=e^2\frac{F(Q^2,{Q^{\prime }}^2)}{4{\pi }^2{f}_{\pi }}{\epsilon }^{\mu \nu \alpha \beta }{Q^{\prime }}^{\alpha }{Q}^{\beta },$$

$$\Sigma =\frac{Q+Q^{\prime }}{2},\Delta =Q-Q^{\prime },$$

$${\eta }_{+}$$

$$\omega$$

$${\eta }_{-}=Q\cdot Q^{\prime }={\Sigma }^2-\frac{{\Delta }^2}{4}={\eta }_{+}+\frac{m_{\pi }^2}{2}.$$

$$\sigma ={\Sigma }^2,t=\frac{{\Delta }^2}{4},Z=\hat{\Sigma }\cdot \hat{\Delta }=\frac{\Sigma \cdot \Delta }{2\sqrt{\sigma t}}$$

$${\Lambda }^{\mu \nu }$$

$$\times {\Gamma }^{\mu }(r_{-},-Q)S(k+\Sigma ){\Gamma }^{\nu }(r_{+},Q^{\prime }).$$

$$k_{\pm }=k\pm \frac{\Delta }{2},r_{\pm }=k+\frac{\Sigma }{2}\pm \frac{\Delta }{4},$$

$$S(p)=Z_f(p^2)\frac{-i/p+M(p^2)}{p^2+M^2(p^2)}$$

$${\Gamma }_{\pi }(k,\Delta )=\left(f_1+f_{2}i/\Delta +f_{3}k\cdot \Delta i/k+f_4\left[/k,/\Delta \right]\right){\gamma }_5,$$

$${\Gamma }^{\mu }(p,Q)=\sum _{i=1}^{12}{\lambda }_i(p^2,p\cdot Q,Q^2){\tau }_i^{\mu }(p,Q)$$

$$F(Q^2,{Q^{\prime }}^2)=\frac{\mathcal{A}(w)+w(1-z^2){\mathcal{B}}_1(w)(1+{\mathcal{B}}_2(w)z^2)}{(1+w{)}^2-w^2{z}^2}.$$

$$\begin{array}{rl}\mathcal{A}(w)& =\frac{a_0+\xi (a_1{b}_{1}w+a_2{b}_2{w}^2+a_3{b}_3{w}^3)}{1+b_{1}w+b_2{w}^2+b_3{w}^3},\\ {\mathcal{B}}_i(w)& =\frac{c_i{e}_i{w}^2}{1+d_{i}w+e_i{w}^2}\end{array}$$

$$\begin{array}[]{rl}a_{1}&=0.735\,,\\ a_{2}&=1.214\,,\\ a_{3}&=1.547\,,\\[2.84526pt] c_{1}&=0.384\,,\\ d_{1}&=2.010\,,\\ e_{1}&=1.540\,,\end{array}\qquad\begin{array}[]{rl}b_{1}&=0.089\,,\\ b_{2}&=0.133\,,\\ b_{3}&=0.0002\,,\\[2.84526pt] c_{2}&=0.430\,,\\ d_{2}&=0.024\,,\\ e_{2}&=0.00005\,.\end{array}$$

$${\Gamma }_{{\pi }^0\to e^{+}{e}^{-}\gamma }=\frac{e^6{m}_{\pi }^3}{6(4\pi {)}^3}{\int }_{4m^2}^{m_{\pi }^2}\frac{dx}{x}{\left|\frac{F(Q^2,0)}{4{\pi }^2{f}_{\pi }}\right|}^2$$

$$\times \sqrt{1-\frac{4m^2}{x}}\left(1+\frac{2m^2}{x}\right){\left(1-\frac{x}{m_{\pi }^2}\right)}^3,$$

$$F(Q^2,Q^{\prime 2})=0.996-3.55(10){\eta }_{+},$$

$${\Gamma }_{{\pi }^0\to 2e^{+}2e^{-}}$$

$${\left|\mathcal{M}\right|}^2={\left|{\mathcal{D}}_1\right|}^2+{\left|{\mathcal{D}}_2\right|}^2+2\mathrm{Re}[{\mathcal{D}}_1{\mathcal{D}}_2^{\ast }].$$

$${\Gamma }_{{\pi }^0\to 2e^{+}2e^{-}}={\Gamma }_{{\pi }^0\to 2e^{+}2e^{-}}^{(\text{direct})}+{\Gamma }_{{\pi }^0\to 2e^{+}2e^{-}}^{(\text{indirect})},$$

$$|{\mathcal{D}}_1{|}^2+|{\mathcal{D}}_2{|}^2=2e^4\frac{{\Lambda }^{\mu \nu }(Q,Q^{\prime }){\Lambda }^{\alpha \beta }(Q,Q^{\prime })}{Q^4{Q^{\prime }}^4}\text{Tr}{G}_{34}^{\mu \alpha }\text{Tr}{G}_{12}^{\nu \beta },$$

$${\Gamma }_{{\pi }^0\to 2e^{+}2e^{-}}^{(\text{direct})}=\frac{e^8}{36(2\pi {)}^5{m}_{\pi }^3}\underset{4m^2}{\overset{(m_{\pi }-2m{)}^2}{\int }}dx\underset{4m^2}{\overset{(m_{\pi }-\sqrt{x}{)}^2}{\int }}dy$$

$$\times \sqrt{x-4m^2}\sqrt{y-4m^2}{\left[\frac{(x+y-m_{\pi }^2{)}^2}{4xy}-1\right]}^{3/2}$$

$$\times \left(1+\frac{2m^2}{x}\right)\left(1+\frac{2m^2}{y}\right){\left|\frac{F(Q^2,{Q^{\prime }}^2)}{4{\pi }^2{f}_{\pi }}\right|}^2.$$

$$\mathrm{Re}[{\mathcal{D}}_1{\mathcal{D}}_2^{\ast }]=e^4\frac{{\Lambda }^{\mu \nu }(Q,Q^{\prime }){\Lambda }^{\alpha \beta }(K,K^{\prime })}{Q^2{Q^{\prime }}^2{K}^2{K^{\prime }}^2}\text{Tr}{G}_{32}^{\mu \beta }{G}_{14}^{\nu \alpha },$$

$${\Gamma }_{{\pi }^0\to 2e^{+}2e^{-}}^{(\text{direct})}$$

$${\Gamma }_{{\pi }^0\to 2e^{+}2e^{-}}^{(\text{indirect})}$$

$$R=\frac{B({\pi }^0\to e^{+}{e}^{-})}{B({\pi }^0\to \gamma \gamma )}=2{\left(\frac{m{\alpha }_{\text{em}}}{\pi m_{\pi }}\right)}^2\beta (t_0)|\mathcal{A}(t_0){|}^2,$$

$$\int \frac{d^4\Sigma }{(2\pi {)}^4}\Lambda (p_{+}){\gamma }^{\mu }S(p+\Sigma ){\gamma }^{\nu }\Lambda (p_{-})\frac{{\Lambda }^{\mu \nu }(Q,Q^{\prime })}{Q^2{Q^{\prime }}^2}$$

$$=\frac{\mathcal{A}(t)}{(4\pi {)}^2}\frac{2im{\alpha }_{\text{em}}}{\pi f_{\pi }}\Lambda (p_{+}){\gamma }_5\Lambda (p_{-}),$$

$$\mathcal{A}(t)$$

$$\text{Im}{\mathcal{A}}^{\text{LO}}(t)=\frac{\pi \ln \gamma (t)}{2\beta (t)}F(0,0),$$

$$R\ge {\left(\frac{m{\alpha }_{\text{em}}}{m_{\pi }}\right)}^2\frac{{\ln }^2\gamma (t_0)}{2\beta (t_0)}=4.75\times 10^{-8}.$$

Full text

I Introduction

Over the last years, low-energy electromagnetic processes in the meson sector have seen continuous interest both in the theoretical and experimental physics communities. Electromagnetic decays of the pion are particularly interesting because they combine the non-perturbative physics of dynamical mass generation and the associated generation of (pseudo-) Goldstone bosons with the Abelian anomaly and its perturbative elements, thus creating an interesting laboratory for theoretical approaches to non-perturbative QCD. Moreover, the rare decay $\pi^{0}\to e^{+}e^{-}$ poses a puzzle since theoretical estimates show a discrepancy with the experimental result from the KTeV E799-II experiment at Fermilab Abouzaid et al. (2007); Dorokhov and Ivanov (2007) of similar magnitude as the muon g-2.

In this work we focus on electromagnetic decays of pseudoscalar mesons such as the rare decay $\pi^{0}\to e^{+}e^{-}$, the Dalitz decay $\pi^{0}\to e^{+}e^{-}\gamma$ and double Dalitz decay $\pi^{0}\to e^{+}e^{-}e^{+}e^{-}$ using a well explored combination of Dyson-Schwinger and Bethe-Salpeter equations. All calculations involve the $\pi^{0}\to\gamma^{(\ast)}\gamma^{(\ast)}$ transition form factor (TFF) with one or two off-shell photons, thus testing it in the region of (very) low momenta. The present work complements a recent evaluation of the TFF at large space-like momenta, see Ref Eichmann et al. (2017).

The paper is organised as follows. In the next section we give a short introduction to the details of our calculations and discuss features of the resulting TFF. In section III we then give results for the leptonic three- and four-body decays of the neutral pion and compare with the experimental values. In section IV we discuss corresponding results for the rare decay of the pion into an electron-positron pair. We conclude in section V.

II Transition form factor $\pi^{0}\to\gamma^{\ast}\gamma^{\ast}$

II.1 Kinematics and definitions

We begin by defining the transition current and the kinematic regions of interest. The $\pi^{0}\to\gamma\gamma$ transition matrix element is given by

[TABLE]

where $Q$ and $Q^{\prime}$ are the photon momenta, $f_{\pi}\approx 92$ MeV is the pion’s electroweak decay constant, and $e^{2}=4\pi\alpha_{\text{em}}$ the squared electromagnetic charge. The pseudoscalar transition is described by a single scalar function, the transition form factor $F(Q^{2},{Q^{\prime}}^{2})$, and the convention of prefactors is such that $F(0,0)=1$ in the chiral limit due to the Abelian anomaly.

In the following it is useful to work with the average photon momentum $\Sigma$ and the pion momentum $\Delta$,

[TABLE]

with $\Delta^{2}=-m_{\pi}^{2}$ for an on-shell pion. The process then depends on two Lorentz invariants,

[TABLE]

or vice versa: $\{Q^{2},\,{Q^{\prime}}^{2}\}=\eta_{+}\pm\omega$, with the third invariant fixed when the pion is on-shell:

[TABLE]

Note that the TFF is symmetric in $Q^{2}$ and ${Q^{\prime}}^{2}$ so it can only depend on $\omega$ quadratically.

For practical calculations it is convenient to introduce the alternative variables

[TABLE]

where $t=-m_{\pi}^{2}/4$ when the pion is on-shell and a hat denotes a normalized four-vector. We will refer to those in section IV and Appendix A.

In the physical processes we study in this work the TFF is tested in both space-like and time-like regions as shown in Fig. 1. The time-like region, where either $Q^{2}<0$ or ${Q^{\prime}}^{2}<0$, contains the physical singularities such as the vector-meson poles in the complex plane of $Q^{2}$ and ${Q^{\prime}}^{2}$. For dilepton decays this region is kinematically restricted such that $Q^{2},Q^{\prime 2}\geq-m_{\pi}^{2}$. The double Dalitz decay $\pi\to 2e^{+}2e^{-}$ is constrained to the light blue shaded area below $\eta_{+}<0$, whilst the Dalitz decay $\pi\to e^{+}e^{-}\gamma$ probes the asymmetric time-like form factor, given by the dark blue lines along the $Q^{2},Q^{\prime 2}$ axes. These decays are discussed in section III.

The space-like region with both photon virtualities positive, $Q^{2}>0$ and $Q^{\prime 2}>0$, is free of any physical singularities. The region that is relevant for the rare decay $\pi^{0}\to e^{+}e^{-}$, discussed in section IV, is the doubly-virtual or symmetric limit (red line in Fig. 1) when $Q^{2}=Q^{\prime 2}$ viz. $\omega=0$. Direct experimental measurements of the spacelike TFF are available in the singly-virtual or asymmetric limit with one of $\left\{Q^{2},Q^{\prime 2}\right\}$ vanishing Behrend et al. (1991); Gronberg et al. (1998); Aubert et al. (2009); Uehara et al. (2012).

In addition, very different kinematic regions to these can be tested where the pion is ‘off-shell’ corresponding to space-like momentum transfer $\Delta^{2}>0$. This and various applications are discussed in Eichmann et al. (2017).

II.2 Triangle diagram

In the impulse approximation, the transition form factor $\pi^{0}\rightarrow\gamma^{(*)}\gamma^{(*)}$ is displayed in Fig. 2 and given by

[TABLE]

The photon momenta were defined in the previous subsection. In addition, $k$ is the loop momentum and

[TABLE]

are the internal quark momenta and the relative momenta appearing in the quark-photon vertices, respectively. The trace in Eq. (II.2) is over Dirac indices only111All quantities are color-diagonal and thus the color trace is 3. The flavor matrix of the $\pi^{0}$ amplitude is $\mathrm{diag}\,(1,-1)$ and that of the quark-photon vertex is $\mathrm{diag}\,(q_{u},q_{d})$, so the flavor trace is $q_{u}^{2}-q_{d}^{2}=1/3$. Since the overall normalization of the pion amplitude is determined by the canonical Bethe-Salpeter norm, which follows from demanding unit residue at the pion pole in the $q\bar{q}$ scattering matrix, a different color-flavor normalization can always be absorbed by the dressing functions in Eq. (10). Our choice above is such that $f_{1}(k^{2})=B(k^{2})/f_{\pi}$ in the chiral limit, with $B(k^{2})=M(k^{2})/Z_{f}(k^{2})$. and the factor 2 in front of the integral stems from the exchange of the photons.

All ingredients in Eq. (II.2) are determined from numerical solutions of their Dyson-Schwinger and Bethe-Salpeter equations. The renormalized quark propagator is given by

[TABLE]

with non-perturbative dressing functions $Z_{f}(p^{2})$ and $M(p^{2})$. The renormalization-group invariant quark mass function $M(p^{2})$ encodes effects of dynamical mass generation due to the dynamical breaking of chiral symmetry. The Dirac structure of the pseudoscalar Bethe-Salpeter amplitude $\Gamma_{\pi}$ is given by

[TABLE]

where the $f_{i}$ are functions of $k^{2}$ and $k\cdot\Delta$ with $\Delta^{2}=-m_{\pi}^{2}$ fixed.

The non-perturbative quark-photon vertex $\Gamma^{\mu}$ describes the coupling of a dressed quark to a photon and is dominated by QCD corrections. It can be decomposed into twelve tensor structures,

[TABLE]

with basis components $\tau_{i}^{\mu}(p,Q)$ and Lorentz-invariant dressing functions $\lambda_{i}$; see App. B of Ref. Eichmann et al. (2016) for details. The argument $p$ denotes the average momentum of the two quark legs and $Q$ is the incoming photon momentum. Due to electromagnetic gauge invariance the vertex can be split into a transverse and non-transverse part, where the latter is the Ball-Chiu vertex Ball and Chiu (1980) and determined by the vector Ward-Takahashi identity. We obtain a numerical solution of the quark-photon vertex from its inhomogeneous Bethe-Salpeter equation; see Maris and Tandy (2000a, b, 2002); Bhagwat and Maris (2008); Goecke et al. (2011); Eichmann (2011). As discussed in detail in Maris and Tandy (2000a); Goecke et al. (2013); Eichmann et al. (2016), the transverse part of the vertex contains poles in the time-like momentum region corresponding to vector-meson states. Thus, the underlying physics of vector-meson dominance is automatically contained in the numerical representation of the vertex without the need for further adjustments.

In this work we restrict ourselves to the rainbow-ladder approach as reviewed in Eichmann et al. (2016). We use the Maris-Tandy model, Eq. (10) of Ref. Maris and Tandy (1999) with parameters $\Lambda=0.74$ GeV and $\eta=1.85\pm 0.2$ (the parameters $\omega$ and $D$ therein are related to the above via $\omega D=\Lambda^{3}$ and $\omega=\Lambda/\eta$). The variation of $\eta$ changes the shape of the quark-gluon interaction in the infrared, cf. Fig. 3.13 in Ref. Eichmann et al. (2016), and we use it in the following to estimate our theoretical error. This construction, with the same respective kernel in the Bethe-Salpeter equation for the pseudoscalar mesons and the one for the quark-photon vertex, preserves chiral symmetry and has the merit of producing reliable results in the pseudoscalar and vector-meson sector as well as for nucleon and $\Delta$ baryons. Our input current-quark mass is $m_{q}=3.57$ MeV at a renormalization point $\mu=19$ GeV; the resulting pion mass and pion decay constant are $m_{\pi^{0}}=135.0(2)$ MeV and $f_{\pi^{0}}=92.4(2)$ MeV.

II.3 Result for the transition form factor

With the ingredients described above we are able to determine the $\pi^{0}\to\gamma\gamma$ transition form factor in the spacelike domain $Q^{2}>0$ and ${Q^{\prime}}^{2}>0$ as well as for small timelike momenta. In practice it turns out that a straightforward calculation is only possible in restricted kinematic regions. As explained in Appendix A, this is due to the singularities of the quark propagator in the integrand, whose nearest singularities correspond to a scale $m_{p}\sim 0.5$ GeV. The symmetric limit is accessible for all $\eta_{+}>0$, whereas in the asymmetric limit one is limited to $Q^{2}_{\text{max}}\approx 4$ GeV2, which is also the domain covered in the calculation of Ref. Maris and Tandy (2002). In addition, also small timelike momenta of the order $Q^{2}$, ${Q^{\prime}}^{2}\gtrsim-m_{p}^{2}$ are accessible directly; cf. Fig. 9 in the appendix.

To determine the TFF in the full spacelike domain, we employ the strategy introduced in Ref. Eichmann et al. (2017): we calculate the form factor for an off-shell pion with $\Delta^{2}>0$ and extrapolate to the on-shell point using the lowest-lying vector-meson pole as a constraint. This allows us to determine the TFF for all space-like momenta, and in the regions that are accessible by direct calculation and extrapolation we have checked that the results of both methods are in perfect agreement.

In Fig. 3 we show the on-shell transition form factor $F(Q^{2},{Q^{\prime}}^{2})$ as a function of the variable $\eta_{+}$, Eq. (3). The plot reveals that the TFF is essentially a function of $\eta_{+}$, which is larger in the asymmetric limit and smaller in the symmetric limit, and this behavior persists up to asymptotically large $\eta_{+}$ Eichmann et al. (2017). Moreover, in the chiral limit the Abelian anomaly entails $F(0,0)=1$ and our numerical result at the physical pion mass is $F(0,0)=0.996$, which provides an important consistency check: replacing both dressed quark-photon vertices by bare ones would only give $F(0,0)\approx 0.29$; and even a Ball-Chiu vertex, which guarantees charge conservation in the pion’s electromagnetic form factor, produces $F(0,0)\approx 0.86$ only. The transverse structure of the vertex is therefore crucial for a quantitative description of the $\pi^{0}\to\gamma\gamma$ transition.

We finally provide a fit function that accurately represents our results in the spacelike region. Abbreviating $w=\eta_{+}/m_{v}^{2}$ and $z=\omega/\eta_{+}$, the TFF is described by

[TABLE]

The denominator implements the lowest-lying vector-meson pole at $wz=\pm(1+w)$, which corresponds to $Q^{2}=-m_{v}^{2}$ and ${Q^{\prime}}^{2}=-m_{v}^{2}$ with $m_{v}=0.77$ GeV. The functions in the numerator ensure that the TFF asymptotically approaches a monopole behavior both in the symmetric ($z=0$) and asymmetric limit ($z=\pm 1$); they are given by

[TABLE]

with fit parameters $a_{0}=0.996$ and

[TABLE]

This fit provides the input for our calculations of the various $\pi^{0}$ decays. The value $\xi=1.0\pm 0.1$ reflects our combined theoretical uncertainty from varying the parameter $\eta=1.85\pm 0.2$ in the effective interaction as well as the uncertainty in the determination of the TFF away from the symmetric limit. These are also the error estimates that we quote in the following results.

Let us finally briefly comment on alternative fit functions available in the literature (see, e.g., Appendix B of Ref. Nyffeler (2016) for a discussion). The simplest vector-meson dominance (VMD) parametrization $F_{\text{VMD}}(Q^{2},{Q^{\prime}}^{2})=1/[(1+w)^{2}-w^{2}z^{2}]$ does not reproduce the monopole behavior in the symmetric limit $z=0$ but instead approaches a dipole at large $Q^{2}$. Its refined version based on lowest-meson dominance, the LMD+V model Knecht and Nyffeler (2001), reproduces both the symmetric and the asymmetric limits and implements two vector-meson poles; an analogous form was recently employed to fit lattice results for the TFF Gerardin et al. (2016). Our fit is practically indistinguishable from the LMD+V parametrization at low $Q^{2}$, i.e., in the momentum range shown in Fig. 3. Also at large $Q^{2}$ the fits in the symmetric limit are almost identical and in the asymmetric limit they are at least qualitatively similar. However, the behavior in between ($0<|z|<1$) differs substantially: at large $\eta_{+}$ the LMD+V form factor develops a sharp peak very close to $|z|=1$, with a turnover to reach the asymmetric point $z=1$ followed by the vector-meson poles at $|z|>1$. By comparison, our fit varies monotonously from $z=0$ to $z=1$ and is therefore better suited for applications where the TFF is tested in the whole spacelike domain.

III Three- and four-body decays

In this subsection we discuss our results for the three- and four-body decays of pseudoscalar mesons shown in Fig. 4(a-b). The Dalitz decay of the neutral pion into a photon and an electron-positron pair has the largest branching ratio $B(\pi^{0}\to e^{+}e^{-}\gamma)=(1.174\pm 0.035)\%$ Patrignani et al. (2016) after that of the two-photon decay. The neutral pion also decays into two dilepton pairs with a branching ratio of $B(\pi^{0}\to e^{+}e^{-}e^{+}e^{-})=(3.34\pm 0.16)\times 10^{-5}$ Patrignani et al. (2016). Both decays depend on the transition form factor discussed above as the only non-trivial input.

III.1 Dalitz Decay: $\pi^{0}\to e^{+}e^{-}\gamma$

The leading-order Feynman diagram for the three-body decay of the neutral pion is shown in Fig. 4(a). The decay rate is easily calculated and given by

[TABLE]

where $m$ is the electron mass, $m_{\pi}$ the pion mass, and $x=-Q^{2}$ is the momentum squared of the virtual photon which evaluates the form factor in the timelike region.

Due to the kinematics of the three-body decay the TFF is probed in the asymmetric region of one photon on-shell and one off-shell and time-like. The leading contribution to the integral in Eq. (III.1) comes from the lower end of the integral. In this region the calculated TFF can in general be described by a simple linear fit with respect to $\eta_{+}=(Q^{2}+Q^{\prime 2})/2$,

[TABLE]

which is also used in the four-body decay below. The TFF is shown in Fig. 5, where the asymmetric limit required for the three-body decay, $F(Q^{2},0)$, is represented by either one of the dark blue lines.

Employing the PDG values for $m$ and $\alpha_{\text{em}}$ in Eq. (III.1), our calculated decay width is $\Gamma_{\pi\to e^{+}e^{-}\gamma}=9.11(4)\times 10^{-11}$ GeV. Since the theoretical uncertainty in Eq. (16) affects this number only at the sub-per-mille level, the error bar in the decay rate comes from the model dependence of $m_{\pi}$ and $f_{\pi}$ which enter in Eq. (III.1). In Table 1 we compare our result with the PDG and the result of a calculation using an effective theory Terschlüsen et al. (2013). Within the quoted errors all results are in good agreement. This is to be expected because the TFF – as the only non-perturbative input – is probed in the kinematic region governed by the anomaly. Since any reasonable construction obeys this constraint, the discriminative potential of the Dalitz decay with regard to different non-perturbative input is very limited.

III.2 Four-body decay: $\pi^{0}\rightarrow e^{+}e^{-}e^{+}e^{-}$

We now proceed with the four-body decay of the neutral pion into two electron-positron pairs. The decay rate is given by

[TABLE]

where $\left|\mathcal{M}\right|^{2}$ is the squared and spin-summed matrix element and the symmetry factor in front accounts for the two pairs of identical final-state particles. $d\Phi_{4}$ is the four-dimensional phase space measure whose detailed derivation is given in Appendix B. Because the amplitude $\mathcal{M}$ with all initial and final particles on-shell depends on five independent variables, $d\Phi_{4}$ involves five non-trivial integrations.

Fig. 4(b) shows the possible diagrams, $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$, where exchange of two leptons (anti-leptons) introduces a relative minus sign between the two contributions. Squaring these we obtain

[TABLE]

As discussed in Refs. Terschlüsen et al. (2013); Escribano and Gonzàlez-Solís (2015), the first two terms are equal and can thus be combined. It follows that the decay rate can be decomposed into

[TABLE]

where the first (direct) contribution comes from the two squared magnitudes and the second (indirect or interference) term comes from the cross-terms.

Abbreviating $G^{\mu\nu}_{ij}=(i\not{p}_{i}+m)\,\gamma^{\mu}\,(i\not{p}_{j}-m)\,\gamma^{\nu}$, the integrand for the direct contribution is given by

[TABLE]

where $\Lambda^{\mu\nu}(Q,Q^{\prime})$ is the $\pi^{0}\to\gamma\gamma$ transition current defined in Eq. (1). Because in this case the integrand only depends on the pairwise sums $Q^{\prime}=-(p_{1}+p_{2})$ and $Q=p_{3}+p_{4}$, the five-dimensional phase-space integral $d\Phi_{4}$ can be reduced to just two integrations, namely

[TABLE]

Here we used the shorthands $x=-Q^{2}$ and $y=-{Q^{\prime}}^{2}$ for the two photon virtualities, where $x,y>0$ are timelike and restricted by the thresholds for the two two-body decays. The direct contribution constitutes the largest fraction of the decay rate.

In contrast, the interference term depends on all possible four-vector combinations of electron and positron pairs. Consequently, the phase-space integral cannot be further reduced and its integrand is given by

[TABLE]

where $Q^{\prime}=-(p_{1}+p_{2})$, $Q=p_{3}+p_{4}$ and $K^{\prime}=-(p_{2}+p_{3})$, $K=p_{1}+p_{4}$ are the possible momenta of the virtual photons. To perform the five-dimensional integral we used various methods, ranging from tensor-product quadrature with a combination of Gauss-Legendre and double exponential rules, to 5-dimensional adaptive cubature as well as standard Monte-Carlo methods Hahn (2005); all agreed perfectly.

For the direct and indirect contributions to the decay rate we obtain

[TABLE]

The sum of these values gives our final result, shown in Table 2 together with the value from experiment and other theoretical calculations. In Ref. Terschlüsen et al. (2013) the authors use an extension of chiral perturbation theory to calculate the form factor, whereas in Ref. Escribano and Gonzàlez-Solís (2015) a data driven approach is employed, which is based on the use of rational approximants applied to the experimental data of the $\pi^{0},\eta$ and $\eta^{\prime}$ transition form factors in the space-like region. All results are compatible with experiment within the quoted error. As with the Dalitz decay, the phase-space restriction of the virtual photons to time-like momenta between $4m^{2}$ and $m_{\pi}^{2}$ entails that the sensitivity to the details of the TFF beyond that dictated by the anomaly is rather small and deviates from its nominal value of $1$ by no more than $3$%.

IV Rare decay: $\pi^{0}\rightarrow e^{+}e^{-}$

Finally we consider the two-body decay of the neutral pion into one electron-positron pair. For the $\pi^{0}$ this is certainly the most interesting decay due to a discrepancy between the KTeV experimental result and theoretical calculations Abouzaid et al. (2007); Dorokhov and Ivanov (2007); Dorokhov (2010); Vasko and Novotny (2011); Husek et al. (2014); Masjuan and Sanchez-Puertas (2015) of the order of $2\sigma$. Using the elaborate reanalysis of radiative corrections Vasko and Novotny (2011); Husek et al. (2014) to the experimental result of the KTeV collaboration Abouzaid et al. (2007) (close to the value given in PDG Patrignani et al. (2016)) one arrives at an extracted experimental value for the branching ratio of $B(\pi^{0}\rightarrow e^{+}e^{-})=(6.87\pm 0.36)\times 10^{-8}$, which is considerably smaller than the decays considered above.

To lowest order in QED the process is described by the one-loop graph in Fig. 4(c), which again includes the transition form factor $F(Q^{2},{Q^{\prime}}^{2})$ as the only non-perturbative input. Defining $t=\Delta^{2}/4$ as in Eq. (6), the corresponding normalized branching ratio is given by

[TABLE]

where $\beta(t)=\sqrt{1+m^{2}/t}$ stems from the two-body phase-space integration and $B(\pi^{0}\to\gamma\gamma)=0.988$. The scalar amplitude $\mathcal{A}(t)$ can be viewed as the pseudoscalar form factor of the electron due to the two-photon coupling, which must be evaluated at the on-shell pion point $t_{0}=-m_{\pi}^{2}/4$.

IV.1 $\mathcal{A}(t)$ with dispersive input

For arbitrary $t$ the amplitude $\mathcal{A}(t)$ can be defined from the matrix element for the $\pi^{0}\to e^{+}e^{-}$ decay:

[TABLE]

where $\Lambda^{\mu\nu}(Q,Q^{\prime})$ is the $\pi^{0}\to\gamma\gamma$ transition current from Eq. (1) and $\Lambda(p_{\pm})=\frac{1}{2}\,(\mathds{1}+\not{p}_{\pm}/(im))$ is the positive-energy projector of the lepton. The kinematics are as discussed in Sec. II.1; in particular, the averaged photon momentum $\Sigma$ becomes the loop momentum and therefore the variables in Eq. (6) take the values $\sigma>0$ and $Z\in[-1,1]$. As a consequence, the photon virtualities $Q^{2}$ and ${Q^{\prime}}^{2}$ are tested at complex values close to the symmetric limit as shown in Fig. 6. In principle the integral depends on the pion momentum $\Delta$ and the average lepton momentum $p$, but since the electron and the positron are on-shell with momenta $p_{\pm}^{2}=(p\pm\Delta/2)^{2}=-m^{2}$ only $t$ remains as an independent variable.

Taking traces yields the following expression for $\mathcal{A}(t)$:

[TABLE]

This integral has poles in the integration domain (which we discuss in more detail in Sec. IV.2) and thus cannot be naively integrated except for the unphysical point $t=\Delta^{2}/4=0$. A standard way to circumvent the problem uses dispersive methods (see e.g.  Bergstrom et al. (9839); Donoghue (1996) ). In that case the imaginary part of the amplitude along its cut at $t<0$ is given by

[TABLE]

with $\gamma(t)=(1-\beta(t))/(1+\beta(t))$, which follows from cutting the two photon lines. The imaginary part gives the well-known unitary bound for the branching ratio through the inequality $|\mathcal{A}(t_{0})|^{2}\geq|\text{Im}\mathcal{A}(t_{0})|^{2}$:

[TABLE]

Using a once-subtracted dispersion relation one then obtains the real part of the amplitude via

[TABLE]

where $\text{Li}_{2}(z)$ is the dilogarithm or Spence function. In particular, this implies $\text{Re}\ \mathcal{A}(t_{0})=\mathcal{A}(0)+31.92(2)$ so that the only unknown left is the constant $\mathcal{A}(0)$.

In fact, $t=0$ is the only point where Eq. (26) can be integrated directly to yield

[TABLE]

where we temporarily abbreviated $x=Q^{2}/(4m^{2})$. A similar formula can be derived using a Mellin-Barnes representation Ghaderi (2013); Dorokhov and Ivanov (2008); Dorokhov et al. (2009),

[TABLE]

which is however only valid to leading order in an expansion in the electron mass. At the point $t=0$ the transition form factor in both cases is evaluated in the symmetric limit of equal photon momenta, and due to $Q^{2}=4m^{2}x$ it is mainly probed at very low $Q^{2}$ of the order of the electron mass. Implementing our result for $F(Q^{2},{Q^{\prime}}^{2})$, we extract the same value $\mathcal{A}(0)=-21.85(2)$ from both formulas above, where the error comes from varying the $\xi$ parameter in Eq. (13). With Eqs. (27–28) one then arrives at the on-shell value $\mathcal{A}(t_{0})=10.07(4)-17.45(1)\,i$, which corresponds to a branching ratio of

[TABLE]

Our result is compared to other approaches in Table 3. Whereas our calculation represents a top-down approach using a well-tested model for the underlying quark-gluon interaction, Refs. Dorokhov and Ivanov (2007); Dorokhov (2010) use a phenomenological parametrization of the transition form factor that is adapted to reproduce experimental data from CLEO together with additional high-energy QCD constraints. A generalization of LMD+V is the two-hadron saturation model (THS) of Ref. Husek and Leupold (2015). The more recent Ref. Masjuan and Sanchez-Puertas (2015) employs a data-driven approach via Padé Theory and Canterbury approximants. All four theoretical results are in agreement with each other, thus showing consistency between different approaches. Again, it appears that the decay rate is not overly sensitive to different representations of the form factor as long as the QCD constraints are satisfied (as guaranteed in all three approaches). However, we would like to point out that all three calculations rely on dispersion relations and the Mellin-Barnes representation. Thus the only number that influences the final result is the constant $\mathcal{A}(0)$. Although a priori one would deem the dispersive approach reliable for this process, it still remains to be checked via a direct calculation.

IV.2 Direct calculation of $\mathcal{A}(t)$

The integrand in Eq. (26) has poles for vanishing denominators, i.e., if either of the photons or the intermediate lepton go on-shell. Depending on the value of $t$, this may prohibit a straightforward Euclidean integration. Specifically, for $t\in\mathds{C}$ one can draw a kinematically safe region in the complex $t$ plane where such an integration is possible, and a forbidden region where the poles enter in the integration domain and the integration would produce a wrong result. The latter case would usually be interpreted as a failure of the Wick rotation; however, as we demonstrate below, the Euclidean expression Eq. (26) is still valid if the poles are treated correctly. Problems of this kind are frequent in Euclidean bound-state calculations and pose limitations, e.g., in computing excited hadrons or form factors for time-like or large space-like arguments Eichmann et al. (2016) and thus it is desirable to find a general method to deal with them.

In the case of Eq. (26) it is the unfortunate combination of all three external momenta being on-shell that complicates the situation. The analysis in Appendix A shows that the lepton poles lead to a narrow parabola

[TABLE]

around the negative (time-like) $t$ axis which is kinematically safe, whereas the photon poles admit a straightforward integration only for real and positive $t$. Taken in combination, the integration is only possible for $t=0$, which leads to the result for $\mathcal{A}(0)$ quoted above.

To analyze the situation for general $t$, we write the integral in hyperspherical variables defined by

[TABLE]

cf. Eq. (6), where the process becomes particularly simple in the frame $\Delta=2\sqrt{t}\,[0,0,0,1]$, $p=i\sqrt{t+m^{2}}\,[0,0,1,0]$ and

[TABLE]

The innermost $\psi$ integration is trivial and thus the integral $\mathcal{A}(t)$ takes the form

[TABLE]

The denominator under the $dZ$ integral is equal to $Q^{2}\,{Q^{\prime}}^{2}$ whereas the second denominator corresponds to the lepton pole. The integration in $Y$ can be performed analytically:

[TABLE]

for all $a\in\mathds{C}$ except $-1<a<1$, in which case the logarithm develops a branch cut.

After performing all angular integrations, the conditions $Q^{2}\,{Q^{\prime}}^{2}=0$ and $-1<a<1$ produce poles and cuts in the complex $\sigma$ plane which are visualized in Fig. 7. In principle, the $\sigma$ integration goes from zero to infinity but the ‘naive’ Euclidean’ integration path $\sigma\in\mathds{R}_{+}$ would cross a singularity, hence causing the problems described above. The vanishing denominator for the photons produces a cut along

[TABLE]

with $-1<Z<1$, which describes a circle with radius $|t|$. The circle opens at $\sigma=t$ since in this case the remaining $(1-Z^{2})$ factor in the the denominator of Eq. (IV.2) cancels with the numerator. On the other hand, the lepton denominator leads to a cut

[TABLE]

which does not open at $\sigma=t$ as this would correspond to $a=0$ in Eq. (39). Instead, it never crosses the arc that passes through $\sigma=-t$ as shown in Fig. 7.

We solve the problem by exploiting Cauchy’s theorem in finding an integration contour that connects $\sigma=0$ and $\sigma\to\infty$ but never crosses any singularity. Such a possible path is shown in Fig. 7: (1) it departs from the origin in the direction opposite to $t$, (2) then navigates between the cuts until it reaches the point $\sigma=t$, (3) returns to the positive real axis at a value $\sigma>|t|$, (4) and from there proceeds to the numerical cutoff at $\sigma\to\infty$. This strategy ensures that for each $\sigma$ along the path the integrand in $Z$ is free of any singularities.

In practice, however, the problem is made worse by the small lepton mass $m\ll m_{\pi}$: in that case both cuts essentially describe the same circle but open on opposite sides, $\sigma=t$ and $\sigma=-t$, so that the path along segment (2) proceeds along a narrow ridge between the two cuts. As a consequence, the integrand in $Z$ is sharply peaked which requires an enormous number of grid points to obtain stable results. To this end we have optimized the procedure so that path (2) is always equally spaced between the cuts, thus minimizing the numerical error. In addition, we implement an adaptive integration in $-1<Z<1$ that accumulates the grid points according to the nearest singularities in the complex $Z$ plane, which are determined by solving $\sigma_{\gamma}^{\pm}(Z)$ and $\sigma_{l}^{\pm}(Z)$ for $Z$. In that way we gain a factor of $\sim 10^{3}$ in CPU time while maintaining the same numerical accuracy.

As a result, we are in principle able to determine $\mathcal{A}(t)$ in the whole complex plane, which is shown in Fig. 8 in the momentum region relevant for the $\pi^{0}\to e^{+}e^{-}$ decay. The real part is sharply peaked at the origin $t=0$, whereas the imaginary part develops the expected branch cut on the timelike axis. The resulting on-shell value is $\mathcal{A}(t_{0})=10.10(3)-17.45(1)\,i$, where the error reflects the uncertainty in the TFF discussed below Eq. (14). The corresponding branching ratio

[TABLE]

nicely agrees with the result in Eq. (31) obtained via dispersion relations.

We have thus established a fast and efficient numerical method to calculate Euclidean integrals, which is applicable even in cases where a naive Wick rotation fails. It can be generalized to arbitrary integrals as long as the singularity structure of the integrand is known, which is not restricted to real poles but accommodates complex poles or cuts as well. Perhaps unsurprisingly, this demonstrates that both Euclidean and Minkowski descriptions are completely equivalent as long as the singularities in the integrand are treated correctly: the correct Euclidean expression, such as Eq. (IV.2), is the one that connects zero with infinity without crossing any singularities.

Exploratory calculations using contour deformations have been performed, e.g., in determining quark, gluon and ghost propagators in the complex plane Maris (1995); Alkofer et al. (2004); Eichmann (2009); Strauss et al. (2012); Windisch et al. (2013). In those cases the singularity structure of the integrands is usually less intertwined. In turn, one has to deal with integral equations and thus self-consistent problems where the singularities that are dynamically generated by the integration enter again into the integrand and must be accounted for as well. In addition to other methods to calculate propagators in the complex plane Fischer et al. (2005, 2009); Krassnigg (2008); Rojas et al. (2014), or a direct inclusion of residues as done for example in Ref. Oettel et al. (2000), contour deformation methods may become an attractive tool for overcoming several long-standing obstacles in the Dyson-Schwinger/Bethe-Salpeter approach, e.g., the calculation of (highly) excited states, form factors at timelike or large spacelike momentum transfer, or the treatment of genuine resonances.

V Conclusions

In this work we determined the branching ratios of various leptonic decays of the neutral pion. The common central element in these calculations is the pion transition form factor. We calculated its momentum dependence in a Dyson-Schwinger and Bethe-Salpeter framework employing an underlying quark-gluon interaction that has been successful elsewhere in describing a range of different static and dynamic meson and baryon properties, see e.g. Maris and Roberts (2003); Eichmann et al. (2016) for reviews. The resulting form factor dynamically generates vector meson dominance and can be shown to comply with the perturbative scaling behavior at large space-like momenta in the symmetric limit, whereas the coefficient of the scaling limit is modified in the asymmetric limit Eichmann et al. (2017).

Due to the kinematic constraints imposed by Dalitz decays – as is well-known – only a very small range of time-like momenta around the zero-momentum limit controlled by the anomaly is probed. Consequently, our results for these pion decays very much resemble those of other theoretical approaches that use, e.g., vector meson dominance models for the form factor supplemented with constraints from the anomaly. The results are subsequently compatible with experiment.

For the rare decay the form factor is probed in the space-like region and we have been able to confirm theoretical calculations using dispersion relations, however through a direct calculation using contour deformations. Thus our result for the rare decay $\pi^{0}\rightarrow e^{+}e^{-}$ still leaves a 2$\sigma$ discrepancy between theory and experiment.

VI Acknowledgements

We thank R. Alkofer, S. Leupold, A. Szczepaniak and M. T. Hansen for enlightening discussions. This work was supported by the DFG collaborative research centre TR 16, the BMBF grant 05H15RGKBA, the DFG Project No. FI 970/11-1, the FCT Investigator Grant IF/00898/2015, the GSI Helmholtzzentrum fuer Schwerionenforschung, and by the Helmholtz International Center for FAIR.

Appendix A Singularity restrictions

Here we return to the problem mentioned in Section II.3, namely, to find the kinematic regions where a Euclidean integral with singularities in the integrand can be calculated directly without any contour deformations. For generality, consider a Lorentz-invariant integral

[TABLE]

It depends on a collection of external momenta $p_{i}^{\mu}$ and one integrates over the loop momenta $k_{j}^{\mu}$. The integrand consists of Lorentz-invariant functions such as $f(q^{2})$, where $q^{\mu}$ is a linear combination of the external and loop momenta. The loop momenta $k_{j}^{\mu}$ are always real; however, if an external momentum is time-like ($p_{i}^{2}<0$) it will inject imaginary components into $q^{\mu}$, which is therefore a complex four-vector. Splitting $q^{\mu}$ into its real and imaginary parts, where $A$, $B\in\mathds{R}$ and $e_{\mu}$, $e^{\prime}_{\mu}$ are Euclidean unit vectors, yields:

[TABLE]

Because $e\cdot e^{\prime}\in[-1,1]$, the Lorentz invariant $q^{2}$ is tested inside a parabola $(A\pm iB)^{2}$ with apex $-B^{2}$ on the real axis; for $B=0$ it becomes the positive real axis. In other words, if $q^{\mu}$ has imaginary components then $q^{2}$ is sampled inside a complex parabola.

Now, if $f(q^{2})$ has singularities in the complex plane, the corresponding parabola passing through the first singularity (e.g., a real pole, complex conjugate poles, or the onset of a cut) defines the ‘contour mass’ $q^{2}=-m_{p}^{2}$. The kinematically safe region is then subject to the restriction

[TABLE]

Because the imaginary part of $q^{\mu}$ can only come from the external momenta $p_{i}^{\mu}$, this imposes restrictions on the domain of the external Lorentz invariants.

We specifically consider the transition matrix element in Eq. (II.2), where the quark momenta in the loop are given by $k\pm\Delta/2$ and $k+\Sigma$. The dressing functions of the quark propagator in Eq. (9),

[TABLE]

develop a certain singularity structure in the complex plane. In a rainbow-ladder truncation the nearest singularities are complex conjugate poles with a typical contour mass $m_{p}\sim 0.5$ GeV for light quarks (see Windisch (2016) for a detailed investigation). The loop momentum $k$ is always real and thus one arrives at the conditions

[TABLE]

In general these restrictions depend on the chosen frame for $\Sigma$ and $\Delta$. The components of $\Sigma$ and $\Delta$ can be arranged arbitrarily as long as the Lorentz invariants defined in Eqs. (3–6), $\Sigma^{2}=\sigma=\eta_{+}+m_{\pi}^{2}/4$, $\Delta^{2}=-m_{\pi}^{2}$ and $\Sigma\cdot\Delta=\omega$, remain unchanged. For any possible choice, however, one can find a linear combination $\Sigma+\alpha\,\Delta$ that has a four-component only, with an arbitrary parameter $\alpha\in\mathds{R}$. The general arrangement satisfying these constraints is

[TABLE]

with $\mathcal{N}=\sqrt{\sigma+2\alpha\omega-\alpha^{2}m_{\pi}^{2}}$.

Take for example the pion’s rest frame, which corresponds to $\alpha\to\infty$:

[TABLE]

If $\sigma>0$ and $Z\in[-1,1]$, $\Sigma$ is real and the imaginary part only comes from $\Delta$, so the resulting condition $m_{\pi}<2m_{p}$ is always satisfied. However, the situation is different in the quadrant $\eta_{+}>0$, $|\omega|<\eta_{+}$ shown in Fig. 1, because in this case $Z$ is imaginary. Here both $\Delta$ and $\Sigma$ have imaginary four-components and the resulting condition becomes $\left|\omega\right|<m_{\pi}m_{p}$, which defines a narrow strip around the symmetric limit $\omega=0$.

The arbitrariness of $\alpha$ can be exploited to optimize the frame, i.e. to reach kinematic regions for the form factor that are not accessible in the pion rest frame. The resulting domains are plotted in Fig. 9 for different values of $\alpha$. The leftmost plot shows $\alpha=0$ and the rightmost plot $\alpha=8$; for $\alpha\to\infty$ one recovers the pion rest frame. For example, with $\alpha=0$ the momentum $\Sigma$ is the one with a four-component only, whereas $\Delta$ has an imaginary three-component and, for $\sigma>0$, leads to the condition

[TABLE]

For the singly-virtual form factor $F(Q^{2},0)$ the optimal choice is $\alpha=1/2$ (second plot in Fig. 9). In that case it is the photon momentum $Q=\Sigma+\Delta/2$ that has a four-component (resembling the Breit frame in elastic form factor calculations). The denominator $\mathcal{N}=\sqrt{\eta_{+}+\omega}$ is always real if $\eta_{+}>0$ and $|\omega|<\eta_{+}$, and therefore the three-components of $\Sigma$ and $\Delta$ are imaginary. The resulting condition is

[TABLE]

This region crosses the line $\eta_{+}=\omega$ at

[TABLE]

with $\varepsilon=m_{\pi}^{2}/(4m_{p}^{2})$. Hence, in the asymmetric limit the form factor can be calculated up to $Q^{2}_{\text{max}}\approx 16\,m_{p}^{2}\approx 4$ GeV2 without crossing any quark singularities.

As a second example we consider the integral $\mathcal{A}(t)$ in Eq. (26) which describes the $\pi^{0}\to e^{+}e^{-}$ decay. In that case the external momenta are $\Delta$ and $p$, with Lorentz invariants $\Delta^{2}=4t$, $p^{2}=-(t+m^{2})$ and $p\cdot\Delta=0$, whereas $\Sigma$ is the real loop momentum. The internal photon momenta are $\Sigma\pm\Delta/2$ and the lepton momentum is $\Sigma+p$, so the singularity conditions become

[TABLE]

where $m$ is the lepton mass. The analogous arrangement in the general frame is

[TABLE]

with $\mathcal{N}=\sqrt{-(t+m^{2})+4\alpha^{2}t}$. In this case the maximal domains correspond to the limits $\alpha=0$ or $\alpha\to\infty$ and one arrives at the condition $(\text{Im}\,t)^{2}<4m^{2}\,\text{Re}\,(-t)$ from the lepton pole and $t>0$ from the photon poles, which are quoted in Eq. (32). Although in principle these regions are disjoint, the analysis of Eq. (IV.2) shows that the point $t=0$ contains an integrable singularity and thus $\mathcal{A}(0)$ is well-defined if $m>0$.

Appendix B Four-body phase space

In this appendix we work out the four-body phase-space integral $d\Phi_{4}$ that enters in the $\pi^{0}\to e^{+}e^{-}e^{+}e^{-}$ decay of Eq. (17). The decay width of a particle with momentum $P$ and mass $M$ decaying into $n$ daughter particles with momenta $p_{i}$ and masses $m_{i}$ is given by

[TABLE]

with $S$ the symmetry factor, $|\mathcal{M}|^{2}$ the spin summed and squared matrix element, and $d\Phi_{n}$ the phase-space integral for an $n$-particle final state given by

[TABLE]

Following Ref. Byckling and Kajantie (1971), we rewrite the integration in terms of invariant mass variables. For $n=4$ one obtains

[TABLE]

where, for degenerate decay products with $m_{i}=m$, the two- and three-particle Mandelstam variables read

[TABLE]

and the four-dimensional Gram determinant $\Delta_{(4)}$ contains all possible dot products of four-vectors:

[TABLE]

In contrast to Escribano and Gonzàlez-Solís (2015); Byckling and Kajantie (1971); Nyborg et al. (1965) we employ a Euclidean signature, however with $s_{ij}$ and $s_{ijk}$ defined such that they have the same meaning in Minkowski and Euclidean conventions. To work with the invariant mass variables of Eq. (50), one replaces the $p_{i}\cdot p_{j}$ in the Gram determinant according to Eq. (51) together with

[TABLE]

The physical region of integration is bounded by the surface $\Delta_{(4)}=0$. Following the derivation of Refs. Escribano and Gonzàlez-Solís (2015); Byers and Yang (1964); Nyborg et al. (1965), we impose this relation on the invariant mass variables and begin by solving $\Delta_{(4)}=0$ for $s_{14}$. It yields

[TABLE]

where $\lambda(u,v,w)=u^{2}+v^{2}+w^{2}-2uv-2uw-2vw$ is the Källen function, the $G$ functions are given by

[TABLE]

and $b$ is the coefficient of $-8\Delta_{(4)}$ linear in $s_{14}$:

[TABLE]

with $c=s_{124}-m^{2}$ and $d=s_{134}-m^{2}$.

The regions of the $s_{124}$ and $s_{134}$ integrations are bounded by the surfaces satisfying $s_{14}^{+}=s_{14}^{-}$, which is fulfilled when either of the $G$ functions in Eq. (54) vanishes $G(s_{124},s_{34},s_{12})=0$ or $G(s_{134},s_{12},s_{34})=0$. Solving $G(u,v,w)=0$ for $u$ yields

[TABLE]

and thus determines the integration boundaries $s_{124}^{\pm}$ and $s_{134}^{\pm}$ as functions of the two dilepton invariant masses $s_{12}$ and $s_{34}$.

Finally, $s_{34}$ and $s_{12}$ range from the threshold at $4m^{2}$ to $(M-\sqrt{s_{12}})^{2}$ and $(M-2m)^{2}$, respectively; here the ordering is arbitrary and could also be exchanged.

A valuable check when rewriting the phase space integral in different coordinates is the massless limit. For massless daughter particles the $n$-body phase space is given by

[TABLE]

Integrating over the phase space volume only, as given in Eq. (50) with the borders as suggested for $m=0$, reproduces the limit exactly as it should.

The final decay rate for the decay of the neutral pion ($M=m_{\pi}$) into two dileptons ($m=m_{e}$) is then given by

[TABLE]