Precise determination of the branching ratio of the neutral-pion Dalitz decay

We provide a new value for the ratio $R={\Gamma(\pi^0\to e^+e^-\gamma(\gamma))}/{\Gamma(\pi^0\to\gamma\gamma)}=11.978(6)\times10^{-3}$, which is by two orders of magnitude more precise than the current Particle Data Group average. It is obtained using the complete set of the next-to-leading-order radiative corrections in the QED sector, and incorporates up-to-date values of the $\pi^0$-transition-form-factor slope. The ratio $R$ translates into the branching ratios of the two main $\pi^0$ decay modes: $\mathcal{B}(\pi^0\to\gamma\gamma)=98.8131(6)\,\%$ and $\mathcal{B}(\pi^0\to e^+e^-\gamma(\gamma))=1.1836(6)\,\%$.


I. INTRODUCTION
The present 3% experimental precision on B(π 0 → e + e − γ) [1] represents a limitation for rare π 0 decay measurements, which commonly use the π 0 Dalitz decay for normalization. In particular, it dominates the uncertainty on B(π 0 → e + e − e + e − ) [2] and is the largest source of uncertainty on B(π 0 → e + e − ) [3]. For the same reason, the precision on B(π 0 → e + e − γ) is becoming a limiting factor for rare kaon decay measurements. An example is the K + → π + e + e − decay [4]: accurate knowledge of B(π 0 → e + e − γ) would improve the precision on the rate measurement by 30 %, and the precision on the low-energy parameter a + [5] by 10 %. The uncertainty on B(π 0 → e + e − γ) also dominates the precision on the K ± → π ± π 0 e + e − rate measurement [6], and is among the principal contributions to the uncertainties on the measured K L,S → π + π − e + e − rates [7]. In these circumstances, considering the improving precision on rare decay measurements, and the recent progress on the π 0form-factor measurement [8] and radiative corrections to the π 0 Dalitz decay [9], both a precision measurement of B(π 0 → e + e − γ) and an updated theoretical evaluation of this branching ratio are becoming more important.
Besides the direct extraction of the ratio R from the experiment, the shape of the singly-virtual π 0 transition form factor can be measured. This can be expanded in the transfered momentum squared, with the linear coefficient called the (form-factor) slope a π . Since the slope a π embodies the most relevant input to the ratio R regarding the (non-perturbative) low-energy QCD sector, its knowledge from the experiment is crucial to obtain a model-independent prediction of R. Recently, it was measured in the NA62 experiment, which analyzed 1.1×10 6 fully reconstructed π 0 Dalitz decays with the result a π = 0.0368(57) [8], taking into account the complete set of NLO radiative corrections in the QED sector [9]. The current PDG value is dominated by two inputs: the above NA62 result and the value provided by the CELLO collaboration [14] by extrapolation from the space-like region.
Our calculation combines a wide range of theoretical models and available experimental results on the formfactor shape and the well-established QED calculation including the complete set of NLO corrections, taking into account higher orders in the QED expansion in an conservative uncertainty estimate. As such, it represents a precise and reliable improvement (by two orders of magnitude) to the current PDG-based value of R, which might be further used in various theoretical predictions and experimental analyses. Moreover, for the first time, the slope corrections were not neglected in the bremsstrahlung contribution. Finally, we present the ratio R for the full as well as partial kinematic regions.
Measurements of the form-factor shape or the ratio R require significant theoretical input and crucially depend on a proper incorporation of radiative corrections. Consequently, a statement that the experiment itself provides a more relevant value of R than our theoretical prediction is by its nature imprecise. However, the computation would not be possible without the experimental evidence that the form-factor slope lies within a certain range of values. An example of how theoretical inputs influence the experimental values in this sector is the well-known discrepancy in the rare decay π 0 → e + e − driven most probably by the approximate radiative corrections [15] which do not agree with the exact calculation [16]; for details and discussion see [17,18].

II. THEORETICAL FRAMEWORK
Considering the QED expansion, the leading-order differential decay width for the π 0 Dalitz decay reads where the two-photon decay width is parametrized as Above, M π is the neutral-pion mass and m e is the electron mass. The definition (3) holds to all orders in the QED and Chiral Perturbation Theory expansions and covers also possible physics from beyond the SM, putting these effects into the form-factor normalization F (0). As usual, we have introduced kinematical variables x and y defined as The slope a π and curvature b π of the singly-virtual pion transition form factor are defined in terms of the Taylor expansion in the invariant mass of the vector current [19]: With particular theoretical models at hand, one can explore the properties of the singly-virtual pion transition form factor F (q 2 ) and calculate a π and b π . As examples we are going to discuss the vector-meson dominance (VMD) ansatz, the lowest-meson dominance (LMD) model [20] and the two-hadron saturation (THS) model  Table I. Slope and curvature of the singly-virtual pion transition form factor evaluated in various approaches. The THS model [18] is numerically consistent with VMD. On the other hand, the value predicted by LMD [20] is unsatisfying. Regarding the values for the slope aπ, none of the models (VMD, LMD, THS) fully (within 1 σ) agrees with the results of the dispersive calculation [23], the data driven method of Padé approximants [24], the PDG average [1] or the recent measurement performed by the NA62 experiment [8], which though agree with each other well.
proposed in [18]. These models belong to the family of large-N c motivated analytic resonance-saturation models and as such they can be used in a straightforward way in the calculation of radiative corrections. Let us mention that the simplest models do not satisfy some of the major high-energy constraints: VMD violates operator-product-expansion constraints [18,21] and the LMD model cannot satisfy the Brodsky-Lepage scaling limit [22]. Finally, the THS model by construction satisfies all the mentioned constraints. By means of the first and second derivatives it is straightforward to find the analytic expressions for a π and b π within these form-factor models; for details see Section 5 of [18]. Numerical results are shown in Table I, together with other theoretical approaches and experimental results. From (5) it follows that In [8], for the estimation of the slope the prescription f 2 (x) NA62 ≡ (1 + a π x) 2 was used, which is accurate up to order O(x), but introduces an error of 2b π x 2 . If we consider b π ≃ a 2 π , as suggested by the models (cf .  Table I), we introduce an error of σ(a π ) = a 2 π + O(a 3 π ), i.e. a relative error η(a π ) ≃ a π ≃ 3 %, on the estimate of a π . To avoid such incompatibility and an introduction of an additional uncertainty, we will use in the following the same prescription as was used in [8] to evaluate R, i.e. f 2 (x) NA62 . Finally, we can use the expansion (6) to obtain the formula for the integrated decay width. Inserting (6) into (2) and taking into account that x ∈ (4m 2 e /M 2 π , 1) and Our aim is also to address the NLO radiative corrections in the QED sector. It is convenient to introduce the NLO correction δ to the LO differential decay width (and thus to write schematically dΓ = (1+δ +. . .) dΓ LO ), which can be in general (for the two-fold and one-fold differential case, respectively) defined as (8) One can obtain δ(x) from δ(x, y) using the following prescription: To calculate the NLO radiative corrections we use the approach documented in [9,25], which reviewed and extended the classical work of [26]. Hence, together with the bremsstrahlung (BS) beyond the soft-photon approximation, we take into account in the following calculations the one-photon-irreducible (1γIR) contribution. Let us discuss the case when the 1γIR contribution to the NLO radiative corrections is not considered in the analysis to extract the form-factor slope from the data. If we start with the equation among the measured spectral shapes (one-fold differential decay widths) and eliminate dΓ LO /dx| f (x)=0 from both sides, take the expansion (6) up to the linear order and neglect the corrections of order αa π , we find Numerically, ∆a π ≃ 0.005 [27]. This is the value to be added to the experimental value a exp π extracted neglecting the 1γIR contribution in order to find an estimate of the physical parameter a phys π . Above, δ NLO 1γIR (x) is calculated from δ 1γIR (x, y) stated in Section IV of [9] using prescription (9). Note that the 1γIR contribution was already taken into account in the NA62 analysis [8].
Finally, we can define the ratio R as where we use the NLO radiative corrections in the QED sector to approximate the exact differential width of the Dalitz decay, using at the same time the conservative error estimate. Taking the NA62 prescription f 2 (x) NA62 we finally arrive to

III. CALCULATION AND UNCERTAINTY ESTIMATION
Our aim now is to precisely and reliably (using conservative error estimates) determine R. Based on (12) we obtain During the estimation of the above uncertainty, higher orders in the QED expansion were considered, surpassing in size the uncertainty stemming from the form-factor dependence. In this regard, we can take the absolute value of the dominant correction (R NLO BS,div ; see Table II) as the typical expected maximum value appearing in the NLO order and anticipate that the NNLO correction is suppressed compared to the NLO one in the similar manner as is the NLO suppressed with respect to LO: circa on the level of 3 %. This uncertainty is already conservative: the total NLO correction is on the level of 1 %. Together with the LO value R LO = 11.879(5) × 10 −3 (14) based on the slope a univ π ≡ 0.0355(70) chosen by stretching the uncertainty band over the whole interval of the numerical values suggested by different approaches (see Table II), we finally obtain This is the main result of the presented work. The most significant part of the quoted uncertainty comes from R LO and in particular from a π , since m e , M π and α are known very precisely. Let us point out here that this calculation also includes all the contributions from the decays where additional photon(s) with arbitrarily high (kinematically allowed) energies are present in the final state. Indeed, the bremsstrahlung correction at NLO (calculated a la [9,26]) takes into account an additional photon in the final state and integrates over its energy and angle of emission without any additional cuts. On the other hand, the used approach allows for recalculation  To be suitable for interpolation, higher precision is used. The quoted uncertainties are dominated by the knowledge of the form-factor slope (we assume aπ = 0.0355(70)); the additional 3% uncertainty covering the higher-order corrections is also taken into account. Note the different additional multiplicative factors depending on xcut.
of the prediction for R, if needed, using such a cut in the bremsstrahlung correction δ BS (x, y). In other words, the results stated in this paper are meant to be used for inclusive process and this is how they should be understood. However, similar quantities for exclusive processes can be obtained in a similar way while introducing some specific cut-off. A combined approach was used in the analysis of the recent NA62 measurement [8], when an additional photon was simulated above the cut-off given by the detector sensitivity. Note that we have explicitly stated that the results include an additional photon in the final state, denoting it as (γ). We also take this tacitly into account in the results for R. To conclude, for each experiment the specific approach for including radiative corrections must be used. In experiments, some specific kinematic regions might be considered. The sample values for are listed in Table III. One can obtain the values for any  intermediate region using the entries of this table. As an example, in the π 0 -rare-decay measurement performed by KTeV [3] the region x > x cut = 0.232 was used for the Dalitz decay, which served as the normalization channel in this search. Direct calculation based on this work leads to R(0.232) = 38.0(2) × 10 −5 and the interpolation based on Table III gives R(0.232)| intpol. = 37.9(2)×10 −5 , which is compatible within uncertainties. To compare with KTeV, in [29] the value [R(0.2319)/R]| KTeV = 0.0319 was used, which is compatible with our calculation: R(0.232)/R = 0.0317 (2). Based on R, we can also predict B(π 0 → γγ) and B(π 0 → e + e − γ) very precisely. Taking into account the uncertainty of R, we can write 1 − B(π 0 → e + e − e + e − ) ≃ B(π 0 → γγ) + B(π 0 → e + e − γ(γ)) , since the branching ratios of other decay modes are smaller than 10 −6 . Using B(π 0 → e + e − e + e − ) = 3.3(2) × 10 −5 [1,30], we find B(π 0 → γγ) ≃ 1 − B(π 0 → e + e − e + e − ) 1 + R = 98.8131(6) % .
(18) Let us note that taking tacitly into account inclusive Dalitz decays in (17) is justified and contributes to the relevant decay modes. Finally the branching ratio for the Dalitz decay reads The above results are compatible with the PDG averages, exhibiting much higher precision. Let us see, as a simple example, how the new result on the Dalitz decay branching ratio (19) influences the measurements of the K + → π + e + e − decay. The lowenergy parameters a + and b + have been measured by the NA48/2 collaboration to be a + = −0.578 (16) and b + = −0.779(66), leading to the model-dependent branching ratio B(K + → π + e + e − ) = 3.11(12) × 10 −7 , using the 2008 PDG average B(π 0 → e + e − γ) = 1.198(32) % [31] for normalization [4]. The central value of our result (19) is 1.2 % lower than the quoted PDG average and has a negligible error. The remaining external uncertainty on the measurement [4] related to the normalization comes from B(K + → π + π 0 ) known to 0.4% precision. The corrected values are a + = −0.575 (14), b + = −0.771(64) and B(K + → π + e + e − ) = 3.07(10) × 10 −7 .

IV. COMPARISON WITH OTHER WORKS
Radiative corrections to the integral decay width were first addressed by Joseph [32], who numerically arrived to the result ∆R = 1.05 × 10 −4 neglecting, among others, the pion form-factor slope. A simple analytical prescription in the limit of vanishing electron mass was later found by Lautrup and Smith [33]: Numerically, ∆R aπ =0 L&S = 1.0378 × 10 −4 and ∆R a univ π L&S = 1.0414(7) × 10 −4 . The two approaches are compatible and should be compared with our result (13). However, the symmetrization with respect to the two photons in the bremsstrahlung contribution [32] or the interference terms therein were neglected [33], respectively. The interference is indeed negligible and corresponds (for a π = 0) to ∆R BS interf = 3.60 × 10 −7 . Similarly, the 1γIR contribution was left out from these calculations. However, even though the 1γIR contribution was considered to be negligible, it is significant compared to ∆R BS interf (see Table II) and embodies the main source of the difference between our result and the previous works [32,33]. Let us stress again that the prediction (15) is based on the complete calculation which includes the entire bremsstrahlung and the 1γIR contributions. Here, form-factor effects were taken into account also in the bremsstrahlung correction [28] and the mass of the final-state leptons was not neglected.
To mention other results on R in literature, let us conclude with the value published in [19]: R = 11.8754(5) × 10 −3 . It is based on the form factor calculated using dispersive formalism, while no radiative corrections were considered; thus it needs to be compatible with (14).

V. DISCUSSION
As the neutral pion represents the lightest hadron, its properties are important for the whole hadronic physics. The main subject of the presented work is connected to the branching ratios of π 0 decay modes. These can further serve to translate the lifetime into decay widths and vice versa. Together with the pion lifetime, the pion decay constant is another important parameter. Although both values can be found in PDG with relatively small errors, a certain caution and prudence are of major importance in this respect.
The pion decay constant F π represents a fundamental order parameter of the spontaneous chiral-symmetry breaking in QCD. Its value can be obtained from the π l2 decay via the combination of F π |V ud | [34]. To obtain the latest prediction one should include both the decay width Γ(π − → µ −ν (γ)) and the V ud parameter as experimental inputs. Note that in pure QCD the difference between F π 0 and F π − is quadratic in m d − m u and can be neglected [35]. Any deviation from the standard V − A coupling of quarks to W should be visible in tension between F π determined from the π ± weak decays and the π 0 lifetime. In other words, the value of F π reflects also the weak sector of SM [36].
The π 0 lifetime is important for the normalization of all π 0 decay modes and consequently other hadronic processes. The first method how to obtain it is connected with the direct measurement of the average distance of the highly-relativistic pion. The second method uses the conserved-vector-current hypothesis which connects the vector form factor (i.e. charged pions) to the π 0 lifetime [37]. Finally, the Primakoff effect can be used. Historically, the π 0 lifetime value has been very stable. However, this is not due to the compatibility of the measurements but rather the lack of the new ones. For 20 years until 2008, the PDG average read the same: 8.4(6) × 10 −17 s [31]. Since 2012 its value settled to 8.52(18) × 10 −17 s [38]. This number is mainly given by the measurement by PrimEx experiment at JLab [39]. The reason for error reduction is the exclusion of two Primakoff-effect-based experiments from 1970. As summarized by PDG, those experiments suffered modelrelated systematics unknown at the time (cf. [40]). To conclude, let us summarize the values of the most precise π 0 -lifetime measurements (×10 −17 s) given by two different methods: 8.32(23) (PrimEx) and 8.97(28) (CERN direct [41]). It is clear that the situation is not satisfactory and a new independent measurement is desirable. Our main result (15) together with the value (18) should decrease considerably the uncertainties in future studies.