%Title: Heavy-Flavor Contribution to Bhabha Scattering
%Author: R. Bonciani, A. Ferroglia, A. A. Penin
%Published: <EM> Phys.Rev.Lett. </EM> <B> 100 </B> (2008) 131601. 
\documentclass[twocolumn,showpacs,preprintnumbers,amsmath,amssymb,floatfix]{revtex4}
%\documentclass[11pt]{article}
\usepackage{epsfig,amsfonts}
\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math
%
\DeclareMathSymbol{\lesssim}      {\mathrel}{AMSa}{"2E}
\DeclareMathSymbol{\gtrsim}       {\mathrel}{AMSa}{"26}
%
\topmargin=-40pt


\begin{document}

\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bfm}[1]{\mbox{\boldmath$#1$}}
\newcommand{\bff}[1]{\mbox{\scriptsize\boldmath${#1}$}}
\newcommand{\al}{\alpha}
\newcommand{\bt}{\beta}
\newcommand{\lm}{\lambda}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\gm}{\gamma}
\newcommand{\Gm}{\Gamma}
\newcommand{\dl}{\delta}
\newcommand{\Dl}{\Delta}
\newcommand{\ep}{\epsilon}
\newcommand{\vep}{\varepsilon}
\newcommand{\kp}{\kappa}
\newcommand{\Lm}{\Lambda}
\newcommand{\om}{\omega}
\newcommand{\pa}{\partial}
\newcommand{\nn}{\nonumber}
\newcommand{\dd}{\mbox{d}}
\newcommand{\grtsim}{\mbox{\raisebox{-3pt}{$\stackrel{>}{\sim}$}}}
\newcommand{\lessim}{\mbox{\raisebox{-3pt}{$\stackrel{<}{\sim}$}}}
\newcommand{\uk}{\underline{k}}
\newcommand{\gsim}{\;\rlap{\lower 3.5 pt \hbox{$\mathchar \sim$}} \raise 1pt \hbox {$>$}\;}
\newcommand{\lsim}{\;\rlap{\lower 3.5 pt \hbox{$\mathchar \sim$}} \raise 1pt \hbox {$<$}\;}
\newcommand{\Li}{\mbox{Li}}
\newcommand{\bc}{\begin{center}}
\newcommand{\ec}{\end{center}}


\title{Heavy-Flavor Contribution to Bhabha Scattering}


\author{R.~Bonciani}
\affiliation{Departamento de F\'{\i}sica Te\`orica, IFIC, CSIC -- Universidad
de Valencia, E-46071 Valencia, Spain}
\author{A.~Ferroglia}
\affiliation{Institut f{\"u}r Theoretische Physik, Universit{\"a}t Z\"urich,
CH-8057 Zurich, Switzerland}
\author{A.~A.~Penin}
\affiliation{Department of Physics, University Of Alberta,
Edmonton, AB T6G 2J1, Canada\\
Institute for Nuclear Research of Russian Academy of Sciences,
117312 Moscow, Russia}


\begin{abstract}
  We evaluate the last missing piece of the two-loop QED corrections to the
  high-energy electron-positron scattering cross section originated from
  the vacuum polarization by heavy fermions.
  The calculation is performed within a new approach applicable
  to a wide class of perturbative problems with  mass hierarchy.
  The result is crucial  for the  high-precision physics program at existing
  and future $e^+e^-$ colliders.

\end{abstract}

\pacs{11.15.Bt, 12.20.Ds}
%\keywords{Feynman diagrams, multi-loop calculation, Bhabha scattering}
%
\maketitle
%


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

Electron-positron  scattering at high energy has been originally
studied  by J.H. Bhabha  almost a century ago \cite{Bha} and
nowadays plays an ultimate role in particle phenomenology. It
provides a very efficient tool for luminosity determination at
electron-positron colliders and thus mediates the process of
extracting physical information from the raw experimental data
\cite{Jad}. The accuracy of theoretical predictions for the cross
section of {\it Bhabha  scattering} directly affects our quest for
``new physics'' at operating and planned $e^+e^-$ colliders. For
example, the future International Linear Collider (ILC) requires the
accuracy  of 0.1 permille for the small angle  scattering to provide
the most accurate  tests of the Standard Model of electroweak
interactions \cite{Agu}. At the same time the low-energy experiments
\cite{Alo,Eid} require about one permille accuracy of the cross
section at large scattering angles to  reduce the hadronic
uncertainty in  the theoretical prediction for the muon anomalous
magnetic moment, which is a renowned probe  of the physics beyond
the Standard Model \cite{CzaMar}. Thus, insufficient theoretical
accuracy may jeopardize the whole high-precision physics program at
$e^+e^-$ colliders. To match the  experimental needs, the two-loop
QED corrections must  be included into the theoretical analysis.

The calculation of the radiative corrections to the Bhabha cross
section is among the classical problems of perturbative quantum
field theory with a long history  and remarkable efforts have been
devoted to the study of the two-loop QED contributions.   The
corrections due to the vacuum polarization by virtual
electron-positron pairs have been evaluated in \cite{BFMR} including
the full dependence on the electron mass $m_e$. Pure photonic
corrections have been obtained in Ref.~\cite{Pen} (see also
\cite{GTV}) in the leading order of the small electron mass
expansion through the {\it infrared matching} to the massless case
\cite{BDG}. This result is sufficient for all the phenomenological
applications and has been recently confirmed within a slightly
different framework \cite{BecMel}.
%The above corrections have already been taken into account in the
%Monte-Carlo event generators  \cite{Car}.
For the vacuum polarization by virtual heavy-flavor fermion-antifermion pairs
the result is available only in the limit
when the heavy fermion mass $m_f$  is much less than the
center-of-mass energy $\sqrt{s}$ \cite{BecMel,Act}.
This condition, however, does
not hold for $\tau$-lepton, $c$ and $b$-quark in the practically
interesting energy range of about a few GeV as well as for the top
quark at the typical ILC energies $500~{\rm GeV}\lsim\sqrt{s}\lsim
1000~{\rm GeV}$. The result for an arbitrary value of $m_f^2/s$
constitutes the last missing part of the
QED two-loop corrections. This  two-loop scattering problem
with two different  masses is a real challenge
for modern methods of perturbative
quantum field theory.

In this Letter we present  the  heavy-flavor
contribution retaining the full dependence on $m_f^2/s$. We first outline
a new approach based on the general theory of infrared divergences, which
makes the whole calculation possible.
Then we present the complete numerical result for the
correction in  two phenomenologically most important  cases and
give the  analytical expression  in two limits
$m_f^2 \ll s$ and $m_f^2 \gg s$. For completeness, we include
the mixed QED-QCD two-loop corrections in the numerical analysis
of the heavy-quark contribution. 

The perturbative expansion for the Bhabha cross section in the fine
structure constant $\al$ is defined
as follows ${\rm d}\sigma=\sum_{n=0}^\infty \left({\al/\pi}\right)^n
\delta^{(n)}{\rm d}\sigma^{(0)}$, where $\delta^{(n)}$ is the
correction factor, with $\delta^{(0)}=1$.
In the phenomenologically interesting kinematical region
all the Mandelstam invariants are large relative to the electron mass
$s,~t,~u\gg m_e^2$.  Therefore all the
terms suppressed by $m_e$ can be neglected. The leading order
differential cross section  in the small electron mass approximation
takes the form
\begin{equation}
{{\rm d}\sigma^{(0)}\over{\rm d}\Omega} ={\al^2\over
s}\left({1-x+x^2\over x}\right)^2+{\mathcal O}(m_e^2/s)\,,
\label{losig}
\end{equation}
where $x=(1-\cos\theta)/2$ and $\theta$ is the scattering angle.
Virtual corrections are infrared divergent. These {\it soft} divergencies are
canceled in the physical inclusive cross section when one adds the
photonic bremsstrahlung \cite{Kin}.  We regulate the soft divergencies
by dimensional regularization with $d=4-2\vep$ space-time dimensions.
The standard approach
to deal with the bremsstrahlung is to split it into a soft part,
which accounts for the emission of the photons with the energy
below a given cutoff $\vep_{cut}\ll m_e$, and a hard part
corresponding to the emission of the photons with the energy
larger than $\vep_{cut}$.  The infrared finite hard part is then
computed numerically using Monte-Carlo methods with physical cuts
dictated by the experimental setup.  At the same time, the soft
part is computed analytically and combined with the virtual
corrections ensuring the cancellation of the soft divergencies.
Thus we consider the finite sum of two terms
$\delta^{(2)}=\delta^{(2)}_{vv}+\delta^{(1)}_{s}\delta^{(1)}$
corresponding to the two-loop virtual correction \footnote{We do not
  consider the trivial correction given by two heavy-fermion loop
  insertions.} and the one-loop correction to the single soft photon
emission, which factorizes into the product of the first-order
contributions \cite{YFS}.
\begin{figure}[t]
\bc
\begin{picture}(0,0)%
\includegraphics[height=5.0cm]{CS_muon_60_LAT_log}%
\end{picture}%
\begingroup
\setlength{\unitlength}{0.0130bp}%
\begin{picture}(18000,10800)(0,0)%
\put(2200,3253){\makebox(0,0)[r]{\strut{} 0.1}}%
\put(2200,7485){\makebox(0,0)[r]{\strut{} 1}}%
\put(2475,1429){\makebox(0,0){\strut{} 0.1}}%
\put(14974,1429){\makebox(0,0){\strut{} 1}}%
\put(9825,275){\makebox(0,0){\strut{} $\sqrt{s}$ (GeV)}}%
\put(3465,7485){\makebox(0,0)[l]{\strut{} $\theta=60^{\circ}$ }}%
\put(5766,9170){\makebox(0,0)[l]{\strut{} large mass exp. }}%
\put(8439,5756){\makebox(0,0)[l]{\strut{} exact }}%
\put(10393,3059){\makebox(0,0)[l]{\strut{} small mass exp. }}%
\put(800,6114){\rotatebox{90}{\makebox(0,0){\strut{} $10^4 \cdot d\sigma^{(2)}/d\sigma^{(0)}$ }}}%
\end{picture}%
\endgroup
\caption{\small{ Two-loop corrections to the Bhabha scattering
differential cross section at $\theta=60^{\circ}$ due to a closed
loop of muon ($m_\mu=113$~MeV)} \label{fig6} \vspace*{-6mm}} \ec
\end{figure}
In the small electron mass limit the two-loop correction
can be decomposed according to the asymptotic dependence on $m_e$
\begin{equation}
\delta^{(2)}= \delta^{(2)}_{1}\ln\left({s\over m_e^2}\right)
+\delta^{(2)}_{0}+{\cal O}(m_e^2/s)\, . \label{2lexp}
\end{equation}
The logarithmic term in Eq.~(\ref{2lexp}) is a remnant of the {\it
  collinear} divergence which is regulated by the electron mass.
Now we come to the crucial point of our approach.
The collinear divergencies and, hence, the singular dependence of the
corrections on $m_e$ have a peculiar structure which was
extensively studied in the context of QCD. In particular, the
collinear divergencies factorize and can be absorbed in the
external field renormalization \cite{FreTay}. This property is
crucial for our analysis because it allows us to perform the most
difficult part of the calculation with a strictly massless
electron.  Indeed, due to the factorization, the singular
dependence on $m_e$ is the same for the Bhabha amplitude and (the
square of) the vector form factor \cite{Pen}. This attributes the
total logarithmic correction ot the amplitude
to  the one-particle reducible diagram given by the product
of the one-loop contributions.  In calculating
the cross section one must take into account also the contributions coming
from the interference of the one-loop corrections to the amplitude
and the corresponding soft emission; both contributions have a
factorized form. Thus it is straightforward to obtain the
coefficient of the logarithmic term in Eq.~(\ref{2lexp}) which reads
$\delta^{(2)}_{1} = \left[ 4\ln\left({\vep_{cut}/\vep}\right) + 3
\right] \delta^{(1)}$, where $\vep=\sqrt{s}/2$.
\begin{figure}[t]
\bc
\begin{picture}(0,0)%
\includegraphics[height=5.0cm]{CS_DAPtaucb_LAT}%
\end{picture}%
\begingroup
\setlength{\unitlength}{0.0130bp}%
\begin{picture}(18000,10800)(0,0)%
\put(2200,1979){\makebox(0,0)[r]{\strut{}-2}}%
\put(2200,3554){\makebox(0,0)[r]{\strut{} 0}}%
\put(2200,5130){\makebox(0,0)[r]{\strut{} 2}}%
\put(2200,6705){\makebox(0,0)[r]{\strut{} 4}}%
\put(2200,8281){\makebox(0,0)[r]{\strut{} 6}}%
\put(2200,9856){\makebox(0,0)[r]{\strut{} 8}}%
\put(2475,1429){\makebox(0,0){\strut{} 0}}%
\put(4108,1429){\makebox(0,0){\strut{} 20}}%
\put(5742,1429){\makebox(0,0){\strut{} 40}}%
\put(7375,1429){\makebox(0,0){\strut{} 60}}%
\put(9008,1429){\makebox(0,0){\strut{} 80}}%
\put(10642,1429){\makebox(0,0){\strut{} 100}}%
\put(12275,1429){\makebox(0,0){\strut{} 120}}%
\put(13908,1429){\makebox(0,0){\strut{} 140}}%
\put(15542,1429){\makebox(0,0){\strut{} 160}}%
\put(17175,1429){\makebox(0,0){\strut{} 180}}%
\put(9825,275){\makebox(0,0){\strut{} $\theta$ }}%
\put(3292,9226){\makebox(0,0)[l]{\strut{} $\sqrt{s}=1$~GeV }}%
\put(11305,9226){\makebox(0,0)[l]{\strut{} $c$ }}%
\put(8192,4500){\makebox(0,0)[l]{\strut{} $\tau$ }}%
\put(11458,3003){\makebox(0,0)[l]{\strut{} $b$ }}%
\put(1100,6114){\rotatebox{90}{\makebox(0,0){\strut{} $10^5 \cdot d\sigma^{(2)}/d\sigma^{(0)}$ }}}%
\end{picture}%
\endgroup
%\vspace*{-3mm}
\caption{\small{Two-loop corrections to the Bhabha scattering
differential cross section at $\sqrt{s}=1$~GeV due to a closed loop
of $\tau$-lepton ($m_{\tau}=1.7$~GeV), $c$-quark ($m_c=1.25$~GeV)
and $b$-quark ($m_b=4.7$~GeV).} \label{fig2}\vspace*{-0mm}}
\ec
\end{figure}
At the same time,
the sum of the one-particle irreducible diagrams has a regular
behavior in the small electron mass limit and can be computed with
$m_e=0$.  Let us emphasize that this property in general holds only
for the sum of the diagrams.  The individual diagrams computed in a
covariant gauge do exhibit collinear divergencies for $m_e=0$.
This however does not pose any additional problem since we work in
dimensional regularization. In this case the collinear divergencies
show up as the extra poles in $\vep$ which are not related to the
soft emission and disappear in the sum of the one-particle
irreducible diagrams. Thus all the ``true'' two-loop diagrams
contribute only to the non logarithmic term in Eq.~(\ref{2lexp})
which can be written as follows \bea
 \delta^{(2)}_{0} & = & - 4\left[1+\ln\left({1-x\over x}\right)
\right]\ln\left({\vep_{cut}\over\vep}\right){\delta^{(1)}} \nn\\
& & + \left( {x\over 1-x+x^2} \right)^2 Q_f^2 N_cf(\rho,x)\, ,
\label{del20} \eea where the first term is determined by the soft
emission, $f(\rho,x)$ is a function of two dimensionless variables:
$\rho=m_f^2/s$ and $x=-t/s$, and we explicitly factored out  the
heavy-flavor electric charge $Q_f$ and the number of colors $N_c$,
which  is 1 for leptons and 3 for quarks. The two-loop problem with
massless electron falls it the same complexity class as the one
considered in Ref.~\cite{BFMR}. The most difficult part of the
calculation is given by the two-loop box diagrams. We use the
Laporta Algorithm \cite{Laporta} to identify the set of so-called
master integrals which  are subsequently evaluated by means of the
differential equation method \cite{diffeq}. Finally, the result for
the function $f(\rho,x)$ can be expressed in closed analytic form in
terms of  mathematical functions known as one and two-dimensional
generalized harmonic polylogarithms of maximum weight three. In this
Letter we focus on the numerical impact of the correction in the
phenomenologically important cases and postpone the details of the
analysis and the complete analytical result for a future
publication. However, the  analytical structure of the result
becomes much simpler in the small-mass limit $m_f^2\ll s$ and the
large-mass limit $m_f^2\gg s$ where it is given by the  series
$f(\rho,x)=\sum_{n=0}^\infty\rho^nf_n(\rho,x)$ and
$f(\rho,x)=\sum_{n=0}^\infty\rho^{-n}\bar{f}_n(\rho,x)$,
respectively. The series are asymptotical and the coefficients
$f_n(\rho,x)$ and $\bar{f}_n(\rho,x)$ are in general logarithmic
functions of $\rho$. With the exact result at hand we can compute
any term of these expansions. The leading order coefficient
${f}_0(\rho,x)$ has been computed in \cite{BecMel,Act} and we
confirm this result. On the other hand the leading order coefficient
$\bar{f}_0(\rho,x)$ vanishes due to on-shell renormalization of the
fine structure constant. Thus we present the results for the first
unknown  terms of the expansions which read
\begin{eqnarray}
& & \hspace{-4mm} f_1(\rho,x) \!= \!\frac{1}{x^3}\bigg[\!(4 x^5- 6
x^4+ 8 x^3- 8 x^2+ 6 x-4) \left[
\ln^2 \left(\rho\right) \right.\nn \\
& &  \hspace{0mm} \left.+ 4 \Li_{2}(x) \! + \! 12 \right] \! + \!
\ln\left(\rho\right) \! \bigl[
   \!  (2 x^5\! \! - \! 7 x^4\! \! + \! 10 x^3\! \! - \! 10 x^2\! \!
   + \! 7 x\! -\! 2) \nn\\
& &  \hspace{0mm}
- (4x^5-6x^4+10x^3-14x^2+12x-8) \ln(x) \nn\\
& &  \hspace{0mm} + (4x^5-6x^4+8x^3-8x^2+6x-4) \ln(1 - x)
\bigr] \nn \\
& &  \hspace{0mm} - (40 x^5-54 x^4+50 x^3-17 x^2-12
x+8) \zeta(2)  \nn \\
& &  \hspace{0mm}
+ (12 x^5- 17 x^4+ 18 x^3- 18 x^2+ 17 x-12)\ln(1 - x) \nn \\
& &  \hspace{0mm}
+ \frac{1}{2}(8x^5-14x^4+15x^3-15x^2+14x-8) \ln^2 (1 \!  \! -  \! x) \nn \\
& &  \hspace{0mm}
- (12 x^5-21 x^4+26 x^3-26 x^2+21 x-14) \ln(x)   \nn \\
& &  \hspace{0mm}
+ (8x^5-10x^4+8x^3-x^2-6x+4) \ln(x) \ln(1 \!  \! -  \! x)  \nn \\
& &  \hspace{0mm} - \frac{1}{2}(8x^5-7x^4-x^3+12x^2-15x+8)
\ln^2(x) \nn \\
& &+\left(Q_f^2-1\right)\left[ (6 \!-\! 9x\! +\!
12x^2\! -\! 12x^3 + \!9x^4 \!-\! 6x^5) \ln\left(\rho \right)
\right.\nn\\
&&\left.- \left(6 - 9 x + 9 x^2- 3x^3\right) \ln\left( x\right)
\right] \bigg]\, . \label{f1}
\end{eqnarray}
%
%
\begin{eqnarray}
& & \hspace{-3mm}
\bar{f}_1(\rho,x) =
  \frac{5 x^3 \! \! - \! 22 x^2 \! \! + \! 22 x \! - \! 5}{30 x} \ln\left(\rho \right)
+  \! \frac{955 x^3 \! \! - \! 3926 x^2 }{1350 x}
  \nn\\
& &  \hspace{0mm}
+  \! \frac{3926 x \! - \! 955}{1350 x}
-  \frac{12 x^3 \! \! - \! 19 x^2 \! \! + \! 14 x \! - \! 6}{10 x}\zeta(2)
\nn \\
& &  \hspace{0mm}
+  \! \frac{3 x^3 \! \! + \! x^2 \! \! - \! x \! - \! 3}{30 x}\ln{(1 \! - \! x)}
%\nn\\
%& &  \hspace{0mm}
\! + \frac{2 x^3 \! \! -\! 5 x^2 \! \! +\! 5 x\! -\! 2}{20 x} \ln^2 {(1 \! - \! x)}\nn\\
& &  \hspace{0mm}
-  \frac{20 x^3-78 x^2+93 x-58}{90 x}  \ln{(x)}
%\nn\\
%& &  \hspace{0mm}
-  \frac{6 x^2-x-4}{60} \ln^2 {(x)} \nn\\
& &   \hspace{0mm}
+ \frac{12 x^3-19 x^2+14 x-6}{30 x}  \ln{(x)}\ln{(1- x)} \nn\\
& &  \hspace{0mm} +   \frac{4x^3-8 x^2+8 x-4}{5 x}  \mbox{Li}_2(x)
\nn\\
&&-\left(Q_f^2-1\right) \frac{41x^3 -82 x^2 + 82 x -41}{54 x }\,.
\label{f1bar}
\end{eqnarray}
In Fig.~\ref{fig6}  the second-order muon contribution
$(\al/\pi)^2\delta^{(2)}$ is plotted as function of $\sqrt{s}$
in the threshold region around $2m_\mu\approx 216$~MeV at the
fixed scattering angle $\theta=60^{\circ}$. All the terms involving
a power of the logarithm $\ln(\vep_{cut}/\vep)$ are excluded from
the numerical estimates because the corresponding contribution
critically depends on the event selection algorithm and cannot be
unambiguously estimated without imposing specific cuts on the photon
bremsstrahlung. We present the exact result along with the
first nonvanishing  terms of the small and large mass expansions.
The asymptotic results fail to approximate the exact
one with 10\% accuracy in the interval $1.7\lsim \sqrt{s}/m_f\lsim
10$. In particular they are completely useless for the description
of the top quark contribution at the typical energies of the ILC.
Both expansions break down near the threshold $\sqrt{s}\sim 2m_f$.

Let us  discuss the numerical impact of the corrections on the
luminosity determination at present and future electron-positron
colliders. As the first application, we consider the Bhabha
scattering at $\sqrt{s}=1$~GeV. This is the center-of-mass energy of
the KLOE experiment at DA$\Phi$NE, which plays crucial role in
determination of the hadronic vacuum polarization contribution to
the muon anomalous magnetic moment. The second-order contributions
of $\tau$-lepton, $c$, and $b$-quark to the differential cross
section are plotted as functions of the scattering angle in
Fig.~\ref{fig2}. Note that the energy under consideration is
sufficiently below the quarkonium threshold and the heavy-quarks can
be treated perturbatively. Thus we included into the numerical
estimates the mixed ${\cal O}(\alpha\alpha_s)$  corrections due to
the virtual gluon exchange inside the heavy quark loop.
 The contributions of the $\tau$-lepton, $c$, and
$b$-quark are suppressed with respect to the muon at least by one
order of magnitude ({\it c.f.} Fig.~\ref{fig6}). Thus the total
heavy-flavor contribution is dominated by the muon loop and reaches
about 1/2 permille at large scattering angles. Let us now discuss
the top quark contribution to Bhabha scattering at high energies
characteristic to the ILC. The large electroweak decay width of the
top quark serves as an infrared regulator and suppresses the
hadronization effects. Thus our perturbative result is  applicable
also at the energy about and above the top-antitop threshold. In
Fig.~\ref{fig5} the contribution of top quark is plotted as function
of the scattering angle at $\sqrt{s}=500$~GeV. At large angles it
reaches about 1/4 permille in magnitude. Note that  the corrections
are of the same magnitude  as the nonlogarithmic pure photonic terms
\cite{Pen} if $\theta$ is not close to $0^\circ$ or $180^\circ$ and
vanish for strictly forward scattering.

\begin{figure}[t]
\bc
\begin{picture}(0,0)%
\includegraphics[height=5.0cm]{CS_ILCtop_LAT}%
\end{picture}%
\begingroup
\setlength{\unitlength}{0.0130bp}%
\begin{picture}(18000,10800)(0,0)%
\put(2475,1979){\makebox(0,0)[r]{\strut{} 0}}%
\put(2475,3358){\makebox(0,0)[r]{\strut{} 1}}%
\put(2475,4736){\makebox(0,0)[r]{\strut{} 2}}%
\put(2475,6115){\makebox(0,0)[r]{\strut{} 3}}%
\put(2475,7493){\makebox(0,0)[r]{\strut{} 4}}%
\put(2475,8872){\makebox(0,0)[r]{\strut{} 5}}%
\put(2475,10250){\makebox(0,0)[r]{\strut{} 6}}%
\put(2750,1429){\makebox(0,0){\strut{} 0}}%
\put(4353,1429){\makebox(0,0){\strut{} 20}}%
\put(5956,1429){\makebox(0,0){\strut{} 40}}%
\put(7558,1429){\makebox(0,0){\strut{} 60}}%
\put(9161,1429){\makebox(0,0){\strut{} 80}}%
\put(10764,1429){\makebox(0,0){\strut{} 100}}%
\put(12367,1429){\makebox(0,0){\strut{} 120}}%
\put(13969,1429){\makebox(0,0){\strut{} 140}}%
\put(15572,1429){\makebox(0,0){\strut{} 160}}%
\put(17175,1429){\makebox(0,0){\strut{} 180}}%
\put(9962,275){\makebox(0,0){\strut{} $\theta$ }}%
\put(11565,4047){\makebox(0,0)[l]{\strut{} $\sqrt{s}=500$~GeV }}%
%\put(11565,3068){\makebox(0,0)[l]{\strut{} $m_t=170.9$~GeV }}%
\put(1100,6114){\rotatebox{90}{\makebox(0,0){\strut{}
$10^4 \cdot d\sigma^{(2)}/d\sigma^{(0)}$ }}}%
\end{picture}%
\endgroup
%\vspace*{-3mm}
\caption{\small{Two-loop corrections to the Bhabha scattering
differential cross section at $\sqrt{s}=500$~GeV due to a closed
loop of $t$-quark ($m_t=170.9$~GeV). } \label{fig5}\vspace*{-0mm}}
%\vspace*{-6mm}
\ec
\end{figure}

To this moment we have not mentioned the  contribution
of light quarks.
In principle, it requires a special treatment because of hadronization
effects (see {\it e.g.} \cite{KKKS}). However, this contribution
can be estimated by  na\"ive use of the perturbative result
with the constituent quark masses. Such estimates are normally
in rough agreement with the result of the rigorous analysis
based on the dispersion approach. Adopting $m_d=m_u=130$~MeV dictated
by the hadronic  threshold, we find the contribution of $d$ and
$u$-quark to be  below one permille for $\sqrt{s}=1$~GeV.
It is not clear whether a more refined calculation  is necessary
at the present level of experimental accuracy. Thus, we have full
control over the two-loop corrections to Bhabha scattering at
low energies  where the QED contributions
are dominant. The accuracy of the luminosity determination
is now restricted only by the precision of the
Monte-Carlo event generators for the hard photon
and electron-positron pair emission \cite{Car}.
To achieve such an accuracy for the large angle
scattering at the ILC one has to perform
more careful  analysis of the hadronic contributions and
take into account also yet unknown two-loop electroweak corrections.

To conclude, we have derived the two-loop radiative corrections to
Bhabha scattering due to  vacuum polarization by the virtual heavy-flavor
fermion-antifermion pairs which completes the QED analysis in two loops.
This solves one of  the current front-line  problems
in  perturbative quantum field theory.
The result is crucial for the high-precision physics
at electron-positron colliders. It removes the last
piece of pure theoretical uncertainty in luminosity
determination at low-energy accelerators  and gives the  proper
account for the top quark effects at the ILC.
We believe that our approach will  be very useful
for the  high-order perturbative analysis in a wide class of the
processes with mass hierarchy.

%\vspace*{-5mm}
\begin{acknowledgements}
R.~B. would like to thank the Galileo Galilei Institute for Theoretical
Physics for the hospitality and the INFN for partial support during the
completion of this work.
The work of R.~B. was partially supported by
Ministerio de Educaci\'on y Ciencia (MEC) under grant FPA2004-00996,
Generalitat Valenciana under grant GVACOMP2007-156, European
Commission under the grant MRTN-CT-2006-035482 (FLAVIAnet), and
MEC-INFN agreement. The work of A.~F. was supported in part by the
Swiss National Science Foundation (SNF) under contract
200020-117602. The work of A.~P. was supported in part by BMBF Grant
No.\ 05HT4VKA/3 and Sonderforschungsbereich Transregio 9.
\end{acknowledgements}
%




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



\begin{thebibliography}{99}

\bibitem{Bha}
J.~H.~Bhabha, Proc.Roy.Soc. {\bf 154A}, 195 (1935).

\bibitem{Jad}
S. Jadach {\it et al.} in G. Altarelli, T. Sj\"ostrand and F.
Zwirner (eds.), {\it Physics at LEP2}, CERN-96-01.

\bibitem{Agu}
J.~A.~Aguilar-Saavedra {\it et al.}  [ECFA/DESY LC Physics Working Group],
%``TESLA Technical Design Report Part III: Physics at an e+e- Linear
%Collider,''
arXiv:hep-ph/0106315.
%%CITATION = HEP-PH/0106315;%%

\bibitem{Alo}
A.~Aloisio {\it et al.}  [KLOE Collaboration],
%``Measurement of sigma(e+ e- --> pi+ pi- gamma) and extraction of  sigma(e+
%e- --> pi+ pi-) below 1-GeV with the KLOE detector,''
Phys.\ Lett.\  B {\bf 606}, 12 (2005).
%%CITATION = PHLTA,B606,12;%%

\bibitem{Eid}
S.~Eidelman  [CMD-2 and SND Collaborations],
%``Recent Results From Cmd-2 And Snd Experiments At Vepp-2m,''
PoS {\bf HEP2005}, 114 (2006).
%%CITATION = POSCI,HEP2005,114;%%



\bibitem{CzaMar}
  A.~Czarnecki and W.~J.~Marciano,
  %``The muon anomalous magnetic moment: A harbinger for 'new physics',''
  Phys.\ Rev.\  D {\bf 64}, 013014 (2001).
  %%CITATION = PHRVA,D64,013014;%%
%
\bibitem{BFMR}
R.~Bonciani, A.~Ferroglia, P.~Mastrolia, E.~Remiddi, and J.~J.~van der Bij,
%``Two-loop N(F) = 1 QED Bhabha scattering differential cross section,''
Nucl.\ Phys.\  B {\bf 701} 121 (2004);
%%CITATION = NUPHA,B701,121;%%
%
%\bibitem{Bon}
%R. Bonciani, A. Ferroglia, P. Mastrolia, E. Remiddi, and J.J. van der Bij,
{\it ibid.}\ {\bf B716},  280 (2005).
%%CITATION = NUPHA,B716,280;%%

\bibitem{Pen}
A.~A.~Penin, Phys.\ Rev.\ Lett.\ {\bf 95}, 010408 (2005); Nucl.\
Phys.\  {\bf B734}, 185 (2006).
%%CITATION = NUPHA,B734,185;%%
%%CITATION = PRLTA,95,010408;%%


\bibitem{GTV}
E.~W.~N.~Glover, J.~B.~Tausk, and J.~J.~Van der Bij,
%``Second order contributions to elastic large-angle Bhabha scattering,''
Phys.\ Lett.\  B {\bf 516} 33 (2001).
%%CITATION = PHLTA,B516,33;%%


\bibitem{BDG}
Z.~Bern, L.~Dixon, and A.~Ghinculov, Phys.\ Rev.\ D {\bf 63}, 053007
(2001).
%%CITATION = HEP-PH 0010075;%%


\bibitem{BecMel}  T.~Becher and K.~Melnikov,
%``Two-loop QED corrections to Bhabha scattering,''
JHEP {\bf 0706}, 084 (2007).
%%CITATION = JHEPA,0706,084;%%


\bibitem{Act}
S.~Actis, M.~Czakon, J.~Gluza, and T.~Riemann,
%``Two-Loop Fermionic Corrections to Massive Bhabha Scattering,''
Nucl.\ Phys.\ {\bf B786}, 26 (2007).
%%CITATION = ARXIV:0704.2400;%%


\bibitem{Kin}
T.~Kinoshita,
%``Mass Singularities Of Feynman Amplitudes,''
J.\ Math.\ Phys.\  {\bf 3}, 650 (1962);
%%CITATION = JMAPA,3,650;%%
%
%\bibitem{LeeNau}
T.~D.~Lee and M.~Nauenberg,
%``Degenerate Systems and Mass Singularities,''
Phys.\ Rev.\  {\bf 133}, B1549 (1964).
%%CITATION = PHRVA,133,B1549;%%

\bibitem{YFS}
D.~R.~Yennie, S.~C.~Frautschi, and H.~Suura,
%``The infrared divergence phenomena and high-energy processes,''
Annals Phys.\  {\bf 13}, 379 (1961).
%%CITATION = APNYA,13,379;%%


\bibitem{FreTay}
J.~Frenkel and J.~C.~Taylor,
%``Exponentiation Of Leading Infrared Divergences In Massless Yang-Mills
%Theories,''
Nucl.\ Phys.\  B {\bf 116}, 185 (1976).
%%CITATION = NUPHA,B116,185;%%

\bibitem{Laporta}
S.~Laporta and E.~Remiddi,
%``The analytical value of the electron (g-2) at order alpha~3 in QED,''
Phys.\ Lett.\ B {\bf 379} (1996) 283;
%%CITATION = HEP-PH 9602417;%%
S.~Laporta,
%``High-precision calculation of multi-loop Feynman integrals by  difference
%equations,''
Int.\ J.\ Mod.\ Phys.\  A {\bf 15} (2000) 5087
[arXiv:hep-ph/0102033].
%%CITATION = IMPAE,A15,5087;%%

\bibitem{diffeq}
A.~V.~Kotikov,
%``Differential Equations Method: New Technique For Massive Feynman Diagrams
%Calculation,''
Phys. Lett. B {\bf 254} (1991) 158;
E.~Remiddi,
%``Differential equations for Feynman graph amplitudes,''
Nuovo Cim.\ A {\bf 110} (1997) 1435.

\bibitem{KKKS}
B.~A.~Kniehl, M.~Krawczyk, J.~H.~Kuhn, and R.~G.~Stuart,
%``HADRONIC CONTRIBUTIONS TO O (alpha**2) RADIATIVE CORRECTIONS IN e+ e-
%ANNIHILATION,''
Phys.\ Lett.\  B {\bf 209}, 337 (1988).
%%CITATION = PHLTA,B209,337;%%



\bibitem{Car}
G.~Balossini, C.~M.~Carloni Calame, G.~Montagna, O.~Nicrosini, and F.~Piccinini,
%``Matching perturbative and parton shower corrections to Bhabha process at
%flavour factories,''
Nucl.\ Phys.\  B {\bf 758}, 227 (2006).
%%CITATION = NUPHA,B758,227;%%

\end{thebibliography}


\end{document}
