text:
# Radiative Corrections to Double Dalitz Decays: Effects on Invariant Mass Distributions and Angular Correlations

## I Background

Meson decays to two photons should exhibit interesting correlations between the photon polarizationsYang (1950). Although existing particle detectors cannot measure photon polarizations directly, it has long been known that the polarization correlations can be measured indirectly by studying angular correlations in the related double Dalitz decaysKroll and Wada (1955) in which both photons undergo internal conversion to a lepton pair. More recently, it has been pointed outUy (1991) that a detailed study of these correlations can be used to determine the relative amount of two possible meson–$`\gamma \gamma `$ couplings (one CP–conserving and one CP–violating for mesons that are CP eigenstates) that can contribute to this process.

Dalitz and double Dalitz decays are also of interest because they can be exploited to perform a measurement of the electromagnetic form factor of the decaying meson — that is, how the meson couples to one real and one virtual photon (Dalitz decay) or two virtual photons (double Dalitz) depends on the $`q^2`$ values of the photon(s). An accurate knowledge of this form factor is essential, for example, to calculate the so–called long–distance contribution to the rare decay $`K_L\mu ^+\mu ^{}`$. The short–distance contribution to this process, mediated by loops involving heavy quarks and massive vector bosons, is sensitive to the CKM matrix element $`V_{td}`$, but this contribution cannot be extracted from the accurate experimental measurement of the partial width unless the long–distance amplitude is precisely known.

The tree–level rates for several double Dalitz processes have been published in various formsKroll and Wada (1955); Hunt (1970); Miyazaki and Takasugi (1973); Uy (1991). The first experimental observation of a double Dalitz decay was published in 1962Plano et al. (1959); Samios et al. (1962). A total of 206 examples of the decay $`\pi ^0e^+e^{}e^+e^{}`$ were observed by Samios in a sample of some 800,000 bubble chamber photographs. Based on the observed angular correlations, Samios was able to exclude the possibility that the $`\pi ^0`$ was a scalar particle (with a CP–conserving decay) at the $`3.3\sigma `$ confidence level. His measurement of the branching ratio for that process remains the only one published to date.

Experimental observations of the much rarer kaon double Dalitz decays began to appear in the 1990’s. Several measurements have been made of the decay $`K_Le^+e^{}e^+e^{}`$Barr et al. (1991); Vagins et al. (1993); Akagi et al. (1993); Gu et al. (1994); Barr et al. (1995); Akagi et al. (1995); Alavi-Harati et al. (2001a), the most recent of which are based on several hundred observed events. The still rarer double Dalitz mode $`K_Le^+e^{}\mu ^+\mu ^{}`$ is particularly interesting because it is free of complications arising when there are two identical lepton pairs in the final state, and because it probes only the kinematic region in which one of the virtual photons has $`q^2>4m_\mu ^2`$. The first example of this decay was reported in 1996Gu et al. (1996). In 2001 the KTeV experiment published a branching ratio based on a sample of 43 eventsAlavi-Harati et al. (2001b); most recently, KTeV has reported results from a combined sample of 132 events, including the earlier 43Hamm (2002).

Further experimental results on both $`\pi ^0e^+e^{}e^+e^{}`$ and the two kaon decays are expected in the near future from the NA48 and KTeV experiments. As the statistics available to the experimenters increase, it will be necessary to have a more accurate theoretical description of these decays, incorporating the significant effects of QED radiative corrections. These corrections, discussed in this paper, have a significant impact on the extraction of both form factors and angular correlations from high–statistics data.

## II Tree–Level Amplitudes

The most general CPT invariant interaction governing the transition of a spin–zero meson to two photons is

$$=\frac{i}{4M}\left[_Pϵ_{\mu \nu \rho \sigma }+_S\left(g_{\mu \rho }g_{\nu \sigma }g_{\mu \sigma }g_{\nu \rho }\right)\right]F^{\mu \nu }F^{\rho \sigma }\varphi ,$$
(1)
where $`_P`$ and $`_S`$ are dimensionless form factors for a pseudoscalar and scalar coupling which may be momentum dependent, $`F^{\mu \nu }`$ is the electromagnetic field tensor, and $`\varphi `$ is the field of a meson of mass $`M`$. The factor of $`i`$ defines a phase convention in which the form factors will be real if CPT invariance holds and there are no absorptive decay amplitudes. It will be convenient to decompose the couplings into real and imaginary parts, and define them (up to an arbitrary overall phase $`\mathrm{\Psi }`$) in terms of a mixing angle $`\zeta `$, and a phase difference $`\delta `$,

$`_P`$ $`=g_Pf_P(k_1^2,k_2^2)e^{i\mathrm{\Psi }}=\stackrel{~}{g}f_P(k_1^2,k_2^2)\mathrm{cos}\zeta e^{i\mathrm{\Psi }},`$ (2a)
$`_S`$ $`=g_Sf_S(k_1^2,k_2^2)e^{i\mathrm{\Phi }}=\stackrel{~}{g}f_S(k_1^2,k_2^2)\mathrm{sin}\zeta e^{i\delta }e^{i\mathrm{\Psi }},`$ (2b)
where $`\stackrel{~}{g}^2=g_P^2+g_S^2`$ and $`\delta =\mathrm{\Phi }\mathrm{\Psi }`$. $`f_P`$ and $`f_S`$ are functions of the $`k^2`$ values of the two virtual photons and are normalized such that $`f(0,0)=1`$ for both couplings. In principle it is also possible for $`\delta `$ to depend on $`k_1^2`$ and $`k_2^2`$.

### II.1 Two–Photon Decay

The two–photon partial width provides information about the coupling constant $`\stackrel{~}{g}`$ along with some constraints on the mixing angle $`\zeta `$ and phase difference $`\delta `$. The matrix element for the decay to two real photons with helicities $`\lambda _1`$ and $`\lambda _2`$ is

$$\begin{array}{cc}\hfill _{\lambda _1\lambda _2}=& \frac{2}{M}[_Pϵ_{\mu \nu \rho \sigma }\hfill \\ & +_S(g_{\mu \rho }g_{\nu \sigma }g_{\mu \sigma }g_{\nu \rho })]k_1^\mu ϵ_{\lambda _1}^\nu k_2^\rho ϵ_{\lambda _2}^\sigma ,\hfill \end{array}$$
(3)
where $`k^\mu `$ and $`ϵ^\mu `$ are the momentum and polarization of a photon. The calculation of the two–photon couplings for photons of arbitrary mass is carried out in Appendix B. For real photons, the kinematic factors in the couplings reduce to $`\lambda =z=1`$ and $`w=0`$, and the momentum dependent functions $`f_P`$ and $`f_S`$ reduce to unity. Therefore, one has the following two contributions

$`_{++}`$ $`=M\stackrel{~}{g}\left(\mathrm{sin}\zeta e^{i\delta }i\mathrm{cos}\zeta \right),`$ (4a)
$`_{}`$ $`=M\stackrel{~}{g}\left(\mathrm{sin}\zeta e^{i\delta }+i\mathrm{cos}\zeta \right).`$ (4b)
The partial width is obtained by integrating $`1/(2M)`$ times the squared matrix element over the available phase space. The matrix element is a constant, so the integration results in a factor of $`1/(16\pi )`$. The partial widths for the two allowed helicity states are then

$`\mathrm{\Gamma }_{++}`$ $`={\displaystyle \frac{M\stackrel{~}{g}^2}{32\pi }}\left(12\mathrm{sin}\zeta \mathrm{cos}\zeta \mathrm{sin}\delta \right),`$ (5a)
$`\mathrm{\Gamma }_{}`$ $`={\displaystyle \frac{M\stackrel{~}{g}^2}{32\pi }}\left(1+2\mathrm{sin}\zeta \mathrm{cos}\zeta \mathrm{sin}\delta \right).`$ (5b)
The decay rates to the two final states will be identical if either the phase difference between the two couplings is zero or one of the form factors is zero. In any case, the total rate is equal to $`\mathrm{\Gamma }_{\gamma \gamma }=M\stackrel{~}{g}^2/(16\pi )`$. An experimental measure of the two–photon width then gives a value of $`\stackrel{~}{g}`$,

$$\stackrel{~}{g}=\sqrt{16\pi \mathrm{\Gamma }_{\gamma \gamma }/M}.$$
(6)
### II.2 Application to Neutral Pion and Kaon Decays

One can compute $`\stackrel{~}{g}`$ for the $`\pi ^0\gamma \gamma `$, $`K_L\gamma \gamma `$, and $`K_S\gamma \gamma `$ decays using the current experimental valuesHagiwara et al. (2002) of the two–photon branching ratios along with the meson lifetimes and masses, with the results

$`\stackrel{~}{g}_0`$ $`=(1.70\pm 0.06)10^3,`$ (7a)
$`\stackrel{~}{g}_L`$ $`=(8.75\pm 0.12)10^{10},`$ (7b)
$`\stackrel{~}{g}_S`$ $`=(1.36\pm 0.11)10^9.`$ (7c)
It will be useful to define CP–even and CP–odd $`\gamma \gamma `$ states:

$$|(\gamma \gamma )_\pm =(|++\pm |)/\sqrt{2}.$$
(8)
The matrix elements for $`\pi ^0\gamma \gamma `$ are then given by

$`(\gamma \gamma )_+|T|\pi ^0`$ $`=i\sqrt{2}M_\pi \stackrel{~}{g}_0\mathrm{sin}\zeta _0e^{i\mathrm{\Phi }_0},`$ (9a)
$`(\gamma \gamma )_{}|T|\pi ^0`$ $`=\sqrt{2}M_\pi \stackrel{~}{g}_0\mathrm{cos}\zeta _0e^{i\mathrm{\Psi }_0}.`$ (9b)
Similar expressions can be given for the $`K_1`$ and $`K_2`$ decay matrix elements with the 0 subscripts changed to 1 or 2 respectively and $`M_\pi `$ replaced by $`M_K`$. These expressions may be compared to those in Ref. Sehgal (1971) in which CP–violation observables for the $`K\gamma \gamma `$ decays were calculated. In that article Sehgal defines

$`(\gamma \gamma )_+|T|K_1`$ $`=c_ee^{i\rho _e},`$ (10a)
$`(\gamma \gamma )_{}|T|K_1`$ $`=ic_oe^{i\rho _o},`$ (10b)
$`(\gamma \gamma )_+|T|K_2`$ $`=id_ee^{i\mu _e},`$ (10c)
$`(\gamma \gamma )_{}|T|K_2`$ $`=d_oe^{i\mu _o}.`$ (10d)
In our notation we have

$`c_e`$ $`=\sqrt{2}M_K\stackrel{~}{g}_1\mathrm{sin}\zeta _1,`$ (11a)
$`c_o`$ $`=\sqrt{2}M_K\stackrel{~}{g}_1\mathrm{cos}\zeta _1,`$ (11b)
$`d_e`$ $`=\sqrt{2}M_K\stackrel{~}{g}_2\mathrm{sin}\zeta _2,`$ (11c)
$`d_o`$ $`=\sqrt{2}M_K\stackrel{~}{g}_2\mathrm{cos}\zeta _2.`$ (11d)
The observable phase differences are given by $`\delta _1=\rho _e\rho _o`$ and $`\delta _2=\mu _e\mu _o`$.

The phases $`\delta _i`$ are observable only if CP is violated in $`\gamma \gamma `$ decays. As Sehgal noted in Ref. Sehgal (1971), the phase $`\delta _0`$ in $`\pi ^0`$ decays should be very close to zero if CPT is conserved because the relevant absorptive amplitudes are of order $`\alpha ^2`$. In the case of kaon decays, on–shell intermediate states such as $`\pi \pi `$ couple strongly so that the phases $`\delta _1`$ and $`\delta _2`$ may be large even if CPT is conserved.

### II.3 Four–Lepton Decay

While the total two–photon decay rate provides information about the constant part of the form factors, the four–lepton rate can be utilized to probe both the momentum dependence of the form factors and the mixing of the couplings. The decay to four leptons may contain either two pairs of identical particles (e.g., $`K_Le^+e^{}e^+e^{}`$) or pairs of non–identical particles (e.g., $`K_Le^+e^{}\mu ^+\mu ^{}`$). The tree–level Feynman diagram for the decay is depicted in Fig. 1.

If the final state does contain identical particles, then there is also an exchange diagram, as shown in Fig. 2.

The matrix element is then the sum of the two diagrams $`=_1+_2`$ and its square is $`||^2=|_1|^2+|_2|^2+2\mathrm{Re}(_1^{}_2)`$. The analysis presented in the rest of this section applies only to the direct contribution to the double Dalitz process. Appendix C gives an explicit expression for the interference term.

The matrix element for the direct contribution to the double Dalitz process can be written as

$$_1=H_{\mu \nu \rho \sigma }k_1^\mu \mathrm{\Pi }^{\nu \alpha }\mathrm{\Gamma }_\alpha k_2^\rho \mathrm{\Pi }^{\sigma \beta }\mathrm{\Gamma }_\beta ,$$
(12)
where $`H`$ is a two–photon coupling given by Eq. (94), $`k`$ is a photon momentum, $`\mathrm{\Pi }`$ is the propagator for a photon of momentum $`k`$

$$\mathrm{\Pi }^{\mu \nu }(k)=\frac{i}{k^2}\left(\underset{\lambda }{}ϵ_\lambda ^\mu ϵ_\lambda ^\nu \widehat{k}^\mu \widehat{k}^\nu \right),$$
(13)
and $`\mathrm{\Gamma }`$ is the fermion current for an electron of momentum $`q`$ and spin $`r`$ and positron of momentum $`p`$ and spin $`s`$

$$\mathrm{\Gamma }^\mu (q,p;r,s)=ie\overline{u}_r(q)\gamma ^\mu v_s(p).$$
(14)
The sum in the propagator extends over the three helicity states and the term proportional to $`k^\mu k^\nu `$ vanishes when contracted with the current. The matrix element can then be cast as

$$_1=\underset{\lambda _1\lambda _2}{}_{\lambda _1\lambda _2}L_{\lambda _1}L_{\lambda _2},$$
(15)
where $``$ is one of the two–photon couplings given in Eqs. (96,97). The quantity $`L`$ contains all the lepton information and is equal to $`L_\lambda =iϵ_\lambda ^\mu \mathrm{\Gamma }_\mu /k^2`$.

The squared matrix element, summed over final state helicities can be expressed in terms of the five phase space variables $`x_{12}`$, $`x_{34}`$, $`y_{12}`$, $`y_{34}`$, and $`\varphi `$ (defined in Appendix A). Focusing on the $`\varphi `$ dependence, the squared amplitude is

$$\begin{array}{cc}\hfill |_1|^2=\frac{2^8\pi ^2\alpha ^2\stackrel{~}{g}^2}{M^2w^4}& (A\mathrm{sin}^2\varphi +B\mathrm{cos}^2\varphi \hfill \\ & +C\mathrm{sin}\varphi \mathrm{cos}\varphi +D\mathrm{sin}\varphi \hfill \\ & +E\mathrm{cos}\varphi +F),\hfill \end{array}$$
(16)
where

$`A=`$ $`w^2\{f_P^2\mathrm{cos}^2\zeta \lambda ^2[1+(1\lambda _{12}^2+y_{12}^2)(1\lambda _{34}^2+y_{34}^2)]`$
$`+f_S^2\mathrm{sin}^2\zeta z^2[2\lambda _{12}^2\lambda _{34}^2+y_{12}^2+y_{34}^2]\},`$ (17a)
$`B=`$ $`w^2\{f_S^2\mathrm{sin}^2\zeta z^2[1+(1\lambda _{12}^2+y_{12}^2)(1\lambda _{34}^2+y_{34}^2)]`$
$`+f_P^2\mathrm{cos}^2\zeta \lambda ^2[2\lambda _{12}^2\lambda _{34}^2+y_{12}^2+y_{34}^2]\},`$ (17b)
$`C=`$ $`2f_Pf_S\mathrm{sin}\zeta \mathrm{cos}\zeta \mathrm{cos}\delta \lambda zw^2`$
$`\times (\lambda _{12}^2y_{12}^2)(\lambda _{34}^2y_{34}^2),`$ (17c)
$`D=`$ $`2f_Pf_S\mathrm{sin}\zeta \mathrm{cos}\zeta \mathrm{cos}\delta \lambda w^3y_1y_2`$
$`\times \sqrt{(\lambda _{12}^2y_{12}^2)(\lambda _{34}^2y_{34}^2)},`$ (17d)
$`E=`$ $`2f_S^2\mathrm{sin}^2\zeta zw^3y_{12}y_{34}\sqrt{(\lambda _{12}^2y_{12}^2)(\lambda _{34}^2y_{34}^2)},`$ (17e)
$`F=`$ $`f_S^2\mathrm{sin}^2\zeta w^4(1y_{12}^2)(1y_{34}^2),`$ (17f)
where the kinematic variables $`w`$, $`z`$, and $`\lambda `$ are defined in Appendix A. Terms $`A`$ and $`B`$ arise from the diagonal pieces of the transverse part of the pseudoscalar and scalar couplings while term $`F`$ is due to the diagonal piece of the longitudinal part of the scalar coupling. Term $`C`$ is the interference between the transverse parts of the pseudoscalar and scalar couplings while terms $`D`$ and $`E`$ are due to interference between the longitudinal part of the scalar and the transverse parts of the pseudoscalar and scalar couplings respectively.

The partial width is obtained by integrating $`1/(2M)`$ times the square of the matrix element over the eight dimensional phase space (given in Appendix A). The integrals over $`y_{ij}`$ extend from $`\lambda _{ij}`$ to $`+\lambda _{ij}`$ and can be done analytically. The differential partial width, normalized to the two–photon width and integrated over $`y_{12}`$ and $`y_{34}`$, is

$$\begin{array}{cc}\hfill \frac{1}{\mathrm{\Gamma }_{\gamma \gamma }}\frac{d\mathrm{\Gamma }_1}{d\varphi }=\frac{𝒮\alpha ^2}{3\pi ^3}& [I_1\mathrm{cos}^2\zeta \mathrm{sin}^2\varphi +I_2\mathrm{sin}^2\zeta \mathrm{cos}^2\varphi \hfill \\ & +I_3\mathrm{sin}\zeta \mathrm{cos}\zeta \mathrm{cos}\delta \mathrm{sin}\varphi \mathrm{cos}\varphi \hfill \\ & +I_4\mathrm{cos}^2\zeta +(I_5+I_6)\mathrm{sin}^2\zeta ],\hfill \end{array}$$
(18)
where $`𝒮`$ is a symmetry factor which is $`1/4`$ for modes with identical particles and 1 otherwise.

The $`I_i`$ factors represent the integrals over $`x_{12}`$ and $`x_{34}`$ given below. The factors $`I_1`$ and $`I_4`$ correspond to the pseudoscalar coupling, $`I_2`$ and $`I_5`$ are the analogous terms for the scalar coupling, $`I_6`$ is the additional longitudinal term in the scalar coupling, and $`I_3`$ is the interference term.

$`I_1=`$ $`{\displaystyle \frac{2}{3}}{\displaystyle 𝑑x_{12}𝑑x_{34}f_P^2\frac{\lambda _{12}^3\lambda _{34}^3\lambda ^3}{w^2}},`$ (19a)
$`I_2=`$ $`{\displaystyle \frac{2}{3}}{\displaystyle 𝑑x_{12}𝑑x_{34}f_S^2\frac{\lambda _{12}^3\lambda _{34}^3\lambda z^2}{w^2}},`$ (19b)
$`I_3=`$ $`{\displaystyle \frac{4}{3}}{\displaystyle 𝑑x_{12}𝑑x_{34}f_Pf_S\frac{\lambda _{12}^3\lambda _{34}^3\lambda ^2z}{w^2}},`$ (19c)
$`I_4=`$ $`{\displaystyle 𝑑x_{12}𝑑x_{34}f_P^2\frac{\lambda _{12}\lambda _{34}\lambda ^3}{w^2}(3\lambda _{12}^2\lambda _{34}^2)},`$ (19d)
$`I_5=`$ $`{\displaystyle 𝑑x_{12}𝑑x_{34}f_S^2\frac{\lambda _{12}\lambda _{34}\lambda z^2}{w^2}(3\lambda _{12}^2\lambda _{34}^2)},`$ (19e)
$`I_6=`$ $`{\displaystyle \frac{1}{6}}{\displaystyle 𝑑x_{12}𝑑x_{34}f_S^2\lambda _{12}\lambda _{34}\lambda (3\lambda _{12}^2)(3\lambda _{34}^2)}.`$ (19f)
The double integral is performed by first integrating over $`x_{34}`$ from $`x_{34}^0`$ to $`(1\sqrt{x_{12}})^2`$ and then over $`x_{12}`$ from $`x_{12}^0`$ to $`(1\sqrt{x_{34}^0})^2`$, where $`x_{ij}^0=(m_i+m_j)^2/M^2`$. In order to obtain values for these integrals, $`f_P`$ and $`f_S`$ must first be specified and then the integrals can be done numerically. Table 1 summarizes the values for the different double Dalitz modes assuming that $`f_P=f_S=1`$.

The numerical value of $`I_5+I_6`$ was found to be several orders of magnitude larger in Ref. Uy (1991). Extracting $`I_6`$ from that result yields a value of 3578.0, compared with our value of 0.01503. This discrepancy has been traced to the use in Ref. Uy (1991) of a scalar coupling $`g_{\mu \rho }g_{\nu \sigma }`$, rather than the correctly antisymmetrized tensor $`g_{\mu \rho }g_{\nu \sigma }g_{\mu \sigma }g_{\nu \rho }`$. One consequence of the much smaller value we find for $`I_6`$ is that the total width for the double Dalitz decay is almost completely insensitive to the mix of scalar and pseudoscalar couplings (except for the $`K_L\mu ^+\mu ^{}\mu ^+\mu ^{}`$ decay), in contradistinction to the conclusion of Ref. Uy (1991) but in agreement with the comments near the end of Ref. Kroll and Wada (1955).

The differential rate can also be expressed in a compact form, suitable for experimental fits to the $`\varphi `$ distribution, involving a constant term, a CP conserving $`\mathrm{cos}2\varphi `$ term, and a CP violating $`\mathrm{sin}2\varphi `$ term

$$\frac{1}{\mathrm{\Gamma }_{\gamma \gamma }}\frac{d\mathrm{\Gamma }_1}{d\varphi }=\frac{\alpha ^2}{3\pi ^3}R(1+\kappa _1\mathrm{cos}2\varphi +\kappa _2\mathrm{sin}2\varphi ),$$
(20)
where

$`R=`$ $`𝒮[(I_1/2+I_4)\mathrm{cos}^2\zeta `$
$`+(I_2/2+I_5+I_6)\mathrm{sin}^2\zeta ],`$ (21a)
$`\kappa _1=`$ $`𝒮(I_2\mathrm{sin}^2\zeta I_1\mathrm{cos}^2\zeta )/(2R),`$ (21b)
$`\kappa _2=`$ $`𝒮I_3\mathrm{sin}\zeta \mathrm{cos}\zeta \mathrm{cos}\delta /(2R).`$ (21c)
The values of $`R`$ and $`\kappa _1`$ at $`\zeta =0`$ and $`\zeta =\pi /2`$, along with the maximum value of $`\kappa _2`$ and the angle $`\zeta _0`$ at which it takes on that value, are listed in Table 2.

As expected, for a pure pseudoscalar decay, the amplitude of the $`\mathrm{cos}2\varphi `$ term will be negative while the amplitude of the $`\mathrm{sin}2\varphi `$ vanishes. For a pure scalar decay, the amplitude of the $`\mathrm{cos}2\varphi `$ term is nearly the same magnitude as in the pseudoscalar decay but positive, and the amplitude of $`\mathrm{sin}2\varphi `$ again vanishes. Depending on the mode, the amplitude for the CP–violating $`\mathrm{sin}2\varphi `$ term is maximal for values of the mixing angle between $`\pi /8`$ and $`\pi /4`$.

Alternatively, we could have integrated over $`\varphi `$ before $`x_{12}`$ and $`x_{34}`$, in which case we would have

$$\begin{array}{c}\frac{1}{\mathrm{\Gamma }_{\gamma \gamma }}\frac{d^2\mathrm{\Gamma }_1}{dx_{12}dx_{34}}=\frac{2𝒮\alpha ^2}{9\pi ^2}\frac{\lambda _{12}\lambda _{34}\lambda }{w^2}(3\lambda _{12}^2)(3\lambda _{34}^2)\hfill \\ \hfill \times [f_P^2\mathrm{cos}^2\zeta \lambda ^2+f_S^2\mathrm{sin}^2\zeta (\lambda ^2+3w^2/2)],\end{array}$$
(22)
where we have used $`z^2=\lambda ^2+w^2`$. The interference term integrates to zero and what remains clearly shows the kinematic differences between the contributions of the two couplings.

Going back to Eq. (20), the final integral over $`\varphi `$ from $`0`$ to $`2\pi `$ can be performed to get the direct contribution to the four–lepton decay rate relative to the two–photon rate

$$\frac{\mathrm{\Gamma }_1}{\mathrm{\Gamma }_{\gamma \gamma }}=\frac{2\alpha ^2}{3\pi ^2}R.$$
(23)
The total tree–level rate can now be computed for arbitrary form factors for modes without identical particles in the final state. For the other modes, there is the interference between the direct and exchange graphs that must be included. The decay rate has the form

$$\mathrm{\Gamma }=\mathrm{\Gamma }_1+\mathrm{\Gamma }_2+\mathrm{\Gamma }_{12},$$
(24)
where, for modes without identical particles $`\mathrm{\Gamma }_2=\mathrm{\Gamma }_{12}=0`$, and for modes with identical particles $`\mathrm{\Gamma }_2=\mathrm{\Gamma }_1`$. The expression given in Appendix C for the interference term could in principle be integrated numerically. We choose instead to use a Monte Carlo (MC) simulation to integrate the rate and make histograms of the relevant phase space variables. The decay rates for the various modes, broken into diagonal and interference terms, are listed in Table 3.

Ref. Miyazaki and Takasugi (1973) included a similar table of values, some of which are in disagreement with our results. The most significant discrepancy involves the size of the interference term for the decays $`K_Le^+e^{}e^+e^{}`$ and $`\pi ^0e^+e^{}e^+e^{}`$. We find that the interference in $`K_Le^+e^{}e^+e^{}`$ is roughly 9 times smaller than Ref. Miyazaki and Takasugi (1973) reports, and that the interference in $`\pi ^0e^+e^{}e^+e^{}`$ is about 4 times smaller. We also differ in the total rate for $`K_Le^+e^{}\mu ^+\mu ^{}`$, but the factor of 2 difference is likely due to a typographical error in the previous publication.

The assumption that the form factor is flat contradicts current experimental findings. The two models that have been used to parameterize the kaon form factor are the BMS model Bergström et al. (1990) and the DIP model D’Ambrosio et al. (1998). The BMS model was originally proposed to describe the coupling in the single Dalitz decay $`K_Le^+e^{}\gamma `$ and was therefore written as a function of one $`k^2`$ only

$`f(x,0)`$ $`={\displaystyle \frac{1}{1r_\rho x}}+{\displaystyle \frac{C\alpha _K^{}}{1r_K^{}x}}[{\displaystyle \frac{4}{3}}{\displaystyle \frac{1}{1r_\rho x}}`$
$`{\displaystyle \frac{1}{9}}({\displaystyle \frac{1}{1r_\omega x}}+{\displaystyle \frac{2}{1r_\varphi x}})].`$ (25)
The quantities $`r_i=M^2/M_i^2`$ for $`M_i`$ equal to the $`\rho `$, $`K^{}`$, $`\omega `$, or $`\varphi `$ meson masses. To apply this model to the double Dalitz decay, it is assumed that the coupling factors so that $`f(x_1,x_2)=f(x_1,0)f(x_2,0)`$. In this paper, we will use a simplified form of the DIP form factor, which involves only the $`\rho `$ meson and two parameters

$`f(x_1,x_2)`$ $`=1+\alpha _{DIP}\left({\displaystyle \frac{x_1}{x_1M_\rho ^2/M^2}}+{\displaystyle \frac{x_2}{x_2M_\rho ^2/M^2}}\right)`$
$`+\beta _{DIP}{\displaystyle \frac{x_1x_2}{(x_1M_\rho ^2/M^2)(x_2M_\rho ^2/M^2)}}.`$ (26)
As will be seen in Appendix D, the BMS model can be expressed as a generalized DIP model involving the $`\rho `$, $`\omega `$, and $`\varphi `$ vector mesons.

Experimentally, the form factor has traditionally been linearized in the case of the pion with just a slope parameter measured, while for the kaon, the BMS model has been used and values of $`\alpha _K^{}`$ quoted. The conversion to the DIP parameters is easily done, using the world average Hagiwara et al. (2002) for the kaon we will use $`\alpha _{DIP}=1.5`$ and for the pion, $`\alpha _{DIP}=1.0`$. There is as yet no experimental sensitivity to $`\beta _{DIP}`$ and so we will use $`\beta _{DIP}=0`$. The effect of using the DIP model with these values of $`\alpha `$ is that the $`\pi ^0e^+e^{}e^+e^{}`$ rate increases by less than $`0.4\%`$, the $`K_Le^+e^{}e^+e^{}`$ rate increases by $`6.5\%`$, the $`K_Le^+e^{}\mu ^+\mu ^{}`$ rate increases by $`56\%`$, and the $`K_L\mu ^+\mu ^{}\mu ^+\mu ^{}`$ rate increases by $`68\%`$. It is clear that the assumption of a flat form factor is completely invalid for modes containing muons.

## III Higher Order Processes

The tree–level double Dalitz process is $`𝒪(\alpha ^2)`$ since it contains two electromagnetic vertices. Higher order contributions to the double Dalitz rate contain one or more internal loops. There are three types of graphs that contribute at $`𝒪(e^4)`$: the vacuum polarization, the vertex correction, and the 5–point diagram. A representative diagram from each of these processes is displayed in Figs. 6,7, and 8, respectively. There are two graphs for both the vacuum polarization and the vertex correction, one for each pair, plus four graphs for the 5–point function. If there are identical particles in the final state, there are exchange diagrams and the number of graphs doubles. The interference between the tree–level diagram and the one–loop diagrams is $`𝒪(\alpha ^3)`$ and therefore contributes to the first order radiative correction to the double Dalitz rate.

Both the vertex correction and the 5–point graph contain IR divergences, that is, divergences in the limit that the exchanged photon energy goes to zero. In order to handle this singular behavior, one must also consider the radiative double Dalitz decay $`Xl_1^+l_1^{}l_2^+l_2^{}\gamma `$, in which one of the leptons internally radiates a photon. There are two contributions to this process, shown in Figs. 4 and 5. The radiative process diverges in the opposite manner from the one–loop graphs making the combined decay rate finite.

The combined process will be indicated by $`Xl_1^+l_1^{}l_2^+l_2^{}(\gamma )`$, where the radiated photon may or may not be detectable. The distinction between the non–radiative and radiative decays is an experimental issue and ultimately related to the hardware. We will use the energy of the radiated photon in the CM frame to differentiate between events with a hard photon, $`E_\gamma >E_{\text{cut}}`$, and events without, $`E_\gamma <E_{\text{cut}}`$. The cutoff is chosen such that photons with energies below the cutoff can have no significant effect on the 4–lepton acceptance. The rate for radiative events with soft photons will be added to the rate for non–radiative decays. This contribution is also $`𝒪(\alpha ^3)`$ and therefore must be considered along with the one–loop corrections.

The double Dalitz differential rate to second order can therefore be expressed as

$$d^5\mathrm{\Gamma }_{\text{rad}}=d^5\mathrm{\Gamma }_{\text{tree}}(1+\delta _{\text{brem}}+\delta _{\text{virt}}),$$
(27)
where $`\delta _{\text{brem}}`$ is the bremsstrahlung contribution due to radiative decays with photons below the photon energy cutoff and $`\delta _{\text{virt}}`$ is the virtual correction due to the interference between the tree–level and one–loop diagrams. The virtual correction can be further decomposed into the contributions from the three one–loop diagrams

$$\delta _{\text{virt}}=\delta _{\text{vp}}+\delta _{\text{vc}}+\delta _{\text{5p}},$$
(28)
where $`\delta _{\text{vp}}`$ is the correction from the vacuum polarization diagrams, $`\delta _{\text{vc}}`$ is the correction from the vertex correction diagrams, and $`\delta _{\text{5p}}`$ is the correction from the 5–point diagrams.

## IV Radiative Decays

The radiative double Dalitz decay will only be considered at tree–level. It is straightforward but tedious to write down the expression for the rate. The two contributions to the rate are shown in Figs. 4 and 5. For each process, there exists three additional diagrams where the photon is radiated off of the other leptons, plus four exchange graphs if applicable.

Our results for the radiative decay rates use a MC simulation in which we calculate the amplitudes for each helicity state using explicit representations of the spinors and polarization vector. An photon energy cutoff of $`400\text{keV}`$ in the CM frame is used for kaon decays while for pions, a cutoff of $`100\text{keV}`$ is used. It is useful to define the quantity $`x_{4e}=m_{4e}^2/M^2`$, where $`m_{4e}`$ is the reconstructed four–lepton invariant mass, to distinguish between the radiative and non–radiative processes. In terms of this variable, the cutoff for both kaon and pion decays is at $`x_{4e}^{\text{cut}}0.9985`$. Fig. 3 shows the distribution of $`x_{4e}`$ for $`K_Le^+e^{}e^+e^{}`$ and $`K_Le^+e^{}e^+e^{}\gamma `$ events.

The large peak at $`x_{4e}=1`$ is due to non–radiative events. The part of the distribution which falls away from the peak at 1 is due to radiated photons from the process of Fig. 4. The rising part of the distribution near $`x_{4e}=0`$ is due to hard Dalitz photons from the process of Fig. 5. Lost due to bin size is the low energy cutoff at $`E_\gamma =400\text{keV}`$ and the high energy cutoff at $`x_{4e}=16m_e^2/M^2`$.

Table 4 shows the tree–level radiative decay rates for the four modes, with no form factor.

The rates in the first column include photons of all energies, while the rates in the second column include only photons with energies large enough that $`x_{4e}<0.95`$. This value of $`x_{4e}`$ is chosen to closely match the resolution on the four–lepton mass in current experiments.

## V Radiative Corrections

The next four sections will describe the different contributions to the radiative corrections to the double Dalitz differential rate. The first three contributions are relatively straightforward to determine and we will therefore only summarize the relevant formulas. The last contribution, coming from the 5–point diagram, is considerably more difficult to calculate. In particular, numerical instabilities plague the evaluation of the tensor 5–point integrals involving light leptons. The fourth section, along with much of the appendix, will outline our procedure for obtaining this (usually) small but non–negligible contribution. We will present the full 5–point diagram corrections to the differential rate in closed form. In Ref. van Neerven and Vermaseren (1984), van Neerven and Vermaseren reported a numerical integral of the radiative corrections to the related two–photon process $`e^+e^{}e^+e^{}\pi ^0`$ but did not present the corrections to the differential cross section. As we will discuss in Sec. VI, the effects of radiative corrections are much more important compared to form factor effects when considering double Dalitz decays, because the $`q^2`$ values in the accessible phase space are much smaller than those typically probed in two–photon resonance formation with final state lepton tags.

### V.1 Bremsstrahlung Correction

The contribution to the double Dalitz differential decay rate due to the soft bremsstrahlung part of the radiative decay is defined as

$$\delta _{\text{brem}}(x_{12},x_{34},y_{12},y_{34},\varphi )=\frac{d^5\mathrm{\Gamma }_{\text{brem}}/d^5\mathrm{\Phi }}{d^5\mathrm{\Gamma }_{\text{tree}}/d^5\mathrm{\Phi }},$$
(29)
where $`d^5\mathrm{\Gamma }_{\text{brem}}/d^5\mathrm{\Phi }`$ is the differential decay rate for the soft part of the radiative decay integrated over the photon momentum with the constraint $`E_\gamma <E_{\text{cut}}`$. The full differential rate is

$$d^8\mathrm{\Gamma }_{\text{brem}}=\frac{1}{2M}|_{\text{brem}}|^2d^5\mathrm{\Phi }\frac{d^3k}{(2\pi )^32E_k},$$
(30)
where $`d^5\mathrm{\Phi }`$ is the four–body phase space differential and $`_{\text{brem}}`$ is the matrix element for the soft bremsstrahlung contribution. If the photon energy cutoff is taken small enough, the matrix element can be approximated as

$$_{\text{brem}}=e\left(\frac{p_2ϵ}{p_2k}+\frac{p_4ϵ}{p_4k}\frac{p_1ϵ}{p_1k}\frac{p_3ϵ}{p_3k}\right)_{\text{tree}},$$
(31)
where $`ϵ`$ and $`k`$ are the radiated photon’s polarization and momentum 4–vectors, respectively. There is one contribution from each of the radiative diagrams represented by Fig. 4. The other type of radiative process (Fig. 5) does not contribute in this limit.

If the cutoff is small enough, the lepton momenta can be held fixed while the photon momentum is integrated out, with the result

$$\delta _{\text{brem}}=4\pi \alpha _0^{E_{\text{cut}}}\frac{d^3k}{(2\pi )^32E_k},$$
(32)
where

$``$ $`={\displaystyle \underset{ϵ}{}}\left|\left({\displaystyle \frac{p_2^\mu }{p_2k}}+{\displaystyle \frac{p_4^\mu }{p_4k}}{\displaystyle \frac{p_1^\mu }{p_1k}}{\displaystyle \frac{p_3^\mu }{p_3k}}\right)ϵ_\mu \right|^2,`$ (33)
$`\begin{array}{c}\hfill ={\displaystyle \frac{2p_1p_2}{(p_1k)(p_2k)}}+{\displaystyle \frac{2p_3p_4}{(p_3k)(p_4k)}}+{\displaystyle \frac{2p_1p_4}{(p_1k)(p_4k)}}\\ \hfill +{\displaystyle \frac{2p_2p_3}{(p_2k)(p_3k)}}{\displaystyle \frac{2p_2p_4}{(p_2k)(p_4k)}}{\displaystyle \frac{2p_1p_3}{(p_1k)(p_3k)}}\\ \hfill {\displaystyle \frac{p_1^2}{(p_1k)^2}}{\displaystyle \frac{p_2^2}{(p_2k)^2}}{\displaystyle \frac{p_3^2}{(p_3k)^2}}{\displaystyle \frac{p_4^2}{(p_4k)^2}}.\end{array}`$
The correction can be expressed in terms of a sum of ten integrals which can be done in closed form

$$\begin{array}{c}\delta _{\text{brem}}=4\pi \alpha [2I(p_1,p_2)+2I(p_3,p_4)+2I(p_1,p_4)\hfill \\ \hfill +2I(p_2,p_3)2I(p_2,p_4)2I(p_1,p_3)\\ \hfill I(p_1,p_1)I(p_2,p_2)I(p_3,p_3)I(p_4,p_4)],\end{array}$$
(34)
where

$$I(k_i,k_j)=_0^{E_{\text{cut}}}\frac{d^3k}{(2\pi )^32E_k}\frac{k_ik_j}{(k_ik)(k_jk)}.$$
(35)
Each integral yields both a finite part and an IR divergent part which goes as $`\mathrm{ln}(2E_{\text{cut}}/\mathrm{\Lambda })`$ where $`\mathrm{\Lambda }`$ is the photon mass which will be taken to zero after the divergent terms are canceled against each other. The first two divergent terms will be seen to cancel the divergent parts of the vertex correction, the next four cancel divergences in the 5–point functions, and the last four cancel the electron self–energy divergences (which are included in the renormalized vertex function). For $`k_ik_j`$

$`I(k_i,`$ $`k_j)={\displaystyle \frac{z_{ij}}{8\pi ^2\lambda _{ij}}}\{\mathrm{ln}\left({\displaystyle \frac{z_{ij}+\lambda _{ij}}{z_{ij}\lambda _{ij}}}\right)\mathrm{ln}\left({\displaystyle \frac{2E_{\text{cut}}}{\mathrm{\Lambda }}}\right)`$
$`+{\displaystyle \frac{1}{4}}\mathrm{ln}^2\left({\displaystyle \frac{\mathrm{\Omega }_i^{}}{\mathrm{\Omega }_i^+}}\right){\displaystyle \frac{1}{4}}\mathrm{ln}^2\left({\displaystyle \frac{\mathrm{\Omega }_j^{}}{\mathrm{\Omega }_j^+}}\right)`$
$`+\mathrm{Li}_2\left(1{\displaystyle \frac{\mathrm{{\rm Y}}_{ij}\mathrm{\Omega }_i^+}{x_{ij}\lambda _{ij}}}\right)+\mathrm{Li}_2\left(1{\displaystyle \frac{\mathrm{{\rm Y}}_{ij}\mathrm{\Omega }_i^{}}{x_{ij}\lambda _{ij}}}\right)`$
$`\mathrm{Li}_2(1{\displaystyle \frac{\mathrm{{\rm Y}}_{ij}\mathrm{\Omega }_j^+}{x_{ij}\lambda _{ij}}})\mathrm{Li}_2(1{\displaystyle \frac{\mathrm{{\rm Y}}_{ij}\mathrm{\Omega }_j^{}}{x_{ij}\lambda _{ij}}})\},`$ (36)
where

$`\mathrm{\Omega }_i^\pm `$ $`=(1+\delta _{i,jkl}\pm \lambda _{i,jkl})/2,`$ (37a)
$`\mathrm{\Omega }_j^\pm `$ $`=(1+\delta _{j,ikl}\pm \lambda _{j,ikl})/(2\sigma _{ij}),`$ (37b)
$`\mathrm{{\rm Y}}_{ij}`$ $`=\sigma _{ij}(1+\delta _{i,jkl})(1+\delta _{j,ikl}),`$ (37c)
$`\sigma _{ij}`$ $`=(z_{ij}+\lambda _{ij})/(1z_{ij}+\delta _{ij}),`$ (37d)
and the various $`\delta `$, $`z`$, and $`\lambda `$ symbols are defined in Appendix A. For $`k_i=k_j`$

$$I(k_i,k_i)=\frac{1}{4\pi ^2}\left[\mathrm{ln}\left(\frac{2E_{\text{cut}}}{\mathrm{\Lambda }}\right)\frac{1}{2\lambda _{ii}}\mathrm{ln}\left(\frac{1+\lambda _{ii}}{1\lambda _{ii}}\right)\right],$$
(38)
where $`\lambda _{ii}`$ is again defined in Appendix A.

It will be enlightening to extract the IR divergent part of the soft brem contribution and express it in a way that will make the cancellation obvious. Collecting terms, one has

$$\begin{array}{c}\delta _{\text{brem}}^{IR}=\mathrm{ln}\mathrm{\Lambda }\{\frac{2\alpha }{\pi }[1\frac{z_{12}}{2\lambda _{12}}\mathrm{ln}\left(\frac{z_{12}+\lambda _{12}}{z_{12}\lambda _{12}}\right)]\hfill \\ \hfill +\frac{2\alpha }{\pi }\left[1\frac{z_{34}}{2\lambda _{34}}\mathrm{ln}\left(\frac{z_{34}+\lambda _{34}}{z_{34}\lambda _{34}}\right)\right]\\ \hfill +\frac{\alpha }{\pi }\frac{z_{13}}{\lambda _{13}}\mathrm{ln}\left(\frac{z_{13}+\lambda _{13}}{z_{13}\lambda _{13}}\right)\frac{\alpha }{\pi }\frac{z_{14}}{\lambda _{14}}\mathrm{ln}\left(\frac{z_{14}+\lambda _{14}}{z_{14}\lambda _{14}}\right)\\ \hfill \frac{\alpha }{\pi }\frac{z_{23}}{\lambda _{23}}\mathrm{ln}\left(\frac{z_{23}+\lambda _{23}}{z_{23}\lambda _{23}}\right)+\frac{\alpha }{\pi }\frac{z_{24}}{\lambda _{24}}\mathrm{ln}\left(\frac{z_{24}+\lambda _{24}}{z_{24}\lambda _{24}}\right)\}.\end{array}$$
(39)
As will be seen shortly, the first two terms cancel the divergent part of the vertex correction while the last four cancel the divergent part of the 5–point correction.

### V.2 Virtual Correction

As mentioned above, the virtual correction arises from the interference between the tree–level and one–loop diagrams. If the full matrix element is

$$=_{\text{tree}}+_{\text{virt}}+𝒪(e^6),$$
(40)
then the squared matrix element to $`𝒪(\alpha ^3)`$ is

$`|_{\text{rad}}|^2`$ $`=|_{\text{tree}}|^2\left[1+{\displaystyle \frac{2\mathrm{Re}(_{\text{tree}}^{}_{\text{virt}})}{|_{\text{tree}}|^2}}\right],`$ (41)
$`=|_{\text{tree}}|^2(1+\delta _{\text{virt}}).`$ (42)
This defines $`\delta _{\text{virt}}`$,

$$\delta _{\text{virt}}=\frac{2\mathrm{Re}(_{\text{tree}}^{}_{\text{virt}})}{|_{\text{tree}}|^2}.$$
(43)
Therefore, we must compute the matrix element for each of the one–loop contributions.

#### V.2.1 Vacuum Polarization

The vacuum polarization process involves higher order corrections to the photon propagator and is a function of the square of the photon momentum, or the $`x`$ of that pair. There is one contribution for each photon propagator. One contribution is shown in Fig. 6.

The vacuum polarization diagram is IR finite but UV divergent. The divergence can be handled by renormalization of the photon wavefunction. The vacuum polarization matrix element can be written as the tree–level matrix element times the renormalized polarization insertion

$$_{\text{vp}}=_{\text{tree}}\underset{l}{}\mathrm{\Pi }_l(x_{ij}),$$
(44)
where the sum is over lepton species in the loop and the renormalized polarization insertion is

$$\begin{array}{c}\mathrm{\Pi }_l(x_{ij})=\frac{2\alpha }{\pi }_0^1𝑑zz(1z)\hfill \\ \hfill \times \mathrm{ln}\left[1z(1z)x_{ij}M^2/m_l^2iϵ\right],\end{array}$$
(45)
where $`m_l`$ is the mass of the lepton in the loop. The integration depends on the size of $`x_{ij}`$ compared to $`m_l^2/M^2`$, such that

$$\begin{array}{c}\mathrm{\Pi }_l(x_{ij})=\frac{\alpha }{3\pi }\{\frac{8}{3}\beta _{ij}^2\hfill \\ \hfill +\frac{\beta _{ij}}{2}(3\beta _{ij}^2)[\mathrm{ln}\left(\frac{1\beta _{ij}}{1+\beta _{ij}}\right)+i\pi ]\},\end{array}$$
(46)
for $`x_{ij}>4m_l^2/M^2`$, while

$$\begin{array}{c}\mathrm{\Pi }_l(x_{ij})=\frac{\alpha }{3\pi }\{\frac{8}{3}+\rho _{ij}^2\hfill \\ \hfill \frac{\rho _{ij}}{2}(3+\rho _{ij}^2)[\pi 2\mathrm{tan}^1\rho _{ij}]\},\end{array}$$
(47)
for $`x_{ij}<4m_l^2/M^2`$. The functions $`\rho `$ and $`\beta `$ are related to $`\lambda `$ as defined in Appendix A but are functions of the loop mass

$`\beta _{ij}`$ $`=\sqrt{14m_l^2/(x_{ij}M^2)},`$ (48a)
$`\rho _{ij}`$ $`=\sqrt{4m_l^2/(x_{ij}M^2)1}.`$ (48b)
The correction then is

$$\delta _{\text{vp}}=2\underset{g}{}\left[\underset{l_g}{}\mathrm{Re}\mathrm{\Pi }_{l_g}(x_{ij})\right],$$
(49)
where the first sum is over the number of vacuum polarization graphs and the second sum is over the possible lepton species in the loop.

#### V.2.2 Vertex Function

The vertex function involves higher order corrections to the QED vertex and is a function of the momenta of the pair. One contribution is shown in Fig. 7.

The vertex correction contains both UV and IR divergences. We will also include the self–energy correction to the lepton lines which also are UV and IR divergent. Both UV divergences will be handled simultaneously by renormalization of the electromagnetic coupling and the lepton wavefunction while the IR divergent part will cancel the IR divergence in the soft bremsstrahlung correction. The matrix element corresponding to one of the diagrams is

$$_{\text{vc}}=_{\text{tree}}𝒱(x_{ij},y_{ij}),$$
(50)
where

$$𝒱(x_{ij},y_{ij})=F_1(x_{ij})+F_2(x_{ij})\left[\frac{2}{2\lambda _{ij}^2+y_{ij}^2}\right].$$
(51)
$`F_1`$ and $`F_2`$ are vertex form factors defined by

$`F_1(`$ $`x_{ij})={\displaystyle \frac{\alpha }{\pi }}(\{1{\displaystyle \frac{z_{ij}}{2\lambda _{ij}}}[\mathrm{ln}\left({\displaystyle \frac{z_{ij}+\lambda _{ij}}{z_{ij}\lambda _{ij}}}\right)2i\pi ]\}\mathrm{ln}{\displaystyle \frac{m_i}{\mathrm{\Lambda }}}`$
$`1+{\displaystyle \frac{1+2\lambda _{ij}^2}{8\lambda _{ij}}}\left[\mathrm{ln}\left({\displaystyle \frac{z_{ij}+\lambda _{ij}}{z_{ij}\lambda _{ij}}}\right)2i\pi \right]`$
$`{\displaystyle \frac{z_{ij}}{\lambda _{ij}}}\left[\mathrm{Li}_2\left({\displaystyle \frac{2\lambda _{ij}}{1+\lambda _{ij}}}\right)+{\displaystyle \frac{1}{16}}\mathrm{ln}^2\left({\displaystyle \frac{z_{ij}+\lambda _{ij}}{z_{ij}\lambda _{ij}}}\right){\displaystyle \frac{\pi ^2}{2}}\right]`$
$`+{\displaystyle \frac{i\pi }{2}}{\displaystyle \frac{z_{ij}}{\lambda _{ij}}}\mathrm{ln}\left({\displaystyle \frac{2\lambda _{ij}^2}{1z_{ij}}}\right)),`$ (52)
$`F_2(`$ $`x_{ij})={\displaystyle \frac{\alpha }{\pi }}{\displaystyle \frac{1z_{ij}}{4\lambda _{ij}}}[\mathrm{ln}\left({\displaystyle \frac{z_{ij}+\lambda _{ij}}{z_{ij}\lambda _{ij}}}\right)2i\pi ],`$ (53)
where $`z_{ij}`$ and $`\lambda _{ij}`$ are defined in Appendix A and $`\mathrm{\Lambda }`$ is the photon mass.

The correction is then just

$$\delta _{\text{vc}}=2\underset{g}{}\mathrm{Re}𝒱(x_{ij},y_{ij}),$$
(54)
where the sum over $`g`$ is over the number of vertex function graphs.

The IR divergent part is contained within $`F_1`$. Twice the real part of the $`\mathrm{ln}\mathrm{\Lambda }`$ term is exactly what is necessary to cancel the first two divergent terms (one for each vertex correction graph) in the soft brem contribution. The final four divergent soft brem terms will have to wait for the last piece of the virtual correction, the 5–point diagram.

#### V.2.3 5–Point Diagram

There are four distinct 5–point diagrams that contribute to the direct process, plus four more if there is an exchange process. The diagram shown in Fig. 8 contains a photon exchanged between $`l_1^{}`$ and $`l_2^+`$. The matrix element for that graph is

$$\begin{array}{c}_{\text{5p}}^1=\frac{d^4t}{(2\pi )^4}\frac{2}{M}\left[_Pϵ_{\mu \nu \rho \sigma }+_S(g_{\mu \rho }g_{\nu \sigma }g_{\mu \sigma }g_{\nu \rho })\right](p_{12}+t)^\mu (p_{34}t)^\rho \hfill \\ \hfill \times \left[\frac{ig^{\delta \eta }}{t^2\lambda ^2+iϵ}\right]\left[\frac{ig^{\alpha \nu }}{(p_{12}+t)^2\lambda ^2+iϵ}\right]\left[\frac{ig^{\beta \sigma }}{(p_{34}t)^2\lambda ^2+iϵ}\right]\\ \hfill \times \overline{u}(p_2)(ie\gamma _\delta )\left(\frac{i}{\text{/}t+\text{/}p_2m_1+iϵ}\right)(ie\gamma _\alpha )v(p_1)\overline{u}(p_4)(ie\gamma _\beta )\left(\frac{i}{\text{/}t\text{/}p_3m_2+iϵ}\right)(ie\gamma _\eta )v(p_3),\end{array}$$
(55)
where $`t`$ is the loop momentum and $`p_1`$, $`p_2`$, $`p_3`$, and $`p_4`$ are the momenta of $`l_1^+`$, $`l_1^{}`$, $`l_2^+`$, and $`l_2^{}`$, respectively. This can be re–expressed as

$`_{\text{5p}}^1`$ $`={\displaystyle \frac{2ie^4}{M}}(_Pϵ_{\mu \rho \nu \sigma }\{p_{12}^\mu p_{34}^\rho [I_{50}A^{\nu \sigma }+I_{51}^\alpha B_\alpha ^{\nu \sigma }+I_{52}^{\alpha \beta }C_{\alpha \beta }^{\nu \sigma }]+p_5^\mu [I_{51}^\rho A^{\nu \sigma }+I_{52}^{\rho \alpha }B_\alpha ^{\nu \sigma }+I_{53}^{\rho \alpha \beta }C_{\alpha \beta }^{\nu \sigma }]\}`$
$`+_S\{[(p_{12}p_{34})g^{\mu \nu }p_5^\mu p_5^\nu ][I_{50}A_{\mu \nu }+I_{51}^\alpha B_{\mu \nu \alpha }+I_{52}^{\alpha \beta }C_{\mu \nu \alpha \beta }]+(p_{34}p_{12})_\rho g_{\mu \nu }[I_{51}^\rho A^{\mu \nu }+I_{52}^{\rho \alpha }B_\alpha ^{\mu \nu }+I_{53}^{\rho \alpha \beta }C_{\alpha \beta }^{\mu \nu }]`$
$`g_{\mu \nu }[(I_{52})_\rho ^\rho A^{\mu \nu }+(I_{53})_{\alpha \rho }^\rho B^{\alpha \mu \nu }+(I_{54})_{\alpha \beta \rho }^\rho C^{\alpha \beta \mu \nu }]\}),`$ (56)
where $`A`$, $`B`$, and $`C`$ are combinations of spinors and gamma matrices and $`p_5=p_1+p_2+p_3+p_4`$. The factors of $`I_i`$ are integrals over the loop momentum. There are three basic integral forms from which all the other may be obtained. The notation for the integrals has the following meaning: the first digit in the subscript refers to the number of denominators and the second refers the number of powers of the loop momentum appearing in the integral. The 5–point integrals are defined in Appendix D as a function of four 4–vector arguments $`k_1`$, $`k_2`$, $`k_3`$, $`k_4`$.

For the diagram shown in Fig. 8, the arguments of the 5–point integral functions should take on the values

$`k_1`$ $`=p_2,`$ $`k_2`$ $`=p_1,`$ (57a)
$`k_3`$ $`=p_5,`$ $`k_4`$ $`=p_4,`$ (57b)
so we can write

$$I_{50}=I_{50}(p_2,p_1,p_5,p_4)$$
(58)
for the scalar integral, and analogous expressions for the higher-rank tensor integrals. The spinor terms for this diagram are

$`A^{\mu \nu }`$ $`=4(p_2p_3)\overline{u}(p_2)\gamma ^\mu v(p_1)\overline{u}(p_4)\gamma ^\nu v(p_3),`$ (59a)
$`B^{\mu \nu \rho }`$ $`=2[\overline{u}(p_2)\gamma ^\nu v(p_1)\overline{u}(p_4)\gamma ^\rho \gamma ^\mu \text{/}p_2v(p_3)`$
$`\overline{u}(p_2)\text{/}p_3\gamma ^\mu \gamma ^\nu v(p_1)\overline{u}(p_4)\gamma ^\rho v(p_3)],`$ (59b)
$`C^{\mu \nu \rho \sigma }`$ $`=\overline{u}(p_2)\gamma ^\eta \gamma ^\mu \gamma ^\rho v(p_1)\overline{u}(p_4)\gamma ^\sigma \gamma ^\nu \gamma _\eta v(p_3).`$ (59c)
The spinor terms for the diagram containing a photon exchanged between $`l_1^{}`$ and $`l_2^{}`$ are

$`A^{\mu \nu }`$ $`=4(p_2p_4)\overline{u}(p_2)\gamma ^\mu v(p_1)\overline{u}(p_4)\gamma ^\nu v(p_3),`$ (60a)
$`B^{\mu \nu \rho }`$ $`=2[\overline{u}(p_2)\text{/}p_4\gamma ^\mu \gamma ^\nu v(p_1)\overline{u}(p_4)\gamma ^\rho v(p_3)`$
$`\overline{u}(p_2)\gamma ^\nu v(p_1)\overline{u}(p_4)\text{/}p_2\gamma ^\mu \gamma ^\rho v(p_3)],`$ (60b)
$`C^{\mu \nu \rho \sigma }`$ $`=\overline{u}(p_2)\gamma ^\eta \gamma ^\mu \gamma ^\rho v(p_1)\overline{u}(p_4)\gamma _\eta \gamma ^\nu \gamma ^\sigma v(p_3),`$ (60c)
and the scalar integral for this diagram is

$$I_{50}=I_{50}(p_2,p_1,p_5,p_3).$$
(61)
The spinor terms for the diagram containing a photon exchanged between $`l_1^+`$ and $`l_2^+`$ are

$`A^{\mu \nu }`$ $`=4(p_1p_3)\overline{u}(p_2)\gamma ^\mu v(p_1)\overline{u}(p_4)\gamma ^\nu v(p_3),`$ (62a)
$`B^{\mu \nu \rho }`$ $`=2[\overline{u}(p_2)\gamma ^\nu \gamma ^\mu \text{/}p_3v(p_1)\overline{u}(p_4)\gamma ^\rho v(p_3)`$
$`\overline{u}(p_2)\gamma ^\nu v(p_1)\overline{u}(p_4)\gamma ^\rho \gamma ^\mu \text{/}p_1v(p_3)],`$ (62b)
$`C^{\mu \nu \rho \sigma }`$ $`=\overline{u}(p_2)\gamma ^\rho \gamma ^\mu \gamma ^\eta v(p_1)\overline{u}(p_4)\gamma ^\sigma \gamma ^\nu \gamma _\eta v(p_3),`$ (62c)
and the scalar integral for this diagram is

$$I_{50}=I_{50}(p_1,p_2,p_5,p_4).$$
(63)
The spinor terms for the diagram containing a photon exchanged between $`l_1^+`$ and $`l_2^{}`$ are

$`A^{\mu \nu }`$ $`=4(p_1p_4)\overline{u}(p_2)\gamma ^\mu v(p_1)\overline{u}(p_4)\gamma ^\nu v(p_3),`$ (64a)
$`B^{\mu \nu \rho }`$ $`=2[\overline{u}(p_2)\gamma ^\nu v(p_1)\overline{u}(p_4)\text{/}p_1\gamma ^\mu \gamma ^\rho v(p_3)`$
$`\overline{u}(p_2)\gamma ^\nu \gamma ^\mu \text{/}p_4v(p_1)\overline{u}(p_4)\gamma ^\rho v(p_3)],`$ (64b)
$`C^{\mu \nu \rho \sigma }`$ $`=\overline{u}(p_2)\gamma ^\rho \gamma ^\mu \gamma ^\eta v(p_1)\overline{u}(p_4)\gamma _\eta \gamma ^\nu \gamma ^\sigma v(p_3),`$ (64c)
and the scalar integral for this diagram is

$$I_{50}=I_{50}(p_1,p_2,p_5,p_3).$$
(65)
The tensors $`A`$, $`B`$, and $`C`$ are computed for a given helicity combination and combined with the integrals to yield the matrix element for that helicity state. The correction then involves a sum over the sixteen possible final states

$$\delta _{\text{5p}}=\frac{_{\lambda =1}^{16}2\mathrm{Re}\left[_{\text{tree}}^{}(\lambda )_{\text{5p}}(\lambda )\right]}{_{\lambda =1}^{16}\left|_{\text{tree}}(\lambda )\right|^2},$$
(66)
where $`\lambda `$ here refers to the helicity state and

$`_{\text{tree}}(\lambda )`$ $`={\displaystyle \underset{g}{}}_{\text{tree}}^g(\lambda ),`$ (67)
$`_{\text{5p}}(\lambda )`$ $`={\displaystyle \underset{g}{}}_{\text{5p}}^g(\lambda ).`$ (68)
The sums here are over the number of graphs for each process.

The IR divergent part of the 5–point correction is most easily isolated by looking at the 5–point matrix element in the IR limit. All terms involving tensor integrals vanish leaving only the $`I_{50}`$ term. The divergent part of $`I_{50}`$ is due to the two divergent box integrals, $`I_{40}^{(3)}`$ and $`I_{40}^{(4)}`$. The relevant terms in the scalar 5–point function for Fig. 8, in terms of the divergent 3–point function, are

$`I_{50}^{IR}`$ $`={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{_jS_{3j}^1}{p_{34}^2M_2^2}}+{\displaystyle \frac{_jS_{4j}^1}{p_{12}^2M_1^2}}\right)I_{IR},`$
$`=I_{IR}/\left[(p_{12}^2M_1^2)(p_{34}^2M_2^2)\right].`$ (69)
Extracting the $`\mathrm{ln}\mathrm{\Lambda }`$ piece of $`I_{IR}`$, $`I_{50}^{IR}`$ can be written as

$$\begin{array}{c}I_{50}^{IR}=\frac{i}{16\pi ^2\lambda _{23}p_{23}^2(p_{12}^2M_1^2)(p_{34}^2M_2^2)}\hfill \\ \hfill \times \left[\mathrm{ln}\left(\frac{z_{23}+\lambda _{23}}{z_{23}\lambda _{23}}\right)2i\pi \right]\mathrm{ln}\frac{\mathrm{\Lambda }}{\mu },\end{array}$$
(70)
where $`\mu `$ is a kinematic function with dimensions of mass which is independent of $`\mathrm{\Lambda }`$. The divergent part of the 5–point matrix element, dropping the finite term involving $`\mu `$, is proportional to the tree–level matrix element

$$_{\text{5p}}^{IR}=_{\text{tree}}\frac{\alpha }{2\pi }\frac{z_{23}}{\lambda _{23}}\left[\mathrm{ln}\left(\frac{z_{23}+\lambda _{23}}{z_{23}\lambda _{23}}\right)2i\pi \right]\mathrm{ln}\mathrm{\Lambda }.$$
(71)
The IR divergent part of the 5–point correction coming from all four diagrams is

$$\delta _{\text{5p}}^{IR}=\underset{g}{}s_{ij}\frac{\alpha }{\pi }\frac{z_{ij}}{\lambda _{ij}}\mathrm{ln}\left(\frac{z_{ij}+\lambda _{ij}}{z_{ij}\lambda _{ij}}\right)\mathrm{ln}\mathrm{\Lambda },$$
(72)
where $`s_{ij}`$ is the product of the sign of the charges of $`p_i`$ and $`p_j`$ and the sum is over the four diagrams. Again, it can be seen that this is the necessary form to cancel the remaining four divergent brem terms.

## VI MC Simulation Results

The inclusion of the radiative corrections impacts both the differential rate and the total rate. The total correction to the differential $`\pi ^0e^+e^{}e^+e^{}`$ rate is shown in Fig. 9.

The average size of the correction factor for the four different modes is shown in Table 5.

The total rate for the combined 4–lepton plus photon process is independent of IR cutoff. Table 6 summarizes the tree–level rate and the rate for the combined, cutoff independent process, divided into the rate including all radiation and the rate including only soft radiation ($`x_{4e}>0.95`$), all with $`f_P=f_S=1`$.

It is the last column which should most accurately predict the observed non–radiative 4–lepton rate. It is seen that the non–radiative rate is smaller than the tree–level rate for both 4–electron modes while it is larger for the modes with muons.

The probability of radiation can now be computed as the ratio of the radiative rate to the combined rate. Table 7 lists the probability of radiating a photon ($`x_{4e}<0.9985`$) along with the probability of radiating a hard photon ($`x_{4e}<0.95`$) for each of the four modes.

The probability is highest for $`K_Le^+e^{}e^+e^{}(\gamma )`$ where the $`x`$ values can be the smallest. The probabilities for $`K_Le^+e^{}\mu ^+\mu ^{}(\gamma )`$ are slightly less than half of what they are for the four electron mode, and in $`K_L\mu ^+\mu ^{}\mu ^+\mu ^{}(\gamma )`$, there is very little radiation.

The effect of the radiative corrections on the differential rate can be observed in the distributions of the five phase space variables. The statistics in the following plots reflect the amount of CPU time dedicated to each mode. While the calculation of the radiative corrections is CPU–intensive, it is actually the generation of the radiative decays that takes the most time.

For the modes with identical leptons, it is useful to adopt a method of pairing the electrons with the positrons in order to study the dilepton mass distributions. We choose to use the pairing for which the product of $`x`$’s is minimized. It is this pairing that will contribute the most to the matrix element in general. Therefore, $`x_a`$ and $`x_b`$ are the $`x`$’s belonging to this pairing, with the additional requirement that $`x_a<x_b`$. In addition, $`y_a`$ is the $`y`$ variable defined in the $`a`$–pair CM, and $`y_b`$ is the same quantity in the $`b`$–pair CM. And lastly, $`\varphi _{ab}`$ is the angle between the planes of the $`a`$–pair and $`b`$–pair in the overall CM.

The first variable that we will look at is $`x`$ which is modified by both the existence of a form factor and the inclusion of the radiative corrections. In all cases we set $`\beta =0`$ in the DIP form factor model. Figs. 10 and 11 show the distribution of $`x_a`$ and $`x_b`$, respectively, for $`\pi ^0e^+e^{}e^+e^{}`$ events. The plot on the left compares the distribution using the tree–level matrix element with no form factor ($`\alpha =0`$) to that using the same matrix element but with $`\alpha =1.0`$. The plot on the right compares the distribution using the tree–level matrix element with no form factor to that using the radiatively–corrected matrix element also with no form factor.

We have provided a linear fit to the ratio over some reasonable range on a scale appropriate for comparing the two effects. For the form factor comparisons, the dependence should be primarily linear. This is not the case for the radiative corrections in general. The $`\chi ^2`$ per degree of freedom is included as a measure of the linearity. It can be seen that the form factor has a much smaller effect on the $`x`$ distribution than the radiative corrections do. This is not too surprising since the range of accessible $`q^2`$ values for the $`\pi ^0`$ decay is relatively small in addition to being far from our assumed $`\rho `$ pole.

For the kaon modes, we observe that the form factor has a much larger effect on the $`x`$ distribution than the radiative corrections do. Figs. 12 and 13 show the distribution of $`x_a`$ and $`x_b`$, respectively, for $`K_Le^+e^{}e^+e^{}`$ events. The plot on the left compares the distribution using the tree–level matrix element with no form factor ($`\alpha =0`$) to that using the same matrix element but with $`\alpha =1.5`$. The plot on the right compares the distribution using the tree–level matrix element with no form factor to that using the radiatively–corrected matrix element also with no form factor.

The roll off at high $`x`$ in plot (c) of Fig. 13 is due to presence of the exchange diagram in this mode. Figs. 14 and 15 show the same distributions for $`K_Le^+e^{}\mu ^+\mu ^{}`$ events. Here there are no pairing ambiguities and we plot $`x_{ee}`$ and $`x_{\mu \mu }`$.

It can be seen that there is no roll off in plot (c) of Fig. 15, and furthermore, a small quadratic dependence is observable. While the $`x`$ of the $`ee`$ pair is slightly modified by the radiative corrections, the $`x`$ of the $`\mu \mu `$ pair does not change at all. This is as expected for the massive muons.

Fig. 16 shows the distribution of $`y_a`$ and $`y_b`$ for the tree level differential rate and the radiatively corrected differential rate for $`K_Le^+e^{}e^+e^{}`$ events.

The effect here is quite small. Since $`y`$ is a measure of the energy asymmetry of the lepton pair, it is seen that the radiative corrections tend to make the pairs slightly more asymmetric on average.

The effect on the $`\varphi `$ distribution is due entirely to the 5–point diagram. Fig. 17 shows a comparison of the distribution of $`\varphi _{ab}`$ generated with the tree–level matrix element to the same distribution generated with the radiative corrections, for $`K_L\mu ^+\mu ^{}\mu ^+\mu ^{}`$ events.

The enhancement at $`\varphi _{ab}=\pi `$ and the corresponding depletion at $`\varphi _{ab}=0=2\pi `$ can be understood in terms of Coulomb interaction between the final state particles. The configuration at $`\varphi =0`$ has all leptons in a plane with the same sign particles near each other. The effect is only observable in the $`K_L\mu ^+\mu ^{}\mu ^+\mu ^{}`$ decay where the leptons in each pair are usually well separated.

## VII Conclusions

The main conclusion that can be drawn from these distributions is that the radiative corrections are extremely important for extracting a form factor in the $`\pi ^0e^+e^{}e^+e^{}`$ mode. For the kaon modes, the form factor has a larger impact on the $`x`$ distribution and the modification of the distribution due to the radiative corrections is less important. The only published result for the $`K_Le^+e^{}e^+e^{}`$ mode Alavi-Harati et al. (2001a) quotes $`\alpha _{DIP}=1.1\pm 0.6(\text{stat})`$. Fig. 13 shows that at most the radiative corrections would change the slope by $`0.22`$, which while significant, is smaller than the current experimental error. Likewise, for the $`K_Le^+e^{}\mu ^+\mu ^{}`$ mode, the latest resultHamm (2002) based on the invariant mass shape is $`\alpha _{DIP}=4.53_{2.70}^{+1.81}`$. The present experimental error is again larger than the impact of the radiative corrections on the mass distribution.

As for the extraction of the mixing angle $`\zeta `$ from the observed $`\varphi `$ distribution, the radiative corrections can be safely neglected at present.

The two publications above quote an integrated rate (normalized to the two–photon rate) of $`(6.24\pm 0.34)10^5`$ for $`K_Le^+e^{}e^+e^{}`$ and $`(4.51\pm 0.42)10^6`$ for $`K_Le^+e^{}\mu ^+\mu ^{}`$, where the errors are purely statistical. These results are in good agreement with our predictions when both the radiative corrections and a form factor with $`\alpha =1.5`$ are included. For $`K_Le^+e^{}e^+e^{}`$, the two effects offset and the net result is an increase of just less than $`2\%`$ over the tree–level rate with no form factor. In $`K_Le^+e^{}\mu ^+\mu ^{}`$, the form factor is the dominant effect.

###### Acknowledgements.

This work was supported by Department of Energy grant DE-FG03-95ER40894 and by NSF/REU grant PHY-0097381. The authors also acknowledge useful comments from Prof. L. M. Sehgal.

## Appendix A Kinematics

The four particle final state can be kinematically described by considering subsystems containing only two particles. Consider the system composed of two particles with momenta $`p_i`$ and $`p_j`$ and total momentum $`p_{ij}=p_i+p_j`$ and mass squared $`m_{ij}^2=p_{ij}^2`$. We will define a dimensionless dot product of any two vectors $`p_i`$ and $`p_j`$ as

$$z_{ij}=2(p_ip_j)/p_{ij}^2=1x_ix_j,$$
(73)
where

$$x_i=p_i^2/p_{ij}^2,x_j=p_j^2/p_{ij}^2.$$
(74)
The energy and momentum of each particle in the two–particle CM frame are

$$E_i^{}=m_{ij}(1+\delta _{ij})/2,E_j^{}=m_{ij}(1\delta _{ij})/2,$$
(75)
$$p^{}=m_{ij}\lambda _{ij}/2,$$
(76)
where

$`\delta _{ij}`$ $`=x_ix_j,`$ (77)
$`\lambda _{ij}`$ $`=\sqrt{z_{ij}^2w_{ij}^2},`$ (78)
$`w_{ij}`$ $`=2\sqrt{x_ix_j}.`$ (79)
Occasionally we will use symbols like $`z_{i,jk}`$ whose meaning is interpreted as $`z(p_i,p_j+p_k)`$.

Now consider a three body system composed of momenta $`p_i`$, $`p_j`$, and $`p_k`$. There are two phase space variables needed to describe the system. The first one will be $`x_{ij}=p_{ij}^2/p_{ijk}^2`$. The other one is defined in the $`ij`$ CM frame as the cosine of the angle between the direction of particle $`i`$ and particle $`k`$,

$$\mathrm{cos}\theta _{ij,k}=\frac{2p_k(p_ip_j)p_{ijk}^2\delta _{ij}z_{ij,k}}{p_{ijk}^2\lambda _{ij}\lambda _{ij,k}},$$
(80)
A more convenient variable that will be used in place of $`\mathrm{cos}\theta _{ij,k}`$ is

$$y_{ij}=\lambda _{ij}\mathrm{cos}\theta _{ij,k},$$
(81)
where $`k`$ will always refer to the total momentum minus the $`ij`$ momentum.

Finally, the four body final state requires five phase space variables to uniquely describe it. We will use the $`x`$ and $`y`$ values for the two lepton pairs plus the angle between the normals of the planes defined by each pair in the overall CM frame. The first four variables are

$`x_{12}`$ $`=p_{12}^2/M^2,`$ (82)
$`x_{34}`$ $`=p_{34}^2/M^2,`$ (83)
$`y_{12}`$ $`={\displaystyle \frac{2p_{34}(p_1p_2)}{M^2\lambda }},`$ (84)
$`y_{34}`$ $`={\displaystyle \frac{2p_{12}(p_3p_4)}{M^2\lambda }},`$ (85)
where the second term in the numerator of the $`y`$’s vanishes since $`\delta _{12}=\delta _{34}=0`$. Any use of $`z`$, $`w`$, or $`\lambda `$ without subscripts will refer to the functions of $`x_{12}`$ and $`x_{34}`$, so $`\lambda =\lambda _{12,34}`$ for instance. The last phase space variable is defined as

$$\varphi =\mathrm{tan}^1(\mathrm{sin}\varphi /\mathrm{cos}\varphi ),$$
(86)
where

$`\mathrm{sin}\varphi `$ $`={\displaystyle \frac{16ϵ_{\mu \nu \rho \sigma }p_1^\mu p_2^\nu p_3^\rho p_4^\sigma }{M^4\lambda w\sqrt{(\lambda _{12}^2y_{12}^2)(\lambda _{34}^2y_{34}^2)}}},`$ (87)
$`\mathrm{cos}\varphi `$ $`={\displaystyle \frac{M^2zy_{12}y_{34}2(p_1p_2)(p_3p_4)}{M^2w\sqrt{(\lambda _{12}^2y_{12}^2)(\lambda _{34}^2y_{34}^2)}}}.`$ (88)
The angle $`\varphi `$ is defined so that at $`\varphi =0`$ the two pairs lie in a plane with the like–sign particles adjacent to each other. The orientation of $`\varphi =\pi `$ again has both pairs in a plane, but with the opposite signed particles adjacent.

The general expression for the phase space integral is

$$d^8\mathrm{\Phi }=\frac{1}{(2\pi )^8}\frac{d^3p_1d^3p_2d^3p_3d^3p_4}{16E_{p_1}E_{p_2}E_{p_3}E_{p_4}}\delta ^4(Pp_1p_2p_3p_4).$$
(89)
Upon integrating out the $`\delta `$–functions, integrating over the Euler angles, and changing variables to those listed above, the phase space reduces to

$$d^5\mathrm{\Phi }=𝒮\frac{M^4}{2^{14}\pi ^6}\lambda dx_{12}dx_{34}dy_{12}dy_{34}d\varphi ,$$
(90)
where the factor $`𝒮`$ is a symmetry factor which is required for modes containing identical particles in the final state. The double Dalitz modes with identical particles contain two sets, thus requiring two factors of $`1/2`$. So $`𝒮=1/4`$ if the final state contains identical particles, and $`𝒮=1`$ otherwise.

When there are identical leptons in the final state, the amplitudes for the exchange diagrams have the same algebraic form as for the non–exchange diagrams except that the kinematic variables $`x_{12}`$, $`x_{34}`$, $`y_{12}`$, $`y_{34}`$, and $`\varphi `$ are replaced by $`x_{14}`$, $`x_{23}`$, $`y_{14}`$, $`y_{23}`$, and $`\varphi _{14,23}`$. These exchange variables will in general be functions of all five of the non–exchange variables. As is turns out, we only need explicit representations for $`x_{14}`$ and $`x_{23}`$. These are given by

$`x_{14}`$ $`=[1\lambda _{12}^2x_{12}\lambda _{34}^2x_{34}\lambda (y_{12}y_{34})zy_{12}y_{34}`$
$`+w\sqrt{(\lambda _{12}^2y_{12}^2)(\lambda _{34}^2y_{34}^2)}\mathrm{cos}\varphi ]/4,`$ (91)
$`x_{23}`$ $`=[1\lambda _{12}^2x_{12}\lambda _{34}^2x_{34}+\lambda (y_{12}y_{34})zy_{12}y_{34}`$
$`+w\sqrt{(\lambda _{12}^2y_{12}^2)(\lambda _{34}^2y_{34}^2)}\mathrm{cos}\varphi ]/4.`$ (92)
## Appendix B Meson–$`\gamma \gamma `$ Couplings

In this section we will work out the explicit form of the two–photon couplings, allowing for photons of arbitrary mass, using the polarization vectors in the helicity basis. The general form of the coupling is

$$_{\lambda _1\lambda _2}=H_{\mu \nu \rho \sigma }k_1^\mu ϵ_{\lambda _1}^\nu k_2^\rho ϵ_{\lambda _2}^\sigma ,$$
(93)
where $`H`$ is either

$`H_{\mu \nu \rho \sigma }^P`$ $`={\displaystyle \frac{2}{M}}_Pϵ_{\mu \nu \rho \sigma },`$ (94a)
$`H_{\mu \nu \rho \sigma }^S`$ $`={\displaystyle \frac{2}{M}}_S(g_{\mu \rho }g_{\nu \sigma }g_{\mu \sigma }g_{\nu \rho }).`$ (94b)
The three polarization vectors for a massive photon in the helicity basis are chosen to be

$`ϵ^+(\pm \widehat{𝒛})`$ $`=(0,1,\pm i,0)/\sqrt{2},`$ (95a)
$`ϵ^{}(\pm \widehat{𝒛})`$ $`=(0,1,i,0)/\sqrt{2},`$ (95b)
$`ϵ^0(\pm \widehat{𝒛})`$ $`=(k,0,0,\pm E)/\sqrt{k^2},`$ (95c)
for a photon traveling in the $`\pm \widehat{𝒛}`$ direction. With these polarization vectors, one finds three couplings for the scalar case

$$_{\lambda _1\lambda _2}^S=\{\begin{array}{cc}M_Sz\hfill & ,\lambda _1=\lambda _2=+,\hfill \\ M_Sz\hfill & ,\lambda _1=\lambda _2=,\hfill \\ +M_Sw\hfill & ,\lambda _1=\lambda _2=0,\hfill \end{array}$$
(96)
where $`z`$ and $`w`$ are defined in Appendix A. The longitudinal contribution vanishes for the pseudoscalar case, and one finds only two couplings

$$_{\lambda _1\lambda _2}^P=\{\begin{array}{cc}+iM_P\lambda \hfill & ,\lambda _1=\lambda _2=+,\hfill \\ +iM_P\lambda \hfill & ,\lambda _1=\lambda _2=,\hfill \\ 0\hfill & ,\lambda _1=\lambda _2=0,\hfill \end{array}$$
(97)
where $`\lambda `$ is also defined in Appendix A. There are three interesting differences between the scalar and the pseudoscalar couplings. First, assuming that $`\delta =0`$, there is a relative phase between them. Additionally, the transverse couplings differ in the kinematic factor. And lastly, there is the additional scalar coupling due to the contribution from longitudinally polarized photons. Where as the transverse couplings go like $`\lambda `$ or $`z`$, both of which are $`𝒪(1)`$ on average, the longitudinal coupling goes like $`w`$, which is $`𝒪(x)`$, making its contribution less significant.

## Appendix C Double Dalitz Interference

The interference between the tree–level direct and exchange contributions for modes with identical leptons is a sum of three terms

$`2\mathrm{Re}`$ $`(_1^{}_2)={\displaystyle \frac{2^3\pi ^2\alpha ^2\stackrel{~}{g}^2}{M^2x_{12}x_{34}x_{14}x_{23}}}`$
$`\times \{Af_P(x_{12},x_{34})f_P(x_{14},x_{23})\mathrm{cos}^2\zeta `$
$`+B[f_P(x_{12},x_{34})f_S(x_{14},x_{23})`$
$`+f_P(x_{14},x_{23})f_S(x_{12},x_{34})]\mathrm{sin}\zeta \mathrm{cos}\zeta `$
$`+Cf_S(x_{12},x_{34})f_S(x_{14},x_{23})\mathrm{sin}^2\zeta \},`$ (98)
where

$`A`$ $`=\lambda ^2\{2\mathrm{\Xi }^2+8\eta ^4+\mathrm{\Xi }[8\eta ^2z(y_{12}+y_{34})^2]`$
$`w^2(1+y_{12}y_{34})(2y_{12}^2y_{34}^2)`$
$`+4\eta ^2(y_{12}+y_{34})(x_{12}y_{12}+x_{34}y_{34})\},`$ (99)
$`B`$ $`=2\lambda [\eta ^2(x_{12}+x_{34}x_{14}x_{23})`$
$`+4(x_{12}x_{34}x_{14}x_{23})]\mathrm{\Xi }\mathrm{tan}\varphi ,`$ (100)
$`C`$ $`=\mathrm{\Xi }^3z\mathrm{\Xi }^2[6\eta ^2zw^2(2z^2+w^2)y_{12}y_{34}]`$
$`\mathrm{\Xi }[z(z^22w^2)+z(z^2+2w^2)y_{12}^2y_{34}^2+2z^3y_{12}y_{34}`$
$`2\eta ^2(z^2+w^2)(1+3y_{12}y_{34})+8\eta ^4z]`$
$`+w^2z^2y_{12}^3y_{34}^3w^2z^2y_{12}y_{34}+w^2(3z^22w^2)y_{12}^2y_{34}^2`$
$`2w^2(z^2w^2)(y_{12}^2+y_{34}^2)+w^2(z^22w^2)`$
$`2\eta ^2[3w^2zy_{12}^2y_{34}^22(z^2z^3w^2)y_{12}y_{34}`$
$`+2z(z^23w^2/2)2(z^2+w^2)]`$
$`8\eta ^4[(zz^2)w^2y_{12}y_{34}],`$ (101)
where $`\mathrm{\Xi }=w\sqrt{(\lambda _{12}^2y_{12}^2)(\lambda _{34}^2y_{34}^2)}\mathrm{cos}\varphi `$ and $`\eta ^2=4m^2/M^2`$. The exchange variables $`x_{14}`$ and $`x_{23}`$ are defined in Appendix A in terms of the five non–exchange phase space variables. The term proportional to $`\mathrm{cos}^2\zeta `$ results from interference between a pseudoscalar coupling in both the direct and exchange graphs, while the one proportional to $`\mathrm{sin}^2\zeta `$ is due to scalar couplings in both graphs, and the one proportional to $`\mathrm{sin}\zeta \mathrm{cos}\zeta `$ is due to a pseudoscalar coupling in one graph and a scalar coupling in the other.

## Appendix D 5–Point Function

The matrix element for the 5–point diagram is composed of tensor integrals with five propagators in the denominator. One can express tensor integrals in terms of lower rank tensor integrals with the same number of propagators and lower rank tensors with fewer propagators van Oldenborgh and Vermaseren (1990). In the end, every tensor integral can be decomposed into scalar 2–, 3–, 4–, and 5–point functions. The scalar 5–point function is not independent and can itself be expressed in terms of scalar 4–point functions.

This appendix will outline our procedure for first reducing the tensor integrals to scalar integrals, and then computing the scalar integrals in closed form.

#### D.0.1 Tensor Integrals

Begin by defining

$`I_{50}(k_1,k_2,k_3,k_4)`$ $`={\displaystyle \frac{d^4t}{(2\pi )^4}\frac{1}{N_1N_2N_3N_4N_5}},`$ $`I_{52}^{\mu \nu }(k_1,k_2,k_3,k_4)`$ $`={\displaystyle \frac{d^4t}{(2\pi )^4}\frac{t^\mu t^\nu }{N_1N_2N_3N_4N_5}},`$ (102)
$`I_{51}^\mu (k_1,k_2,k_3,k_4)`$ $`={\displaystyle \frac{d^4t}{(2\pi )^4}\frac{t^\mu }{N_1N_2N_3N_4N_5}},`$ $`I_{53}^{\mu \nu \rho }(k_1,k_2,k_3,k_4)`$ $`={\displaystyle \frac{d^4t}{(2\pi )^4}\frac{t^\mu t^\nu t^\rho }{N_1N_2N_3N_4N_5}},`$ (103)
where

$$N_1=t^2\mu _1^2,N_2=(t+k_1)^2\mu _2^2,N_3=(t+k_1+k_2)^2\mu _3^2,$$
(104)
$$N_4=(t+k_1+k_2+k_3)^2\mu _4^2,N_5=(t+k_1+k_2+k_3+k_4)^2\mu _5^2,$$
(105)
where $`\mu _i`$ is an internal mass and the $`k_i`$ are external momenta.

The original reduction scheme of Ref. Passarino and Veltman (1979), while theoretically sound, suffers from uncontrollable numerical inaccuracies. To avoid this problem, we follow the procedure suggested in Ref. van Oldenborgh and Vermaseren (1990), and use a reduction scheme based on the Schouten identity which utilizes Gram determinants to express any tensor integral as a sum of integrals, one with the same number of propagators and the rest with one less propagator, and all with the power of the loop momentum reduced by one. The identity has the following form

$`t^\mu ϵ^{k_1k_2k_3k_4}`$ $`=(tk_1)ϵ^{\mu k_2k_3k_4}+(tk_2)ϵ^{k_1\mu k_3k_4}+(tk_3)ϵ^{k_1k_2\mu k_4}+(tk_4)ϵ^{k_1k_2k_3\mu },`$ (106)
$`=v^\mu {\displaystyle \frac{1}{2}}\left(N_1ϵ^{\mu k_2k_3k_4}N_2ϵ^{\mu (k_1+k_2)k_3k_4}+N_3ϵ^{\mu k_1(k_2+k_3)k_4}N_4ϵ^{\mu k_1k_2(k_3+k_4)}+N_5ϵ^{\mu k_1k_2k_3}\right),`$ (107)
where

$$v^\mu =(sk_1)ϵ^{\mu k_2k_3k_4}(sk_2)ϵ^{\mu k_1k_3k_4}+(sk_3)ϵ^{\mu k_1k_2k_4}(sk_4)ϵ^{\mu k_1k_2k_3},$$
(108)
and $`s^\mu `$ is defined in terms of its dot products

$`sk_1`$ $`={\displaystyle \frac{1}{2}}\left[\mu _2^2\mu _1^2k_1^2\right],`$ $`sk_3`$ $`={\displaystyle \frac{1}{2}}\left[\mu _4^2\mu _1^2(k_1+k_2+k_3)^2\right]sk_1sk_2,`$ (109)
$`sk_2`$ $`={\displaystyle \frac{1}{2}}\left[\mu _3^2\mu _1^2(k_1+k_2)^2\right]sk_1,`$ $`sk_4`$ $`={\displaystyle \frac{1}{2}}\left[\mu _5^2\mu _1^2(k_1+k_2+k_3+k_4)^2\right]sk_1sk_2sk_3.`$ (110)
The notation $`ϵ^{k_1k_2k_3k_4}`$ is shorthand for $`ϵ^{\mu \nu \rho \sigma }k_{1\mu }k_{2\nu }k_{3\rho }k_{4\sigma }`$.

The reduction then proceeds as follows:

$`I_{51}^\mu `$ $`={\displaystyle \frac{1}{ϵ^{k_1k_2k_3k_4}}}\{v^\mu I_{50}{\displaystyle \frac{1}{2}}[ϵ^{\mu k_2k_3k_4}I_{40}^{(1)}ϵ^{\mu (k_1+k_2)k_3k_4}I_{40}^{(2)}+ϵ^{\mu k_1(k_2+k_3)k_4}I_{40}^{(3)}ϵ^{\mu k_1k_2(k_3+k_4)}I_{40}^{(4)}`$
$`+ϵ^{\mu k_1k_2k_3}I_{40}^{(5)}]\},`$ (111)
$`I_{52}^{\mu \nu }`$ $`={\displaystyle \frac{1}{ϵ^{k_1k_2k_3k_4}}}\{v^\mu I_{51}^\nu {\displaystyle \frac{1}{2}}[ϵ^{\mu k_2k_3k_4}(I_{41}^{(1)\nu }k_1^\nu I_{40}^{(1)})ϵ^{\mu (k_1+k_2)k_3k_4}I_{41}^{(2)\nu }+ϵ^{\mu k_1(k_2+k_3)k_4}I_{41}^{(3)\nu }`$
$`ϵ^{\mu k_1k_2(k_3+k_4)}I_{41}^{(4)\nu }+ϵ^{\mu k_1k_2k_3}I_{41}^{(5)\nu }]\},`$ (112)
$`I_{53}^{\mu \nu \rho }`$ $`={\displaystyle \frac{1}{ϵ^{k_1k_2k_3k_4}}}\{v^\mu I_{52}^{\nu \rho }{\displaystyle \frac{1}{2}}[ϵ^{\mu k_2k_3k_4}(I_{42}^{(1)\nu \rho }k_1^\nu I_{41}^{(1)\rho }k_1^\rho I_{41}^{(1)\nu }+k_1^\nu k_1^\rho I_{40}^{(1)})ϵ^{\mu (k_1+k_2)k_3k_4}I_{42}^{(2)\nu \rho }`$
$`+ϵ^{\mu k_1(k_2+k_3)k_4}I_{42}^{(3)\nu \rho }ϵ^{\mu k_1k_2(k_3+k_4)}I_{42}^{(4)\nu \rho }+ϵ^{\mu k_1k_2k_3}I_{42}^{(5)\nu \rho }]\},`$ (113)
where for any $`I_{40}=\frac{d^4t}{(2\pi )^4}\frac{1}{N_1N_2N_3N_4}`$,

$`I_{41}^\mu `$ $`={\displaystyle \frac{1}{2\delta _{k_1k_2k_3}^{k_1k_2k_3}}}\left\{\delta _{k_1k_2k_3}^{s\alpha \beta }\delta _{k_1k_2k_3}^{\mu \alpha \beta }I_{40}\left[\delta _{k_1k_2k_3}^{\mu k_2k_3}I_{30}^{(1)}\delta _{k_1k_2k_3}^{\mu (k_1+k_2)k_3}I_{30}^{(2)}+\delta _{k_1k_2k_3}^{\mu k_1(k_2+k_3)}I_{30}^{(3)}\delta _{k_1k_2k_3}^{\mu k_1k_2}I_{30}^{(4)}\right]\right\},`$ (114)
$`I_{42}^{\mu \nu }`$ $`={\displaystyle \frac{1}{2\delta _{k_1k_2k_3}^{k_1k_2k_3}}}\{\delta _{k_1k_2k_3}^{s\alpha \beta }\delta _{k_1k_2k_3}^{\mu \alpha \beta }I_{41}^\nu [\delta _{k_1k_2k_3}^{\mu k_2k_3}(I_{31}^\nu k_1^\nu I_{30}^{(1)})\delta _{k_1k_2k_3}^{\mu (k_1+k_2)k_3}I_{31}^{(2)\nu }`$
$`+\delta _{k_1k_2k_3}^{\mu k_1(k_2+k_3)}I_{31}^{(3)\nu }\delta _{k_1k_2k_3}^{\mu k_1k_2}I_{31}^{(4)\nu }]`$
$`+{\displaystyle \frac{2ϵ^{k_1k_2k_3\mu }ϵ^{k_1k_2k_3\nu }}{\delta _{k_1k_2k_3}^{k_1k_2k_3}}}[\delta _{sk_1k_2k_3}^{sk_1k_2k_3}I_{40}+{\displaystyle \frac{1}{2}}(\delta _{k_1k_2k_3}^{(s+k_1)k_2k_3}I_{30}^{(1)}\delta _{k_1k_2k_3}^{s(k_1+k_2)k_3}I_{30}^{(2)}`$
$`+\delta _{k_1k_2k_3}^{sk_1(k_2+k_3)}I_{30}^{(3)}\delta _{k_1k_2k_3}^{sk_1k_2}I_{30}^{(4)})]\},`$ (115)
and for any $`I_{30}=\frac{d^4t}{(2\pi )^4}\frac{1}{N_1N_2N_3}`$

$$I_{31}^\mu =\frac{1}{\delta _{k_1k_2}^{k_1k_2}}\left\{\delta _{k_1k_2}^{s\alpha }\delta _{k_1k_2}^{\mu \alpha }I_{30}\frac{1}{2}\left[\delta _{k_1k_2}^{\mu k_2}I_{20}^{(1)}\delta _{k_1k_2}^{\mu (k_1+k_2)}I_{20}^{(2)}+\delta _{k_1k_2}^{\mu k_1}I_{20}^{(3)}\right]\right\}.$$
(116)
The reduction notation has the following meaning: $`I_{40}^{(i)}`$ is the 4–point function obtained from $`I_{50}`$ by dropping the $`i`$th propagator, $`I_{30}^{(j)}`$ is the 3–point function obtained from its corresponding 4–point function by dropping the $`j`$th propagator, and so on. The Gram determinants that appear in the reduction are kinematic functions which are defined as

$`\delta _{q_1q_2}^{k_1k_2}`$ $`=\left|\begin{array}{cc}q_1k_1& q_1k_2\\ q_2k_1& q_2k_2\end{array}\right|,`$ $`\delta _{q_1q_2q_3}^{k_1k_2k_3}`$ $`=\left|\begin{array}{ccc}q_1k_1& q_1k_2& q_1k_3\\ q_2k_1& q_2k_2& q_2k_3\\ q_3k_1& q_3k_2& q_3k_3\end{array}\right|,`$ $`\delta _{q_1q_2q_3q_4}^{k_1k_2k_3k_4}`$ $`=\left|\begin{array}{cccc}q_1k_1& q_1k_2& q_1k_3& q_1k_4\\ q_2k_1& q_2k_2& q_2k_3& q_2k_4\\ q_3k_1& q_3k_2& q_3k_3& q_3k_4\\ q_4k_1& q_4k_2& q_4k_3& q_4k_4\end{array}\right|.`$ (117)
The traces that appear in Eq. 56 can also be reduced

$`(I_{52})_\mu ^\mu `$ $`=I_{40}^{(1)},`$ (118a)
$`(I_{53})_\mu ^{\alpha \mu }`$ $`=I_{41}^{(1)\alpha }k_1^\alpha I_{40}^{(1)},`$ (118b)
$`(I_{54})_\mu ^{\alpha \beta \mu }`$ $`=I_{42}^{(1)\alpha \beta }k_1^\alpha I_{41}^{(1)\beta }k_1^\beta I_{41}^{(1)\alpha }+k_1^\alpha k_1^\beta I_{40}^{(1)}.`$ (118c)
While this procedure is generally much more reliable then that of Ref. Passarino and Veltman (1979), there are still problems that occur when $`v^\mu `$ as defined in Eq. (108) is not an independent combination of the final state momenta. This happens when all the momenta in the parent particle’s rest frame lie in a plane. Even though these configurations form a subspace of zero measure in the final state phase space, finite numerical precision dictates that they will be generated with non–zero probability by the MC program. In this case, we use the identity

$$t^\mu \delta _{k_1k_2k_3}^{k_1k_2k_3}=u^\mu (ut)+\frac{1}{2}\left[\delta _{k_1k_2k_3}^{s\alpha \beta }\delta _{k_1k_2k_3}^{\mu \alpha \beta }N_1\delta _{k_1k_2k_3}^{\mu k_2k_3}+N_2\delta _{k_1k_2k_3}^{\mu k_3k_4}N_3\delta _{k_1k_2k_3}^{\mu k_4k_1}+N_4\delta _{k_1k_2k_3}^{\mu k_1k_2}\right],$$
(119)
where $`u^\mu =ϵ^{k_1k_2k_3\mu }`$. The first term on the right vanishes upon integration, allowing $`I_{51}^\mu `$ to be written as

$$I_{51}^\mu =\frac{1}{2\delta _{k_1k_2k_3}^{k_1k_2k_3}}\left\{\delta _{k_1k_2k_3}^{s\alpha \beta }\delta _{k_1k_2k_3}^{\mu \alpha \beta }I_{50}\left[\delta _{k_1k_2k_3}^{\mu k_2k_3}I_{40}^{(1)}\delta _{k_1k_2k_3}^{\mu k_3k_4}I_{40}^{(2)}+\delta _{k_1k_2k_3}^{\mu k_4k_1}I_{40}^{(3)}\delta _{k_1k_2k_3}^{\mu k_1k_2}I_{40}^{(4)}\right]\right\}.$$
(120)
The other tensor integrals can be expanded in a similar manner. The same problem can arise in the reduction of the 4–point tensor integrals if the three momenta $`k_1`$, $`k_2`$, and $`k_3`$ are linearly dependent, in which case this same procedure is reproduced at one lower rank.

In these degenerate cases we have a choice between numerical inaccuracies resulting from antisymmetric invariants such as $`ϵ^{k_1k_2k_3k_4}`$ being very small or inaccuracies resulting from assuming exact linear dependence. To decide which approximation is better, we do the calculation of the tensor integrals both ways and check whether identities such as those in Eqs. (118) are satisfied. More than $`99\%`$ of the time, one of the two methods yields good agreement for all these “trace checks”.

#### D.0.2 Scalar Integrals

The most general 5–point function we will need to consider is

$`I_{50}=\mu ^{2ϵ}{\displaystyle }{\displaystyle \frac{d^Dt}{(2\pi )^D}}\{`$ $`{\displaystyle \frac{1}{[t^2\mu _1^2+iϵ][(t+k_1)^2\mu _2^2+iϵ][(t+k_1+k_2)^2\mu _3^2+iϵ]}}`$
$`\times {\displaystyle \frac{1}{[(t+k_1+k_2+k_3)^2\mu _4^2+iϵ][(t+k_1+k_2+k_3+k_4)^2\mu _5^2+iϵ]}}\}.`$ (121)
In the case at hand $`\mu _1=0`$, $`\mu _2=m_1`$, $`\mu _3=M_1`$, $`\mu _4=M_2`$, and $`\mu _5=m_2`$ where the $`m`$’s are lepton masses and the $`M`$’s are boson masses. In addition, $`k_1=p_2`$, $`k_2=p_1`$, $`k_3=p_5`$, and $`k_4=p_4`$, where the external momenta satisfy the relation $`p_1+p_2+p_3+p_4=p_5`$. The diagram representing this integral is shown in Fig. 18.

We have allowed the boson propagators to have arbitrary masses $`M_1`$ and $`M_2`$ in order to include a form factor in the calculation of the 5–point function. This is necessary since the form factor becomes a function of the loop momentum. We will use a generalized DIP form factor which, with the appropriate choice of coefficients, can reproduce both the DIP and the BMS form factor models. The generalized DIP form factor is

$$f(x_1,x_2)=1+\underset{i}{}\alpha _i\left(\frac{x_1}{x_1M_i^2/M^2}+\frac{x_2}{x_2M_i^2/M^2}\right)+\underset{ij}{}\beta _{ij}\frac{x_1x_2}{(x_1M_i^2/M^2)(x_2M_j^2/M^2)},$$
(122)
where $`M`$ is the total mass and the sum is over propagator masses $`M_i`$. The diagram containing two photon propagators and a form factor is then replaced by a sum of $`n^2`$ diagrams of four different types, one containing two photon propagators, two containing one photon and one massive boson propagator, and one containing two massive boson propagators,

$$\frac{f(x_1,x_2)}{x_1x_2}=\frac{1}{x_1x_2}+\underset{i}{}\alpha _i\left(\frac{1}{x_1(x_2M_i^2/M^2)}+\frac{1}{x_2(x_1M_i^2/M^2)}\right)+\underset{ij}{}\beta _{ij}\frac{1}{(x_1M_i^2/M^2)(x_2M_j^2/M^2)}.$$
(123)
The general form of Eq. (121) permits the evaluation of all four of these contributions. To use a specific model, the parameters $`\alpha _i`$ and $`\beta _{ij}`$ must be adjusted. For the simple DIP model, only a $`\rho `$ meson term is included with $`\alpha _1=\alpha _{DIP}`$ and $`\beta _{11}=\beta _{DIP}`$. The DIP model requires the inclusion of four 5–point diagrams involving all combinations of photons and $`\rho `$ mesons. The BMS model is more complicated, requiring 25 different diagrams. To simplify this we have let $`M_\omega =M_\rho `$, which reduces the number of diagrams to 16. The values of the generalized parameters, in terms of $`\alpha _K^{}`$ are

$`\alpha _1`$ $`=1+{\displaystyle \frac{10C\alpha _K^{}}{9(1M_\rho ^2/M_K^{}^2)}},`$ (124a)
$`\alpha _2`$ $`={\displaystyle \frac{2C\alpha _K^{}}{9}}\left({\displaystyle \frac{1}{1M_\varphi ^2/M_K^{}^2}}{\displaystyle \frac{1}{1M_\rho ^2/M_K^{}^2}}\right),`$ (124b)
$`\alpha _3`$ $`={\displaystyle \frac{2C\alpha _K^{}}{1M_\varphi ^2/M_K^{}^2}},`$ (124c)
and $`\beta _{ij}=\alpha _i\alpha _j`$. A flat form factor is obtained by setting all of the $`\alpha _i=0`$ and $`\beta _{ij}=0`$.

We will write the 5–point function in Eq. (121) as a sum of 4–point functions using the following relationship between $`n`$–point functions and $`(n1)`$–point functions Bern et al. (1993),

$$I_n=\frac{1}{2}\left[\underset{i=1}{\overset{n}{}}c_iI_{n1}^{(i)}+(n5+2ϵ)c_0I_n^{D=62ϵ}\right],$$
(125)
where

$`c_i`$ $`={\displaystyle \underset{j=1}{\overset{n}{}}}𝒮_{ij}^1,`$ $`c_0`$ $`={\displaystyle \underset{i=1}{\overset{n}{}}}c_i={\displaystyle \underset{i,j=1}{\overset{n}{}}}𝒮_{ij}^1,`$ $`S_{ij}`$ $`=(\mu _i^2+\mu _j^2k_{ij}^2)/2,`$ (126)
and

$`k_{ii}`$ $`=0,`$ $`k_{ij}`$ $`=k_i+k_{i+1}+\mathrm{}+k_{j1},`$ for $`i<j.`$ (127)
We will use $`\mu `$ to refer to the propagator mass when the distinction between vector bosons and leptons is irrelevant. In the case at hand $`n=5`$ so the second term in Eq. (125) is $`𝒪(ϵ)`$, and since the 5–point function in $`D=62ϵ`$ dimensions is finite in the limit $`ϵ0`$, the scalar 5–point function can be written as a sum of five scalar 4–point functions

$$I_{50}=\frac{1}{2}\underset{i=1}{\overset{5}{}}\left(\underset{j=1}{\overset{5}{}}S_{ij}^1\right)I_{40}^{(i)}+𝒪(ϵ).$$
(128)
##### The 4–Point Function

Two of the five 4–point functions contain IR divergences due to the presence of the 3–point functions where both of the vector boson propagators have been removed. Therefore there are two distinct 4–point functions that we will need. The first has one zero mass propagator

$$\begin{array}{c}I_{40}(k_1,k_2,k_3,\mu _1,\mu _2,\mu _3)=\hfill \\ \hfill \mu ^{2ϵ}\frac{d^Dt}{(2\pi )^D}\frac{1}{[t^2+iϵ][(t+k_1)^2\mu _1^2+iϵ][(t+k_1+k_2)^2\mu _2^2+iϵ][(t+k_1+k_2+k_3)^2\mu _3^2+iϵ]},\end{array}$$
(129)
and the other has two zero mass propagators and two lepton propagators

$$\begin{array}{c}I_{40}^{}(k_1,k_2,k_3,m_1,m_2)=\hfill \\ \hfill \mu ^{2ϵ}\frac{d^Dt}{(2\pi )^D}\frac{1}{[t^2+iϵ][(t+k_1)^2+iϵ][(t+k_1+k_2)^2m_1^2+iϵ][(t+k_1+k_2+k_3)^2m_2^2+iϵ]},\end{array}$$
(130)
where in $`I_{40}^{}`$, $`(k_1+k_2)^2=m_1^2`$ and $`(k_1+k_2+k_3)^2=m_2^2`$.

In order to use Eq. (129) for $`I_{40}^{(1)}`$ where all four propagators have non–zero masses, and to extract the divergent part of $`I_{40}^{(3)}`$ and $`I_{40}^{(4)}`$, we make use of the following propagator identity

$$\frac{1}{[(q+p)^2\mu _1^2][(q+p+k)^2\mu _2^2]}=\frac{\alpha }{[(q+p)^2\mu _1^2](q+p+\alpha k)^2}+\frac{1\alpha }{[(q+p+k)^2\mu _2^2](q+p+\alpha k)^2},$$
(131)
where $`\alpha `$ is chosen to be the positive root of the equation

$$\alpha (1\alpha )k^2(1\alpha )\mu _1^2\alpha \mu _2^2=0.$$
(132)
This identity allows us to write the five 4–point functions as

$`I_{40}^{(1)}`$ $`=\alpha _1I_{40}(\alpha _1p_{23},p_1,p_5,m_1,M_1,M_2)+(1\alpha _1)I_{40}((1\alpha _1)p_{23},p_4p_5,p_5,m_2,M_1,M_2),`$ (133a)
$`I_{40}^{(2)}`$ $`=I_{40}(p_{12},p_5,p_4,M_1,M_2,m_2),`$ (133b)
$`I_{40}^{(3)}`$ $`=\alpha _3I_{40}^{}(\alpha _3p_{34},p_2,p_{23},m_1,m_2)+(1\alpha _3)I_{40}((1\alpha _3)p_{34},p_1p_5,p_{23},M_2,m_1,m_2),`$ (133c)
$`I_{40}^{(4)}`$ $`=\alpha _4I_{40}^{}(\alpha _4p_{12},p_3,p_{23},m_2,m_1)+(1\alpha _4)I_{40}((1\alpha _4)p_{12},p_1,p_{23},M_1,m_1,m_2),`$ (133d)
$`I_{40}^{(5)}`$ $`=I_{40}(p_2,p_1,p_5,m_1,M_1,M_2),`$ (133e)
where

$`\alpha _1`$ $`=(1+\delta _{23}+\lambda _{23})/2,`$ $`\alpha _3`$ $`=1M_2^2/p_{34}^2,`$ $`\alpha _4`$ $`=1M_1^2/p_{12}^2.`$ (134)
The finite 4–point function $`I_{40}`$ defined in Eq. (129) can be expressed in closed form as a sum of 36 dilogarithms. We will define it in terms of the function

$`𝒥(A,B)`$ $`={\displaystyle _0^1}{\displaystyle \frac{dz}{zA}}\left[\mathrm{ln}(zB\pm iϵ)\mathrm{ln}(AB\pm iϵ)\right],`$
$`=\mathrm{Li}_2\left({\displaystyle \frac{A}{AB\pm iϵ}}\right)\mathrm{Li}_2\left({\displaystyle \frac{A1}{AB\pm iϵ}}\right).`$ (135)
For arbitrary complex arguments, $`A`$ and $`B`$, the integration would also produce additional logarithms with prefactors which depend on the relative difference between the signs of the imaginary parts of $`A`$ and $`B`$ t’Hooft and Veltman (1979), however if $`A`$ is real these additional terms vanish. In the case at hand, $`A`$ will always be real and we will only need the dilogarithms. In terms of these new functions,

$`I_{40}`$ $`(k_1,k_2,k_3,\mu _1,\mu _2,\mu _3)={\displaystyle \frac{i}{16\pi ^2\gamma (\eta ^+\eta ^{})}}`$
$`\times [𝒥(\beta +\eta ^+,{\displaystyle \frac{ed+iϵ}{k}})𝒥(\beta +\eta ^{},{\displaystyle \frac{ed+iϵ}{k}})𝒥({\displaystyle \frac{\eta ^+}{1\beta }},{\displaystyle \frac{d+iϵ}{e+k}})`$
$`+𝒥({\displaystyle \frac{\eta ^{}}{1\beta }},{\displaystyle \frac{d+iϵ}{e+k}})+𝒥({\displaystyle \frac{\eta ^+}{\beta }},{\displaystyle \frac{d+iϵ}{e}})𝒥({\displaystyle \frac{\eta ^{}}{\beta }},{\displaystyle \frac{d+iϵ}{e}})`$
$`𝒥(\beta +\eta ^+,z_1^++iϵ)𝒥(\beta +\eta ^+,z_1^{}iϵ)+𝒥(\beta +\eta ^{},z_1^++iϵ)+𝒥(\beta +\eta ^{},z_1^{}iϵ)`$
$`+𝒥({\displaystyle \frac{\eta ^+}{1\beta }},z_2^++iϵ)+𝒥({\displaystyle \frac{\eta ^+}{1\beta }},z_2^{}iϵ)𝒥({\displaystyle \frac{\eta ^{}}{1\beta }},z_2^++iϵ)𝒥({\displaystyle \frac{\eta ^{}}{1\beta }},z_2^{}iϵ)`$
$`𝒥({\displaystyle \frac{\eta ^+}{\beta }},z_3^++iϵ)𝒥({\displaystyle \frac{\eta ^+}{\beta }},z_3^{}iϵ)+𝒥({\displaystyle \frac{\eta ^{}}{\beta }},z_3^++iϵ)+𝒥({\displaystyle \frac{\eta ^{}}{\beta }},z_3^{}iϵ)],`$ (136)
where $`\beta `$ is either root of the equation $`g\beta ^2+j\beta +b=0`$, $`\eta ^\pm `$ are the roots of the equation

$$[egjk\beta gk]\eta ^2+[ehckdj2\beta dg]\eta +[aecd+\beta (akdh)]=0,$$
(137)
and $`\gamma =egjk\beta gk`$. The quantities $`z_i^\pm `$ are the roots of the following equations

$`0`$ $`=gz_1^2+(h+j+k)z_1+(a+b+c+d+e),`$ (138a)
$`0`$ $`=(b+g+j)z_2^2+(c+e+h+k)z_2+(a+d),`$ (138b)
$`0`$ $`=bz_3^2+(c+e)z_3+(a+d).`$ (138c)
The lower case variables are combinations of the elements of the relevant $`4\times 4`$ matrix $`S_{ij}`$ defined in Eq. (126),

$`a`$ $`=S_{33}+S_{44}2S_{34},`$ $`f`$ $`=S_{44},`$ (139a)
$`b`$ $`=S_{22}+S_{33}2S_{23},`$ $`g`$ $`=S_{11}+S_{22}2S_{12},`$ (139b)
$`c`$ $`=2(S_{23}S_{24}S_{33}+S_{34}),`$ $`h`$ $`=2(S_{13}S_{14}S_{23}+S_{24}),`$ (139c)
$`d`$ $`=2(S_{34}S_{44}),`$ $`j`$ $`=2(S_{12}S_{13}S_{22}+S_{23}),`$ (139d)
$`e`$ $`=2(S_{24}S_{34}),`$ $`k`$ $`=2(S_{14}S_{24}).`$ (139e)
The divergent 4–point function $`I_{40}^{}`$ defined in Eq. (130) can also be written in closed from. The divergent part is just the divergent 3–point function

$`I_{IR}`$ $`={\displaystyle \frac{d^4t}{(2\pi )^4}\frac{1}{[t^2\mathrm{\Lambda }^2+iϵ][(tp_2)^2m_1^2+iϵ][(t+p_3)^2m_2^2+iϵ]}},`$
$`={\displaystyle \frac{i}{32\pi ^2p_{23}^2\lambda _{23}}}\{[\mathrm{ln}\left({\displaystyle \frac{z_{23}+\lambda _{23}}{z_{23}\lambda _{23}}}\right)2i\pi ]\mathrm{ln}{\displaystyle \frac{p_{23}^2}{\mathrm{\Lambda }^2}}2\pi ^2+2\mathrm{Li}_2\left({\displaystyle \frac{2\lambda _{23}}{1+\delta _{23}+\lambda _{23}}}\right)+2\mathrm{Li}_2\left({\displaystyle \frac{2\lambda _{23}}{1\delta _{23}+\lambda _{23}}}\right)`$
$`4i\pi \mathrm{ln}\lambda _{23}\mathrm{ln}\left({\displaystyle \frac{1+\delta _{23}\lambda _{23}}{2\lambda _{23}}}\right)\mathrm{ln}\left({\displaystyle \frac{1+\delta _{23}+\lambda _{23}}{2}}\right)\mathrm{ln}\left({\displaystyle \frac{1\delta _{23}\lambda _{23}}{2\lambda _{23}}}\right)\mathrm{ln}\left({\displaystyle \frac{1\delta _{23}+\lambda _{23}}{2}}\right)`$
$`+\mathrm{ln}^2\left({\displaystyle \frac{1\delta _{23}+\lambda _{23}}{2\lambda _{23}}}\right)+\mathrm{ln}^2\left({\displaystyle \frac{1+\delta _{23}+\lambda _{23}}{2\lambda _{23}}}\right)+\mathrm{ln}\lambda _{23}\mathrm{ln}\left({\displaystyle \frac{1\delta _{23}+\lambda _{23}}{2\lambda _{23}}}\right)+\mathrm{ln}\lambda _{23}\mathrm{ln}\left({\displaystyle \frac{1+\delta _{23}+\lambda _{23}}{2\lambda _{23}}}\right)`$
$`+{\displaystyle \frac{1}{2}}[\mathrm{ln}^2\left({\displaystyle \frac{1+\delta _{23}+\lambda _{23}}{2}}\right)+\mathrm{ln}^2\left({\displaystyle \frac{1\delta _{23}+\lambda _{23}}{2}}\right)\mathrm{ln}^2\left({\displaystyle \frac{1+\delta _{23}\lambda _{23}}{2}}\right)\mathrm{ln}^2\left({\displaystyle \frac{1\delta _{23}\lambda _{23}}{2}}\right)]\},`$ (140)
The full expression is then

$`I_{40}^{}(k_1,k_2,k_3,m_1,m_2)`$ $`={\displaystyle \frac{1}{k_1^2}}(I_{IR}+{\displaystyle \frac{i}{16\pi ^2\lambda _{23}k_3^2}}\{\mathrm{ln}\left({\displaystyle \frac{2k_1k_3}{k_1^2}}\right)[\mathrm{ln}\left({\displaystyle \frac{z_{23}\lambda _{23}}{z_{23}+\lambda _{23}}}\right)+2i\pi ]`$
$`+\mathrm{Li}_2\left({\displaystyle \frac{1A}{BA+iϵ}}\right)\mathrm{Li}_2\left({\displaystyle \frac{1A}{CAiϵ}}\right)\mathrm{Li}_2\left({\displaystyle \frac{A}{BA+iϵ}}\right)+\mathrm{Li}_2\left({\displaystyle \frac{A}{CAiϵ}}\right)`$
$`+\mathrm{ln}(1Aiϵ)\left[\mathrm{ln}(1Biϵ)\mathrm{ln}(1C+iϵ)\mathrm{ln}(ABiϵ)+\mathrm{ln}(AC+iϵ)\right]`$
$`\mathrm{ln}(Aiϵ)[\mathrm{ln}(Biϵ)\mathrm{ln}(C+iϵ)\mathrm{ln}(ABiϵ)+\mathrm{ln}(AC+iϵ)]\}),`$ (141)
where

$`A`$ $`=(k_1^2/2+k_1k_2)/(k_1k_3),`$ $`B`$ $`=(1+\delta _{23}+\lambda _{23})/2,`$ $`C`$ $`=(1+\delta _{23}\lambda _{23})/2.`$ (142)
It is worthwhile to consider the special case of two photon propagators, that is $`M_1=M_2=0`$. In that case the 4–point functions simplify considerably. We will write

$`I_{40}^{}(k_1,k_2,k_3,\mu _1,\mu _2)`$ $`={\displaystyle \frac{d^4t}{(2\pi )^4}\frac{1}{[t^2\mu _1^2+iϵ][(t+k_1)^2\mu _2^2+iϵ][(t+k_1+k_2)^2+iϵ][(t+k_1+k_2+k_3)^2+iϵ]}},`$
$`={\displaystyle \frac{i}{\sqrt{\mathrm{\Delta }}}}\{[\mathrm{Li}_2\left({\displaystyle \frac{v^+}{v^+r_1}}\right)\mathrm{Li}_2\left({\displaystyle \frac{v^+1}{v^+r_1}}\right)\mathrm{Li}_2\left({\displaystyle \frac{v^+1}{v^+r_2}}\right)\mathrm{Li}_2\left({\displaystyle \frac{v^+}{v^+\stackrel{~}{v}^+iϵ}}\right)`$
$`+\mathrm{Li}_2\left({\displaystyle \frac{v^+1}{v^+\stackrel{~}{v}^+iϵ}}\right)\mathrm{Li}_2\left({\displaystyle \frac{v^+}{v^+\stackrel{~}{v}^{}+iϵ}}\right)+\mathrm{Li}_2\left({\displaystyle \frac{v^+1}{v^+\stackrel{~}{v}^{}+iϵ}}\right)]`$
$`[\mathrm{Li}_2\left({\displaystyle \frac{v^{}}{v^{}r_1}}\right)\mathrm{Li}_2\left({\displaystyle \frac{v^{}1}{v^{}r_1}}\right)\mathrm{Li}_2\left({\displaystyle \frac{v^{}1}{v^{}r_2}}\right)\mathrm{Li}_2\left({\displaystyle \frac{v^{}}{v^{}\stackrel{~}{v}^+iϵ}}\right)`$
$`+\mathrm{Li}_2\left({\displaystyle \frac{v^{}1}{v^{}\stackrel{~}{v}^+iϵ}}\right)\mathrm{Li}_2\left({\displaystyle \frac{v^{}}{v^{}\stackrel{~}{v}^{}+iϵ}}\right)+\mathrm{Li}_2\left({\displaystyle \frac{v^{}1}{v^{}\stackrel{~}{v}^{}+iϵ}}\right)]\},`$ (143)
where $`\mathrm{\Delta }`$ is the discriminant and $`v^\pm `$ are the roots of the quadratic equation

$$[k(h+k)dg]v^2+[e(h+k)dj]vbd=0,$$
(144)
and $`\stackrel{~}{v}^\pm `$ are the roots of

$$g\stackrel{~}{v}^2+(h+j+k)\stackrel{~}{v}+b=0,$$
(145)
and

$`r_1`$ $`=(de+iϵ)/k,`$ $`r_2`$ $`=iϵ/(h+k).`$ (146)
Therefore, when $`M_1=M_2=0`$, the five 4–point functions are simply

$`I_{40}^{(1)}`$ $`=I_{40}^{}(p_{23},p_4,p_5,m_1,m_2),`$ (147a)
$`I_{40}^{(2)}`$ $`=I_{40}^{}(p_3,p_4,p_5,0,m_2),`$ (147b)
$`I_{40}^{(3)}`$ $`=I_{40}^{}(p_{34},p_2,p_{23},m_1,m_2),`$ (147c)
$`I_{40}^{(4)}`$ $`=I_{40}^{}(p_{12},p_3,p_{23},m_2,m_1),`$ (147d)
$`I_{40}^{(5)}`$ $`=I_{40}^{}(p_2,p_1,p_5,0,m_1).`$ (147e)
Also, in this case, $`I_{40}^{}`$ simplifies somewhat because $`A`$ defined in Eq. (142) becomes one. The first two dilogarithms in the second line of Eq. (141), along with the entire third line, vanish in this case. $`I_{40}^{}`$ then becomes

$`I_{40}^{}(k_1,k_2,k_3,m_1,m_2)`$ $`={\displaystyle \frac{1}{k_1^2}}(I_{IR}+{\displaystyle \frac{i}{16\pi ^2\lambda _{23}k_3^2}}\{\mathrm{ln}\left({\displaystyle \frac{2k_1k_3}{k_1^2}}\right)[\mathrm{ln}\left({\displaystyle \frac{z_{23}\lambda _{23}}{z_{23}+\lambda _{23}}}\right)+2i\pi ]`$
$`\mathrm{Li}_2\left({\displaystyle \frac{1}{B1+iϵ}}\right)+\mathrm{Li}_2\left({\displaystyle \frac{1}{C1iϵ}}\right)`$
$`+i\pi [\mathrm{ln}(Biϵ)\mathrm{ln}(C+iϵ)\mathrm{ln}(1Biϵ)+\mathrm{ln}(1C+iϵ)]\}),`$ (148)
where $`B`$ and $`C`$ are still given by Eq. (142).

##### The 3–Point Function

There are ten 3–point functions that are needed, all of which are finite except one, $`I_{30}^{(34)}=I_{IR}`$ defined in Eq. (D.0.2). The superscripts used in this section denote the two propagators that have been dropped from the original 5–point function to obtain the particular 3–point function. The finite 3–point functions can be generically written as

$`I_{30}(k_1,k_2,\mu _1,\mu _2,\mu _3)`$ $`={\displaystyle \frac{d^Dt}{(2\pi )^D}\frac{1}{[t^2\mu _1^2+iϵ][(t+k_1)^2\mu _2^2+iϵ][(t+k_1+k_2)^2\mu _3^3+iϵ]}},`$ (149)
$`={\displaystyle \frac{i}{16\pi ^2(c+2b\beta )}}\{[\mathrm{Li}_2\left({\displaystyle \frac{v_1}{v_1\stackrel{~}{v}_1^+iϵ}}\right)\mathrm{Li}_2\left({\displaystyle \frac{v_11}{v_1\stackrel{~}{v}_1^+iϵ}}\right)+\mathrm{Li}_2\left({\displaystyle \frac{v_1}{v_1\stackrel{~}{v}_1^{}+iϵ}}\right)`$
$`\mathrm{Li}_2\left({\displaystyle \frac{v_11}{v_1\stackrel{~}{v}_1^{}+iϵ}}\right)\mathrm{Li}_2\left({\displaystyle \frac{v_2}{v_2\stackrel{~}{v}_2^+iϵ}}\right)+\mathrm{Li}_2\left({\displaystyle \frac{v_21}{v_2\stackrel{~}{v}_2^+iϵ}}\right)`$
$`\mathrm{Li}_2\left({\displaystyle \frac{v_2}{v_2\stackrel{~}{v}_2^{}+iϵ}}\right)+\mathrm{Li}_2\left({\displaystyle \frac{v_21}{v_2\stackrel{~}{v}_2^{}+iϵ}}\right)+\mathrm{Li}_2\left({\displaystyle \frac{v_3}{v_3\stackrel{~}{v}_3^+iϵ}}\right)`$
$`\mathrm{Li}_2\left({\displaystyle \frac{v_31}{v_3\stackrel{~}{v}_3^+iϵ}}\right)+\mathrm{Li}_2\left({\displaystyle \frac{v_3}{v_3\stackrel{~}{v}_3^{}+iϵ}}\right)\mathrm{Li}_2\left({\displaystyle \frac{v_31}{v_3\stackrel{~}{v}_3^{}+iϵ}}\right)]\},`$ (150)
where the $`v_i`$ are

$`v_1`$ $`={\displaystyle \frac{2a+d+\beta (c+e)}{c+2b\beta }},`$ $`v_2`$ $`={\displaystyle \frac{d+e\beta }{(1\beta )(c+2b\beta )}},`$ $`v_3`$ $`={\displaystyle \frac{d+e\beta }{\beta (c+2b\beta )}},`$ (151)
and the $`\stackrel{~}{v}_i^\pm `$ are roots of the three quadratic equations

$`0`$ $`=b\stackrel{~}{v}_1^2+(c+e)\stackrel{~}{v}_1+(a+d+f),`$ (152a)
$`0`$ $`=(a+b+c)\stackrel{~}{v}_2^2+(d+e)\stackrel{~}{v}_2+f,`$ (152b)
$`0`$ $`=a\stackrel{~}{v}_3^2+d\stackrel{~}{v}_3+f,`$ (152c)
and $`\beta `$ is either root of the equation $`b\beta ^2+c\beta +a=0`$. The lower case letters are again combination of the relevant matrix $`S_{ij}`$

$`a`$ $`=S_{22}+S_{33}2S_{23},`$ $`c`$ $`=2(S_{12}S_{13}S_{22}+S_{23}),`$ $`e`$ $`=2(S_{13}S_{23}),`$ (153a)
$`b`$ $`=S_{11}+S_{22}2S_{12},`$ $`d`$ $`=2(S_{23}S_{33}),`$ $`f`$ $`=S_{33}.`$ (153b)
The ten 3–point functions that we need can then be expressed as

$`I_{30}^{(12)}`$ $`=I_{30}(p_4,p_5,m_2,M_2,M_1),`$ $`I_{30}^{(24)}`$ $`=I_{30}(p_3,p_5+p_4,0,m_2,M_1),`$ (154a)
$`I_{30}^{(13)}`$ $`=I_{30}(p_{23},p_4,m_1,m_2,M_2),`$ $`I_{30}^{(25)}`$ $`=I_{30}(p_{12},p_5,0,M_1,M_2),`$ (154b)
$`I_{30}^{(14)}`$ $`=I_{30}(p_{23},p_5+p_4,m_1,m_2,M_1),`$ $`I_{30}^{(34)}`$ $`=I_{IR},`$ (154c)
$`I_{30}^{(15)}`$ $`=I_{30}(p_1,p_5,m_1,M_1,M_2),`$ $`I_{30}^{(35)}`$ $`=I_{30}(p_2,p_5p_1,0,m_1,M_2),`$ (154d)
$`I_{30}^{(23)}`$ $`=I_{30}(p_3,p_4,0,m_2,M_2),`$ $`I_{30}^{(45)}`$ $`=I_{30}(p_2,p_1,0,m_1,M_1).`$ (154e)
##### The 2–Point Function

And finally, the general expression for the 2–point function is

$`I_{20}(k_1,\mu _1,\mu _2)`$ $`={\displaystyle \frac{d^Dt}{(2\pi )^D}\frac{1}{[t^2\mu _1^2+iϵ][(t+k_1)^2\mu _2^2+iϵ]}},`$ (155)
$`={\displaystyle \frac{i}{(4\pi )^2}}[\mathrm{\Gamma }(\epsilon )\left({\displaystyle \frac{4\pi \mu ^2}{k_1^2}}\right)^\epsilon +2(1v^+)\mathrm{ln}(1v^+iϵ)(1v^{})\mathrm{ln}(1v^{}+iϵ)`$
$`v^+\mathrm{ln}(v^+iϵ)v^{}\mathrm{ln}(v^{}+iϵ)],`$ (156)
where $`v^\pm `$are roots to the quadratic equation

$$k_1^2v^2+(\mu _1^2\mu _2^2k_1^2)v+\mu _2^2=0.$$
(157)
The UV divergent term containing $`\epsilon `$ cancels when the 2–point functions are combined to form the tensor integrals and can therefore be safely ignored. The ten 2–point functions that we require are then

$`I_{20}^{(123)}`$ $`=I_{20}(p_4,M_2,m_2),`$ $`I_{20}^{(145)}`$ $`=I_{20}(p_1,m_1,M_1),`$ (158a)
$`I_{20}^{(124)}`$ $`=I_{20}(p_5p_4,M_1,m_2),`$ $`I_{20}^{(234)}`$ $`=I_{20}(p_3,0,m_2),`$ (158b)
$`I_{20}^{(125)}`$ $`=I_{20}(p_5,M_1,M_2),`$ $`I_{20}^{(235)}`$ $`=I_{20}(p_{34},0,M_2),`$ (158c)
$`I_{20}^{(134)}`$ $`=I_{20}(p_{23},m_1,m_2),`$ $`I_{20}^{(245)}`$ $`=I_{20}(p_{12},0,M_1),`$ (158d)
$`I_{20}^{(135)}`$ $`=I_{20}(p_5p_1,m_1,M_2),`$ $`I_{20}^{(345)}`$ $`=I_{20}(p_2,0,m_1).`$ (158e)