
\subsection{NNLO corrections to the Bhabha scattering cross section \label{NNLO}}

Beyond the NLO corrections discussed in the previous Section, in recent years a
significant effort was devoted to the calculation of the perturbative corrections  to
the Bhabha process at NNLO in QED.

The calculation of the full NNLO corrections to the Bhabha scattering cross section 
requires three types of ingredients: {\it i)} the two-loop matrix elements for the 
$e^{+}e^{-} \to e^{+}e^{-}$ process; {\it ii)} the one-loop matrix elements for the 
$e^{+}e^{-} \to e^{+}e^{-} \gamma$ process, both in the case in which the additional 
photon is soft or hard; %-collinear\footnote{\Black{By hard-collinear we mean a hard photon which
%is collinear to a charged particle within a given acollinearity angle imposed by the
%experimental cuts.}}; 
 {\it iii)} the tree-level matrix elements for $e^{+}e^{-} \to
e^{+}e^{-}\gamma\gamma$, with two soft or two hard photons, or one soft and one
hard photon. Also the process $e^{+}e^{-} \to e^{+}e^{-} e^{+}e^{-}$, with %the second
one of the two $e^{+}e^{-}$ pairs remaining 
undetected, contributes to the Bhabha signature at NNLO. 
%, because data taking at meson factories typically requires the presence of at least two tracks in the detector. 
 %soft or hard-collinear.
%
Depending on the kinematics, other final states like, e.g., $e^{+}e^{-} \mu^{+}\mu^{-}$
or those with hadrons are also possible.

%Moreover, the fixed-order corrections should be matched with a parton shower generator
%that takes into account the logarithmically enhanced contributions of soft and
%collinear photons at all orders in perturbation theory (see Section~\ref{HO}).

The advent of new calculational techniques and a deeper understanding of the IR
structure of unbroken gauge theories, such as QED or QCD, made the calculation of
the complete set of two-loop QED corrections possible. The history of this
calculation will be presented in Section~\ref{virtual}.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%We consider then the one-loop  matrix elements with three particles in the final state.
Some remarks on the one-loop matrix elements with
three particles in the final state are in order now. 
The diagrams involving the emission of a soft photon are known and they were included
in the calculations of the two-loop matrix elements, in order to remove the IR soft
divergences. However, although the contributions due to a hard collinear  photon are
taken  into account in logarithmic accuracy by the MC generators, a full  calculation
of the diagrams involving a hard photon in a general phase-space configuration is still
missing. In Section~\ref{Penta}, we shall comment on the possible strategies which can
be  adopted in order to calculate these corrections.\footnote{As emphasised in Section \ref{TH} and 
Section \ref{CONCLUSIONS}, the complete calculation of this class of corrections became available 
\cite{Actis:2009uq} during the completion of the present work.}

As a general comment, it must be noticed that the fixed-order corrections 
calculated up to NNLO are taken into account at the LL, and, 
partially, next-to-leading-log (NLL) level in the most precise MC generators, which include, as 
will be discussed in Section 
\ref{MP} and Section \ref{MC}, the logarithmically enhanced contributions of soft and
collinear photons at all orders in perturbation theory.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Concerning the tree level graphs with four particles in the final state, the production of a 
soft $e^{+}e^{-}$ pair was considered in the literature by the authors of \cite{Arbuzov:1995vj} 
by following the evaluation of pair production \cite{Arbuzov:1995cn,Arbuzov:1995vi} within the calculation of the 
$O (\alpha^2 L)$ single-logarithmic accurate small-angle Bhabha cross section 
\cite{Arbuzov:1995qd}, 
and it is included in the two-loop calculation (see
Section~\ref{ElLoop}). New results on lepton and hadron pair corrections, which are at present 
approximately included in the available Bhabha codes, are presented in 
Section~\ref{4pfinal}.

%New results obtained during the workshop about lepton and
%hadron pair corrections, which are presently neglected into the available 
%Bhabha codes, are presented  in Section~\ref{4pfinal}.



%we can
%distinguish once more between the soft-photon contributions, which are included in the
%two-loop calculations,  and the contributions involving one hard or two hard photons,
%which are taken into  account in logarithmic accuracy by the %parton shower in the Monte
%MC generators, as discussed in Section 
%\ref{MP} and Section \ref{MC}. 
%The production of a soft $e^{+}e^{-}$ pair was also studied,
%but only with logarithmic accuracy and it is included in the two-loop calculation (see
%Section~\ref{ElLoop}). However, the production of a hard-collinear pair of $e^{+}e^{-}$
%is not included in the current calculation of the NNLO corrections.
%Further comments on final states with four particles will be given

%Finally, in section~\ref{Conclusions} we draw our conclusions.

\subsubsection{Virtual corrections for the $e^{+}e^{-} \to e^{+}e^{-}$ process \label{virtual}}

The calculation of the virtual two-loop QED corrections to the Bhabha scattering
differential cross section was carried out in the last 10 years.  This calculation was
made possible by an improvement of the techniques employed in  the evaluation of
multi-loop Feynman diagrams. An essential tool used to manage the calculation is the 
Laporta algorithm \cite{Laporta:1996mq,Laporta:2001dd,Tkachov:1981wb,Chetyrkin:1981qh}, which enables one to reduce a generic combination
of dimensionally-regularised scalar integrals to a combination of a small set of
independent integrals called the ``Master Integrals'' (MIs) of the problem under
consideration. The calculation of the MIs is then pursued by means of a variety of 
methods. Particularly important are the differential equations method \cite{Kotikov:1991kg,Kotikov:1991hm,Kotikov:1991pm,Remiddi:1997ny,Caffo:1998yd,Caffo:1998du,Argeri:2007up} and
the Mellin-Barnes techniques 
\cite{Smirnov:2004,Friot:2005cu,Usyukina:1975yg,Smirnov:1999gc,Tausk:1999vh,%
Smirnov:1999wz,Smirnov:2001cm,Heinrich:2004iq,Czakon:2005rk,Gluza:2007rt}.
Both methods proved to be very useful in the evaluation of virtual corrections to Bhabha 
scattering because they are especially effective in problems with a small number of different 
kinematic parameters. They both allow one to obtain an analytic expression for the 
integrals, which must be written in terms of a  suitable functional basis. A basis
which was extensively employed in the calculation of multi-loop Feynman diagrams of the 
type discussed here is represented by the Harmonic Polylogarithms \cite{Goncharov:1998,Broadhurst:1998rz,Remiddi:1999ew,Gehrmann:2001pz,Gehrmann:2001jv,Maitre:2005uu,Maitre:2007kp,Vollinga:2004sn,Weinzierl:2007cx} and their 
generalisations.
%
Another fundamental achievement which enabled one to complete the calculation of  the QED
two-loop corrections was an improved  understanding of the IR structure of QED. In
particular, the relation between the collinear logarithms in which the electron mass
$m_e$ plays  the role of a natural cut-off and the corresponding  poles  in the
dimensionally regularised  massless theory was extensively investigated in
\cite{Penin:2005kf,Penin:2005eh,Mitov:2006xs,Becher:2007cu}.

The first complete diagrammatic calculation of the two-loop QED virtual corrections  to
Bhabha scattering can be found in \cite{Bern:2000ie}. However, this result was obtained  in
the fully massless approximation ($m_e = 0$) by employing dimensional regularisation
(DR) to regulate  both soft and collinear divergences.
%
Today, the complete set of two-loop corrections to Bhabha scattering in pure QED have
been evaluated using $m_e$ as a collinear regulator, as required in order to include
these fixed-order calculations in available Monte Carlo event generators.
%
The Feynman diagrams involved in the calculation can be divided in three gauge-independent 
sets: {\em i)} diagrams without fermion loops (``photonic'' diagrams), {\em
ii)} diagrams involving a closed electron loop, and {\em iii)} diagrams involving a
closed loop of hadrons or a fermion heavier than the electron. Some of the diagrams
belonging to the aforementioned sets are shown in Figs.~\ref{fig0}--\ref{fig3}. 
These three sets are discussed in more detail below. \\

\noindent {\em{Photonic corrections \label{Phot}}} \\
\vspace*{-2mm}

\noindent A large part of the NNLO photonic corrections can be evaluated in a closed
analytic form, retaining  the full dependence on $m_e$ \cite{Bonciani:2005im}, by using the Laporta
algorithm for the reduction of the Feynman diagrams to a combination of MIs, and then
the differential equations method for their analytic evaluation.
%
With this technique it is possible to calculate, for instance, the NNLO corrections to
the form factors \cite{Bonciani:2003te,Bonciani:2003hc,Bonciani:2003ai,Czakon:2004wm}.
%
However, a calculation of the two-loop photonic boxes retaining the full dependence on
$m_e$ seems to be beyond the reach of this method. This is due to the fact that the number of
MIs belonging to the same topology is, in some cases, large.
Therefore, one must solve analytically large systems of first-order
ordinary linear differential equations; this is not possible in general.
Alternatively, in order to calculate the different MIs involved, one could use the
Mellin-Barnes techniques, as shown in \cite{Smirnov:2001cm,Heinrich:2004iq,Czakon:2004wm,Czakon:2006pa,Broadhurst:1993mw,Davydychev:2003mv}, or a combination 
of both methods. The calculation is very complicated and a full result is not 
available yet.\footnote{For the planar double box diagrams, all the MIs are known 
\cite{Czakon:2006pa} for small $m_e$, while the MIs for the non-planar double box diagrams 
are not completed.} However, the full dependence on $m_e$ is not phenomenologically relevant.
%
In fact, the physical problem exhibits a well defined mass hierarchy. The mass of the
electron is always very small compared to the other kinematic invariants and can be
safely neglected everywhere, with the exception of the terms in which it acts as a
collinear regulator. The ratio of the photonic NNLO corrections to the Born cross
section is given by 
%\footnote{Infrared logarithms have been omitted in 
%Eq.~(\ref{exp-photonic})  and following ones, because they disappear when summing up the
%contribution of hard photons and inclusive photon energy conditions are considered.}
%
\be
\frac{{\rm d} \sigma^{(2,\mbox{\tiny{PH}})}}{{\rm d} \sigma^{(\mbox{\tiny{Born}})}}  =
\left( \frac{\alpha}{\pi} \right)^2 \sum_{i=0}^{2} \delta^{(\mbox{\tiny{PH}},i)} 
\left(L_e \right)^i +
O \left(\frac{m_e^2}{s}, \frac{m_e^2}{t}\right) \, ,
\label{exp-photonic}
\ee
%
where $L_e=\ln{(s/m_e^2)}$ and the coefficients $\delta^{(\mbox{\tiny{PH}},i)}$
contain infrared logarithms and are functions of the scattering angle $\theta$. The approximation given by
Eq.~(\ref{exp-photonic}) is sufficient for a phenomenological description of the
process.\footnote{It can be shown that the terms suppressed by a positive power of
$m_e^2/s$ do not play any phenomenological role already at very low c.m. energies,
$\sqrt{s} \sim 10$ MeV. Moreover, the terms $m_e^2/t$ (or $m_e^2/u$) become important
in the extremely forward (backward) region, unreachable for the experimental setup.}
%
% Actis 1
%%%
\begin{figure}
\begin{center}
\resizebox{0.45\textwidth}{!}{%
\includegraphics{Actis_fig1}
}
%\vspace*{4mm}
\caption{\small{ Some of the diagrams belonging to the class of the ``photonic'' NNLO
corrections to the Bhabha scattering differential cross section. The additional
	photons in the final state are soft. \label{fig0}}}
\end{center}	
\end{figure}
%
The coefficients of the double and single collinear logarithm in
Eq.~(\ref{exp-photonic}), $\delta^{(\mbox{\tiny{PH}},2)}$ and
$\delta^{(\mbox{\tiny{PH}},1)}$,  were obtained in \cite{Arbuzov:1998du,Glover:2001ev}. However, the
precision required for luminosity  measurements at $e^+ e^-$ colliders demands the
calculation of the  non-logarithmic coefficient, $\delta^{(\mbox{\tiny{PH}},0)}$.  The
latter  was  obtained in \cite{Penin:2005kf,Penin:2005eh} by reconstructing the  differential cross section
in the $s \gg m_e^2 \neq 0$ limit from the dimensionally regularised massless 
approximation \cite{Bern:2000ie}. The main idea of the method developed in \cite{Penin:2005kf,Penin:2005eh} is
outlined below: As far as the leading term in the small electron mass expansion is
considered, the difference between the massive and the dimensionally regularised
massless Bhabha scattering can be viewed as a difference between two regularisation
schemes for the infrared divergences. With the known massless two-loop result at hand,
the calculation of the massive one is reduced to constructing the {\it infrared
matching} term which relates the two above  mentioned regularisation schemes. To
perform the matching an auxiliary amplitude is constructed, which has the same
structure of the infrared singularities but is sufficiently simple to be evaluated at
least at the leading order in the small mass expansion. The particular form of the auxiliary
amplitude is dictated by the general theory of infrared singularities in QED and
involves the exponent of the one-loop correction as well as the two-loop corrections to
the logarithm of the electron form factor. The difference between the full and
the auxiliary amplitudes is infrared finite.  It can be evaluated by using dimensional
regularisation for each amplitude and then taking the limit of four space-time
dimensions. The infrared divergences, which induce the asymptotic dependence of the
virtual corrections on the electron and photon masses, are absorbed into the auxiliary
amplitude while the technically most nontrivial calculation of the full amplitude  is
performed in the massless approximation. The matching of the massive and massless
results is then necessary only for the auxiliary amplitude and is straightforward.
Thus the two-loop massless result for the scattering amplitude along  with the two-loop
massive electron form factor \cite{Burgers:1985qg} are sufficient to obtain the two-loop photonic
correction to the differential cross section in the small electron mass limit.

A method based on a similar principle was subsequently developed in
\cite{Mitov:2006xs,Becher:2007cu}; the authors of \cite{Becher:2007cu} confirmed  the result of 
\cite{Penin:2005kf,Penin:2005eh} for the NNLO photonic corrections to the Bhabha  scattering differential cross section.  \\

%
% Actis 2
%%%
\begin{figure}%[b]
\begin{center}
%\begin{figure}
\resizebox{0.45\textwidth}{!}{%
\includegraphics{Actis_fig2.ps}
}
%\vspace*{4mm}
\caption{\small{ Some of the diagrams belonging to the class of the ``electron loop'' 
NNLO corrections. The additional photons or electron-positron pair in the final state 
are soft. \label{fig1}}}
\end{center}
\end{figure}
%
%\be
%\frac{d \sigma^{(2,\mbox{\tiny{EL}})}}{d \sigma^{(\mbox{\tiny{Born}})}}  =
%\frac{\alpha^2}{\pi^2} \sum_{i=0}^{3} \delta^{(\mbox{\tiny{EL},i})} L_e^i +
%\mathcal{O}\left(\frac{m_e^2}{s}, \frac{m_e^2}{t}\right) \, .
%\label{exp-electron}
%\ee
%

\noindent {\em{Electron loop corrections \label{ElLoop}}} \\
\vspace*{-2mm}

\noindent The NNLO electron loop corrections arise from the interference of two-loop 
Feynman diagrams with the tree-level amplitude as well as from the interference of
one-loop diagrams, as long as one of the diagrams contributing to each term involves a
closed electron loop. This set of corrections presents a single two-loop box
topology and is therefore technically less challenging to evaluate  with respect to
the photonic correction set.  The calculation of the electron loop corrections was
completed a few years ago \cite{Bonciani:2003cj0,Bonciani:2004gi,Bonciani:2004qt,Actis:2007gi}; the final result retains the full
dependence of the differential cross section on the electron mass $m_e$.  The MIs
involved in the calculation were identified by means of the Laporta algorithm and
evaluated with the differential equation method. As expected, after UV renormalisation
%%%TT
the differential cross section contained only  residual IR  poles which were removed by
adding the contribution of the  soft photon emission diagrams.  The resulting NNLO
differential cross section could be conveniently written in terms of  1- and
2-dimensional Harmonic Polylogarithms (HPLs) of maximum weight three. Expanding the
cross section in the limit $s, |t| \gg m_e^2$, the ratio of the NNLO corrections to the
Born cross section can be  written as in Eq.~(\ref{exp-photonic}):
\be
\frac{{\rm d} \sigma^{(2,\mbox{\tiny{EL}})}}{{\rm d} \sigma^{(\mbox{\tiny{Born}})}}  =
\left( \frac{\alpha}{\pi} \right)^2 \sum_{i=0}^{3} \delta^{(\mbox{\tiny{EL},i})} 
\left(L_e\right)^i +
O \left(\frac{m_e^2}{s}, \frac{m_e^2}{t}\right) \, .
\label{exp-electron}
\ee
Note that the series now contains a cubic collinear logarithm. This logarithm
appears, with an opposite sign, in the corrections due to the production of an
electron-positron pair (the soft-pair production was considered in
\cite{Arbuzov:1995vj}). When the two contributions are considered together in the full
NNLO, the cubic collinear logarithms cancel. Therefore, the physical cross section includes 
at most a double logarithm, as in Eq.~(\ref{exp-photonic}).

The explicit expression of all the coefficients $\delta^{(\mbox{\tiny{EL},i})}$,
obtained by expanding the results of \cite{Bonciani:2003cj0,Bonciani:2004gi,Bonciani:2004qt}, was confirmed by two different
groups \cite{Becher:2007cu,Actis:2007gi}. In \cite{Becher:2007cu} the small electron mass expansion was
performed within the soft-collinear effective theory (SCET)  framework, while  the
analysis in \cite{Actis:2007gi} employed  the asymptotic expansion of the  MIs.  \\


\noindent {\em{Heavy-flavor and hadronic corrections \label{HFLoop}}} \\
\vspace*{-2mm}

\noindent Finally, we consider the corrections originating from two-loop Feynman
diagrams involving a heavy flavour fermion loop.\footnote{Here by ``heavy flavour'' we
mean a muon or a $\tau$-lepton, as well as a heavy quark, like the top, the $b$- or the
$c$-quark, depending on the c.m. energy range that we are considering.} Since this set
of corrections  involves one more mass scale with respect to the corrections analysed
in the previous sections, a direct diagrammatic calculation is in principle a more
challenging task.
%
Recently, in \cite{Becher:2007cu} the authors applied their technique based on SCET to Bhabha
scattering and obtained the heavy flavour NNLO corrections in the limit in which $s,
|t|, |u| \gg m_f^2 \gg m_e^2$, where $m_f^2$ is the mass of the heavy fermion running
in the loop. Their result was very soon confirmed in \cite{Actis:2007gi} by means of a method
based on the asymptotic expansion of Mellin-Barnes representations of the MIs involved in the calculation.
%
However, the results obtained in the approximation $s, |t|, |u| \gg m_f^2 \gg m_e^2$
cannot be applied to the case in which $\sqrt{s} < m_f$ (as in the case of a tau
loop at $\sqrt{s} \sim 1$ GeV), and they apply only to a relatively narrow angular
region perpendicular to the beam direction when $\sqrt{s}$ is not very much larger than
$m_f$ (as in the case of top-quark loops at the ILC). It was therefore necessary to
calculate the heavy flavour corrections to Bhabha scattering assuming only that the
electron mass is much smaller than the other scales in the process, but retaining the
full dependence on  the heavy mass, $s, |t|, |u|, m_f^2 \gg m_e^2$.

The calculation was carried out in two different ways: in \cite{Bonciani:2007eh,Bonciani:2008ep} it was
done analytically, while in \cite{Actis:2007,Actis:2008br} it was done numerically with dispersion relations.

%%%%%%%%%%%%%%%%%%%

The technical problem of the diagrammatic calculation of Feynman integrals with four
scales can be simplified by considering carefully, once more, the structure of the
collinear singularities of the  heavy-flavour  corrections. The ratio of the NNLO heavy
flavour corrections to the Born cross section is given by
%
\be \label{exp-heavyflavor}
\frac{{\rm d} \sigma^{(2,\mbox{\tiny{HF}})}}{{\rm d} \sigma^{(\mbox{\tiny{Born}})}}  =
\left( \frac{\alpha}{\pi} \right)^2 \sum_{i=0}^{1} \delta^{(\mbox{\tiny{HF},i})} 
\left(L_e\right)^i +
O \left(\frac{m_e^2}{s}, \frac{m_e^2}{t}\right) \, ,
\ee
where now the coefficients $\delta^{(i)}$ are functions of the scattering angle
$\theta$ and, in general, of the mass of the heavy fermions involved in the virtual
corrections.
%
It is possible to prove that, in a physical gauge, all the collinear singularities
factorise and can be absorbed in the external field renormalisation \cite{Frenkel:1976bj}. This
observation has two consequences in the case at hand. The first one is that box
diagrams are free of collinear divergences in a physical gauge; since the sum of all
boxes forms a gauge independent block, it can be concluded that the sum of all box
diagrams is free of collinear divergences in any gauge. The second consequence is that
the single collinear logarithm in Eq.~(\ref{exp-heavyflavor}) arises from vertex
corrections only. Moreover, if one chooses on-shell UV renormalisation conditions, the
irreducible two-loop vertex graphs are free of collinear singularities.
%
% Actis 3
%%
%%%\begin{figure}
\begin{figure}
\begin{center}
\resizebox{0.45\textwidth}{!}{%
\includegraphics{Actis_fig3.ps}
}
\vspace*{4mm}
\caption{\small{ Some of the diagrams belonging to the class of the ``heavy fermion'' 
NNLO corrections. The additional photons in the final state are soft.
\label{fig2}}}
\end{center}
\end{figure}
%
Therefore, among all the two-loop diagrams contributing to the NNLO heavy flavour
corrections to Bhabha scattering, only the reducible vertex corrections are
logarithmically divergent in the $m_e \to 0$ limit.\footnote{Additional collinear
logarithms arise also from the interference of one-loop diagrams in which at least one
vertex is present.} The latter are easily evaluated even if they depend on two
different masses.
%
By exploiting these two facts, one can obtain the NNLO heavy-flavour corrections to the
Bhabha scattering differential cross section assuming only that $s, |t|, |u|, m_f^2 \gg
m_e^2$. In particular, one can set $m_e =0$ from the 
beginning in all the two-loop diagrams
with the exception of the reducible ones. This procedure allows one  to effectively
eliminate one mass scale from the two-loop boxes, so that these graphs can be  evaluated
with the techniques already employed in the diagrammatic calculation of the electron
loop corrections.\footnote{{The necessary MIs can be found in \cite{Bonciani:2008ep,Fleischer:1998nb,Aglietti:2003yc,Aglietti:2004tq}.}}
%
In the case in which the heavy flavour fermion is a quark, it is straightforward to
modify the calculation of the two-loop self-energy diagrams to obtain the mixed QED-QCD
corrections to Bhabha scattering \cite{Bonciani:2008ep}.

%%%%%%%%%%%%

An alternative approach to the calculation of the heavy flavour corrections to
Bhabha scattering is based on dispersion relations. 
%\GR{%
This method also applies to hadronic corrections.
The  hadronic and heavy fermion corrections to the Bhabha-scattering cross section
can be obtained by appropriately 
inserting the renormalised irreducible photon
vacuum-polarisation function  $\Pi$  in the photon propagator:
\begin{equation}
\label{1stReplace}
\frac{g_{\mu\nu}}{q^2+i\, \delta} \, \to \,
\frac{g_{\mu\alpha}}{q^2+i\, \delta} \,
\left( q^2\, g^{\alpha\beta} - q^\alpha\, q^\beta \right) \,\Pi(q^2)\,
\frac{g_{\beta\nu}}{q^2+i\, \delta}.
\end{equation}
The vacuum polarisation $\Pi$ 
can be represented by a once-subtracted dispersion integral \cite{Cabibbo:1961sz},
\begin{equation}
\label{DispInt}
\Pi(q^2) =
- \frac{q^2}{\pi} \, 
  \int_{4 M^2}^{\infty} \, {\rm d} z \, 
  \frac{\textrm{Im} \, \Pi(z)}{z} \, 
  \frac{1}{q^2-z+i\, \delta} .
\end{equation}
The contributions to $\Pi$ may then be determined from a (properly normalised) production cross 
section by the optical theorem~\cite{Cutkosky:1960sp},
\begin{eqnarray}
\label{Rhad0}
\textrm{Im} \, \Pi_{\rm had}(z)&=& 
- \frac{\alpha}{3} \, R(z).
\end{eqnarray}
In this way, the hadronic vacuum polarisation may be obtained 
from the experimental data for  $R$: 
\begin{eqnarray}
\label{Rhad}
R(z) &=& 
\frac{\sigma^0_{\rm had}(z)}
     {(4 \pi \alpha^2)\slash (3z)} \, ,
\end{eqnarray}
where $\sigma^0_{\rm had}(z)\equiv \sigma(\{e^+e^-\to\gamma^\star\to \textrm{hadrons}\};z)$.
In the low-energy region the inclusive experimental data may be used \cite{Davier:2002dy,rintpl:2008AA}.
Around a narrow hadronic resonance with mass $M_{\rm res}$ and width $\Gamma^{e^+ e^-}_{\rm res}$ one may use the relation
\begin{equation}
R_{\rm res}(z)= \frac{9 \pi}{\alpha^2} M_{\rm res} \Gamma^{e^+ e^-}_{\rm res}
\delta(z-M^2_{\rm res}),
\end{equation}
and in the remaining regions the perturbative QCD prediction \cite{Harlander:2002ur}.
Contributions to $\Pi$ arising from leptons and heavy  quarks with mass $m_f$, charge $Q_f$ and colour $C_f$ can be 
computed directly in perturbation theory. In the lowest order it reads
\begin{eqnarray}
R_f(z;m_f)&=&
 Q_f^2\, C_f \,\left(1+2\, \frac{m_f^2}{z}\right)\,
        \sqrt{1-4\,\frac{m_f^2}{z}}.
\end{eqnarray}
As a result of the above formulas, the massless photon propagator gets replaced by a massive 
propagator, whose effective mass $z$ is subsequently integrated over:
\begin{equation}
\label{PropReplace}
\frac{g_{\mu\nu}}{q^2+i \delta}
\to
\frac{\alpha}{3 \pi}
\int_{4 M^2}^{\infty} 
\frac{{\rm d}z~R_{\rm tot}(z)}{z(q^2-z+i \delta)}
\left(
g_{\mu\nu} - \frac{q_\mu q_\nu}{q^2+i \delta}
\right)\, ,
\end{equation}
where $R_{\rm tot}(z)$ contains hadronic and leptonic contributions.

For self-energy corrections to Bhabha scattering at one-loop order,  the dispersion relation approach  was first employed 
in \cite{Berends:1976zn}.
Two-loop applications of this technique, prior to Bhabha scattering, are the evaluation of the had\-ronic 
vertex correction \cite{Kniehl:1988id} and of two-loop had\-ronic corrections to the lifetime of the muon \cite{vanRitbergen:1998hn}. The approach was also applied to the evaluation of the two-loop form factors in QED in \cite{Barbieri:1972as,Barbieri:1972hn,Mastrolia:2003yz}.


The fermionic and hadronic corrections to Bhabha scattering at one-loop accuracy come only from  the self-energy diagram; see for details Section~\ref{sec:4}.
%At two loop accuracy, there are a lot of \emph{reducible diagrams} composed of the one-loop insertions of  Fig.~\ref{OneLoop}(c) and of all the one-loop topologies Fig.~\ref{OneLoop}(a-c).
At two-loop level  there are  
reducible and  irreducible  self-energy contributions, vertices and boxes.
The reducible corrections are easily treat\-ed. 
For the evaluation of the irreducible two-loop diagrams, it is advantageous that they are one-loop diagrams with self-energy insertions because the application of the dispersion technique as described here is possible.

The kernel function for the irreducible two-loop vertex was derived in \cite{Kniehl:1988id} and verified e.g. in \cite{Actis:2008br}. The three kernel functions for the two-loop box functions were first obtained in \cite{Actis:2007pn,Actis:2007,Actis:2008br} and verified in \cite{Kuhn:2008zs}. 
A complete collection of all the relevant formulae may be found  in \cite{Actis:2008br}, and the corresponding Fortran code bhbhnnlohf is publicly available 
at the web page \cite{ACGR:bhbhnnlohf} \\ 
%\href{http://www-zeuthen.desy.de/theory/research/bhabha/bhabha.html}.
{\tt www-zeuthen.desy.de/theory/research/bhabha/} .

In \cite{Actis:2008br}, the dependence of the various heavy fermion NNLO corrections on $\ln(s/m_f^2)$ for $s,|t|,|u| \gg m_f^2$ was studied. The irreducible vertex behaves (before a combination with real pair emission terms) like $\ln^3(s/m_f^2)$ \cite{Kniehl:1988id}, while the sum of the various infrared divergent diagrams as a whole behaves like $\ln(s/m_f^2)\ln(s/m_e^2)$.
This is in accordance with 
Eq.~(\ref{exp-heavyflavor}), but the limit plays no effective role at the energies studied here.

As a result of the efforts of recent years we now have at least two completely independent calculations for all the non-photonic virtual two-loop contributions. 
The net result, as a ratio of the NNLO corrections to the Born cross section in per mill, is shown in 
Fig.~\ref{1gevnoSE} for KLOE and in Fig.~\ref{10gevnoSE} 
for BaBar/Belle.\footnote{The pure self-energy corrections deserve a special 
discussion and are thus omitted in the plots.}
While the non-photonic corrections stay at one per mill or less for  KLOE, they reach 
a few per mill at the BaBar/Belle energy range. The NNLO 
photonic corrections are the dominant contributions and amount 
to some per mill, both at $\mathrm{\phi}$ and $B$ factories. However, as already emphasised, 
the bulk of both photonic and non-photonic 
corrections is incorporated into the generators used by the experimental collaborations. Hence, 
the consistent comparison between the results of NNLO calculations and the MC predictions 
at the same perturbative level enables one to assess the theoretical accuracy of the luminosity tools, 
as will be discussed quantitatively in Section \ref{TH}.


\begin{figure}
\begin{center}
\resizebox{0.45\textwidth}{!}{%
\includegraphics{Actis_fig4.ps}
}
%\vspace*{4mm}
\caption{\small{ Some of the diagrams belonging to the class of the ``hadronic'' 
corrections. The additional photons in the final state are soft.
\label{fig3}}}
\end{center}
\end{figure}



\begin{figure}[t]
\begin{center}
\vspace*{7mm}
%\hspace*{-2mm}
\begin{center}
\scalebox{0.33}{\includegraphics{1gevnoSE-try.eps}}
\end{center}
\caption[]
{Two-loop photonic and non-photonic corrections to Bhabha scattering at $\sqrt{s}=1.02$ GeV, normalised
to the QED tree-level cross section, as a function of the electron polar angle; no cuts; the parameterisations of $R_{\rm had}$ from \cite{Burkhardt:1981jk} and \cite{Davier:2002dy,rintpl:2008AA,Harlander:2002ur} are very close to each other. }
\label{1gevnoSE}
\end{center}
\end{figure}

%--\scalebox{h-scale}[v-scale]{argument}
%\begin{figure}[t]
%\begin{center}
%\vspace*{7mm}
%\hspace*{4mm}
%\begin{center}
%\begin{picture}(0,0)
%\scalebox{0.34}{\includegraphics{1gevnoSE}}
%\end{picture}%
%\end{center}
%\caption[]{Two-loop corrections to Bhabha scattering at $\sqrt{s}=1.02$ GeV, normalised
%to the QED tree-level cross section. No cuts. The parameterisation of $R_{had}$  due to \cite{Burkhardt:1981jk} and \cite{Davier:2002dy,rintpl:2008AA,Harlander:2002ur} are very close to each other.}
%\label{1gevnoSE}
%\end{center}
%\end{figure}

\begin{figure}[t]
\begin{center}
\vspace*{7mm}
%\hspace*{-2mm}
\begin{center}
\scalebox{0.33}{\includegraphics{10gevnoSE-try.eps}}
%\includegraphics[width=0.2\textwidth]{10gevnoSE-try.eps}
\end{center}
\caption[]
{Two-loop photonic and non-photonic corrections to Bhabha scattering 
at $\sqrt{s}=10.56$ GeV, normalised
to the QED tree-level cross section, 
as a function of the electron polar angle; no cuts; the parameterisations of $R_{\rm had}$ is from \cite{Burkhardt:1981jk}.}
\label{10gevnoSE}
\end{center}
\end{figure}


%%%
% Actis 4
%%


\subsubsection{Fixed-order calculation of the hard photon emission at one loop
\label{Penta}}

The one-loop matrix element for the process $e^{+}e^{-} \to e^{+}e^{-} \gamma$ is one
of the contributions to the complete set of NNLO corrections to Bhabha scattering. Its
evaluation requires the nontrivial computation of one-loop tensor integrals associated
with pentagon diagrams.

According to the standard Passarino-Veltman (PV) approach \cite{Passarino:1978jh},
one-loop tensor integrals can be expressed in terms of 
MIs with trivial numerators that are independent of the loop variable,
each multiplied by a Lorentz structure depending only on combinations of the external
momenta and the metric tensor.
%
The achievement of the complete PV-reduction amounts to solving a nontrivial system of
equations. 
%\BL{TR: The next statement is not clear to me. I try to modify, please check if it is ok with you.}
%\RE{\sout{Due to its size, it is reasonable replacing the analytic techniques by
%numerical tools.}}
Due to its size, it is reasonable to replace the analytic techniques by
numerical tools.
It is difficult to implement the PV-reduction numerically, since it gives rise to Gram
determinants. The latter naturally arise in the procedure of inverting a system and they 
can vanish at special phase space points. This fact requires a proper modification of the 
reduction algorithm \cite{Bern:1996je,Dixon:1996wi,Binoth:1999sp,Denner:2002ii,Binoth:2005ff,Denner:2005nn,Denner:2006fy}. A viable solution for the complete algebraic reduction of 
tensor-pentagon (and tensor-hexagon) integrals was formulated in \cite{Davydychev:1991va,Tarasov:1996br,Fleischer:1999hq}, by exploiting the algebra of signed minors \cite{Melrose:1965kb}.
In this approach the cancellation of powers of inverse Gram determinants was performed recently
in \cite{Diakonidis:2008dt,Diakonidis:2008ij}.

%\BL{
%Generally I find the whole text a bit long for having no numerical results ... but why not ...
%}


Alternatively, the computation of the one-loop five-point amplitude $e^{+}e^{-} \to e^{+}e^{-} \gamma$
can be performed by using generalised-unitarity cutting rules (see
\cite{Bern:2007dw} for a detailed compilation of references). In the following we
propose two ways to achieve the result, via an analytical and via a 
semi-numerical method. The application of generalised cutting rules as an on-shell
method of calculation is based on two fundamental properties of scattering amplitudes:
{\it i)} analyticity, according to which any amplitude is determined by its
singularity structure \cite{Landau:1959fi,Mandelstam:1958xc,Mandelstam:1959bc,Cutkosky:1960sp,Eden:1966}; and {\it ii)} unitarity, according to which
the residues at the singularities are determined by products of simpler amplitudes. 
Turning these properties into a tool for computing scattering amplitudes is possible
because of the underlying representation of the amplitude in terms of Feynman integrals
and their PV-reduction, which grants the existence of a representation of any one-loop
amplitudes as linear combination of MIs, each multiplied by a rational coefficient. In
the case of $e^{+}e^{-} \to e^{+}e^{-} \gamma$, pentagon-integrals 
%\RE{\sout{should be ultimately} 
may be
expressed, through PV-reduction, by a linear combination of 17 MIs (including 3 boxes,
8 triangles, 5 bubbles and 1 tadpole).
 Since the required MIs are analytically known
\cite{'tHooft:1978xw,Bern:1992em,Bern:1993kr,Tarasov:1996br,Binoth:1999sp,Duplancic:2003tv,Ellis:2007qk},
the determination of their coefficients is needed for reconstructing the amplitude
as a whole.
% To this aim, one may use the Mathematica program hexagon \cite{Diakonidis:2008dt,Diakonidis:2008ij}. Also 
Matching the generalised cuts of the amplitude with the cuts
of the MIs provides an efficient way to extract their (rational) coefficients from the amplitude itself. 
%
In general the fulfilment of multiple-cut conditions requires loop momenta with
complex components. The effect of the cut conditions is to freeze some or all of its
components, depending on the number of the cuts. With the
quadruple-cut \cite{Britto:2004nc} the loop momentum is completely frozen, yielding
the algebraic determination of the coefficients of $n$-point functions with $n\ge4$. In
cases where fewer than four denominators are cut, like triple-cut
\cite{MastroliaTriple,FordeTriBub,BjerrumBohr:2007vu}, double-cut
\cite{Britto:2005ha,Britto:2006sj,Anastasiou:2006jv,Anastasiou:2006gt,Britto:2006fc,FordeTriBub} and single-cut
\cite{NigelGlover:2008ur}, the loop momentum is not frozen: the free components are
left over as phase-space integration variables. 

For each multiple-cut, the
evaluation of the phase-space integral would generate, in general, logarithms and a
non-logarithmic term. The coefficient of a given $n$-point MI finally appears in the
non-logarithmic term of the corresponding $n$-particle cut, where all the internal lines
are on-shell (while the logarithms correspond to the cuts of higher-point MIs which
share that same cut). Therefore all the coefficients of MIs can be determined in a
{\it top-down} algorithm, starting from the quadruple-cuts for the extraction of the
four-point coefficients, and following with the triple-, double- and single-cuts for the
coefficients of three-, two- and one-point, respectively. 
%
The coefficient of an $n$-point MI ($n \ge 2$) can also be obtained by specialising the generating formulas given in \cite{Britto:2008vq} for general one-loop amplitudes to the case at hands.


Instead of the analytic evaluation of the multiple-cut phase-space integrals, it is
worth considering the feasibility of computing the process $e^{+}e^{-} \to e^{+}e^{-}
\gamma$ with a semi-numerical technique by now known as OPP-reduction
\cite{Ossola:2006us,Ossola:2007bb}, based on the decomposition of the numerator of any
one-loop integrand in terms of its denominators \cite{delAguila:2004nf,Pittau:2004bc,Pittau:1996ez,Pittau:1997mv}. Within this approach
the coefficients of the MIs can be found simply by solving a system of numerical
equations, avoiding any explicit integration. The OPP-reduction algorithm
exploits the polynomial structures of the integrand when evaluated at values of the
loop-mo\-men\-tum fulfilling multiple cut-conditions: {\it i)} for each $n$-point MI one
considers the $n$-particle cut obtained by setting all the propagating lines on-shell;
{\it ii)} such a cut is associated with a polynomial in terms of the free components of
the loop-momentum, which corresponds to the numerator of the integrand evaluated at the
solution of the on-shell conditions; {\it iii)} the constant-term of that polynomial
{\it is} the coefficient of the MI. \\ Hence the difficult task of evaluating one-loop
Feynman integrals is reduced to the much simpler problem of polynomial fitting,
recently optimised by using a projection technique based on the Discrete Fourier
Transform \cite{Mastrolia:2008jb}.

In general the result of a dimensional-regulated amplitude in the 4-dimensional limit,
with $D$ ($=4-2\epsilon$) the regulating parameter, is expected to contain
(poly)logarithms, often referred to as the cut-constructible term, and a pure 
rational term. In a later paper \cite{Ossola:2008xq}, which completed the OPP-method, the rising of the rational term was attributed to two potential
sources (of UV-divergent integrals): one, defined as $R_1$, due to the $D$-dimensional
completion of the 4-dimensional contribution of the numerator; a second one, called
$R_2$, due to the $(-2\epsilon)$-dimensional algebra of Dirac-matrices. Therefore
in the OPP-approach the calculation of the one-loop amplitude $e^{+}e^{-} \to
e^{+}e^{-} \gamma$ can proceed through two computational stages:
\begin{enumerate}
\item
the coefficients of the MIs that are responsible both for the cut-constructible and
for the $R_1$-rational terms can be determined by applying the OPP-reduction
discussed above \cite{Ossola:2006us,Ossola:2007bb,Mastrolia:2008jb};
\item 
the $R_2$-rational term can be computed by using additional tree-level-like
diagrammatic rules, very much resembling the computation of the counter terms needed for
the renormalisation of UV-divergences \cite{Ossola:2008xq}.
\end{enumerate}

The numerical influence of the radiative loop diagrams, including the pentagon diagrams, is expected 
not to be  
particularly large. However, the calculation of such corrections would greatly help to 
assess the physical precision of existing luminosity programs.\footnote{As already remarked, the exact calculation of one-loop
corrections to hard photon emission in Bhabha scattering became available~\cite{Actis:2009uq}
 during the completion of the report, exactly according to the
methods described in the present Section.}

 %and not being logarithmically enhanced (??). Nevertheless, we expect some predictions for them in the near future.}


%\subsubsection{Four Particles in the Final State \label{4pfinal}}
\subsubsection{Pair corrections \label{4pfinal}}
As was mentioned  in the paragraph on virtual heavy fla\-vour and hadronic corrections of Section~\ref{virtual}, these virtual corrections have to be combined with real corrections in order to get physically sensible results.
The virtual NNLO electron, muon, tau and pion corrections have to be combined with the emission of real electron, muon, tau and pion pairs, respectively.
The real pair production cross sections are finite, but cut dependent. 
We consider here the pion pair production as it is the dominant part of the hadronic corrections and 
can serve as an estimate of the role of the whole set of hadronic corrections. The description of all relevant hadronic contributions is a much more involved task and will not be covered in this review. 
As was first explicitly shown for Bhabha scattering in \cite{Arbuzov:1995vj} for electron pairs, and also discussed in \cite{Actis:2008br}, there appear exact cancellations of terms of the order $\ln^3(s/m_e^2)$ or $\ln^3(s/m_f^2)$, so that the leading terms are at most of order $\ln^2(s/m_e^2), \ln^2(s/m_f^2)$.

\begin{table}[t]
%\setlength{\tabcolsep}{0.3pc}
\caption{\label{tab-pairs}%
The NNLO lepton and pion pair corrections to the Bhabha scattering Born cross section $\sigma_B$:
 virtual corrections $\sigma_v$ , soft and hard real photon emissions $\sigma_s, \sigma_h$, and pair emission contributions $\sigma_{pairs}$. The total pair correction cross sections are obtained from the sum $\sigma_{s+v+h} + 
\sigma_{pairs}$. All cross sections, %with
%an angular acceptance cut $|\cos\theta|<0.7$
according to the cuts given in the text, are given in nanobarns.}
% numbers fro Gosia: status 14-05-2009 07:38 PDT 
\begin{center}
\begin{tabular}{cllllll}
\hline %-------------------------------------------------
\multicolumn{6}{c}{\bf{Electron pair corrections}}
\\
       &$\sigma_{B}$ & $\sigma_{h}$ & $\sigma_{v+s}$ & $\sigma_{v+s+h}$ & $\sigma_{pairs}$ 
\\
\hline
% KLOE & \RE{529.055} & 8.632  &  -11.564 & -2.932 & 0.271201 +- 0.6%
  KLOE & 529.469 & 9.502  &  -11.567 & -2.065 & 0.271 
\\
% BABAR & \RE{6.74835} & 0.2963 & -0.2714 & 0.0249 & 0.0166305 +- 0.7% &\\
BaBar & 6.744 & 0.246 & -0.271 & -0.025 & 0.017 
\\ \hline
\multicolumn{6}{c}{
\bf{Muon pair corrections}} %-------------------------------------------------
\\
& $\sigma_{B}$ & $\sigma_{h}$ & $\sigma_{v+s}$ & $\sigma_{v+s+h}$ & $\sigma_{pairs}$ 
\\
\hline
% KLOE & \RE{529.055} & 1.3567  & -1.5983  & -0.2416 &  \RE{--}
KLOE & 529.469 & 1.494  & -1.736  & -0.241 &  --
\\
% BABAR & \RE{6.74835}  & 0.11155 & -0.11370 & -0.00215 & 0.00\RE{0.0005318} \\
BaBar & 6.744  & 0.091 & -0.095 & -0.004 & 0.0005 \\
\hline
\multicolumn{6}{c}{
\bf{Tau pair corrections}} %-------------------------------------------------
\\
& $\sigma_{B}$ & $\sigma_{h}$ & $\sigma_{v+s}$ & $\sigma_{v+s+h}$ & $\sigma_{pairs}$ 
\\
\hline
% KLOE&  \RE{529.055}  & 0.01830  & -0.02155  & -0.00325 & \RE{--}
KLOE&  529.469  & 0.020  & -0.023  & -0.003 & --
%$10^{-6}$ &   &   &  \\
\\
% BABAR&  \RE{6.74835} & 0.01929 & -0.01999 & -0.00070 & \RE{--}\\
BaBar &  6.744 & 0.016 & -0.017 & -0.0007 & $< 10^{-7}$ \\
\hline
\multicolumn{6}{c}{
\bf{Pion pair corrections}} %-------------------------------------------------
\\
& $\sigma_{B}$ & $\sigma_{h}$ & $\sigma_{v+s}$ & $\sigma_{v+s+h}$ & $\sigma_{pairs}$ 
\\
\hline
% KLOE & \RE{529.055}  & 1.0660 &  -1.2517 & -0.1857 & 0\RE{--}
KLOE & 529.469  & 1.174 &  -1.360 & -0.186 & --
\\
% BABAR & \RE{6.74835}  & 0.07497  & -0.07754  & -0.00257 & 0.00029\\
BaBar & 6.744  & 0.062  & -0.065  & -0.003 & 0.00003
\\
\hline
\end{tabular}
\end{center}
%\end{tabular*}
\end{table}


In Table~\ref{tab-pairs} we show NNLO 
lepton and pion pair contributions with typical kinematical cuts for the KLOE and BaBar experiments.
Besides contributions from unresolved pair emissions $\sigma_{pairs}$,  we also add unresolved real hard photon emission contributions $\sigma_{h}$.
The corrections  $\sigma_{pairs}$ from fermions have been calculated with the Fortran 
package HELAC-PHEGAS
 \cite{Kanaki:2000ey,Papadopoulos:2000tt,Papadopoulos:2005ky,Cafarella:2007pc}, 
the real pion corrections with EKHARA %\cite{Czyz:2007ue,Czyz:2008zz}},
\cite{Czyz:2005ab,Czyz:2006dm},  
the NNLO hard photonic corrections  $\sigma_{h}$ with a program \cite{chunp} based 
on the  generator  BHAGEN-1PH \cite{Caffo:1996mi}.
The latter depend, technically, on the soft photon cut-off $E_{\gamma}^{\min} = \omega$.
After adding up with  $\sigma_{v+s}$, the sum of the two $\sigma_{v+s+h}$ is independent of that; in fact here we use $\omega/E_{\rm beam} = 10^{-4}$.
In order to cover also pion pair corrections $\sigma_{v+s}$ is determined with an updated version of the Fortran package bhbhnnlohf \cite{Actis:2008br,ACGR:bhbhnnlohf}. 
The cuts 
%which single out Bhabha events have been 
applied in Table~\ref{tab-pairs} for the  KLOE experiment are
% \begin{itemize}
%for the  $\Phi$ factory/DA$\Phi$NE (Frascati):
\begin{itemize}
\item $\sqrt{s} = 1.02 ~\textnormal{GeV}$\,,
\item  $E_{\rm min} = 0.4 ~\textnormal{GeV}$\,, 
\item $55^{\circ} < \theta_{\pm}<125^{\circ}$\,,
\item  $\xi_{\rm max} = 9^{\circ}$\,,
\end{itemize}
and for the BaBar experiment
\begin{itemize}
\item $\sqrt{s} = 10.56 ~\textnormal{GeV}$\,,
\item  $|\cos(\theta_{\pm})|<0.7 
~~\textnormal{and} 
\\
|\cos(\theta_{+})|<0.65 
~~\textnormal{or} ~~|\cos(\theta_{-})|<0.65$\,,  
\item $|\vec{p}_{+}|/E_{\rm beam} > 0.75 ~~\textnormal{and}~~ 
|\vec{p}_{-}|/E_{\rm beam} > 0.5 ~~\textnormal{or}$ \\
$|\vec{p}_{-}|/E_{\rm beam} > 0.75 ~~\textnormal{and}~~ 
|\vec{p}_{+}|/E_{\rm beam} > 0.5$\,, 
\item $\xi_{\rm max}^{3d} =30^{\circ}$\,.
\end{itemize}
Here $E_{\rm min}$  is the energy threshold for the final-state electron/positron, 
$\theta_{\pm}$ are the electron/positron polar angles and $\xi_{max}$ is
the maximum allowed polar angle acollinearity:
\begin{eqnarray}
\xi = \left|\theta_{+}+\theta_{-} -180^{\circ}\right| ,
\end{eqnarray}
and $\xi_{\rm max}^{3d}$ is the 
maximum allowed three dimensional a\-col\-linea\-ri\-ty: 
\begin{eqnarray}
\xi^{3d} = \left|\arccos\left(\frac{p_{+}\cdot p_{-}}{(|\vec{p}_{-}|
|\vec{p}_{+}|}\right) \times \frac{180^{\circ}}{\pi}  -180^{\circ}\right|.
\end{eqnarray}
For $e^{+}e^{-}\rightarrow e^{+}e^{-}\mu^{+}\mu^{-}$,  
cuts are applied only to the $e^{+}e^{-}$ pair.
In the case of  $e^{+}e^{-}\rightarrow 
e^{+}e^{-}e^{+}e^{-}$, all possible $e^{\pm}e^{\mp}$ combinations are 
checked and if at least one pair fulfils the cuts the 
event is accepted. 
%Finally, the following masses have been used:
%$m_e= 0.510998902 \times 10^{-3} ~\textnormal{GeV}, 
%~m_{\mu}= 0.105658369 ~\textnormal{GeV}, ~m_{\tau}= 1.77699 ~\textnormal{GeV}$.


%The pair corrections are very small.
%%NEW RESULTS. 
At KLOE the electron pair corrections contribute about $3\times 10^{-3}$ and at BaBar
about $1\times 10^{-3}$, while all the other contributions of pair production are even smaller. Like in 
small-angle Bhabha scattering at LEP/SLC the pair corrections \cite{Montagna:1998vb} 
are largely dominated by the electron pair contribution.
%This corresponds to earlier estimates. CITE ??? \cite{earlier-pair-estimates?}.


