%Title: Two-Loop Formfactors in Theories with Mass Gap and Z-Boson Production. 
%Author: A. Kotikov, J.H. K&uuml;hn and O. Veretin
%Published: <EM> Nucl.Phys. </EM> <B> B788 </B> (2008) 47-62. 
%hep-ph/0703013
%/

%\documentstyle[12pt,axodraw]{article}
\documentclass[12pt,a4paper]{article}
\usepackage{axodraw}
\usepackage{epsfig}
\topmargin 0pt
\textwidth  178mm
\textheight 230mm
\hoffset   -15mm
\voffset   -15mm
\oddsidemargin 5mm
\pagestyle{plain}
\def\baselinestretch{1.25}

%\documentstyle[12pt]{article,axodraw}
%\textwidth 16cm
%\textheight 23cm
%\topmargin -1cm
%\oddsidemargin 6mm
%\evensidemargin 6mm

\newcommand{\z}{&\hspace*{-8pt}}
\newcommand{\ep}{\varepsilon}
\newcommand{\vp}{\varphi}
\newcommand{\am}{\overline \alpha_{\mu}}
\newcommand{\as}{ \overline \alpha_{s}}
\newcommand{\MS}{\overline{MS}}
\newcommand{\Li}{{\rm Li}}
\newcommand{\Si}{{\rm S}}
\newcommand{\Ls}{{\rm Ls}}
\newcommand{\T}{{\rm T}}
%\newcommand{\C}{{\rm C}}
\newcommand{\Cl}{{\rm Cl}}
%\newcommand{\G}{{\rm G}}
\newcommand{\au}{{\overline a}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\ps}{\varphi}

\newcommand \lz[1] {\log^{#1}{z}} 
\newcommand \Oz[1] { {\cal O}\left( \frac{1}{z^{#1}} \right) }


\def \LOG{L}

%\input{tcilatex}

\begin{document}
\title{\vskip-3cm{\baselineskip14pt
\centerline{\normalsize\hfill TTP07-05}
\centerline{\normalsize\hfill SFB/CPP-07-07}
\centerline{\normalsize\hfill DESY 07-024}
%\centerline{\normalsize\hfill hep-ph/yymmxxx}
%\centerline{\normalsize\hfill feb-26-2007}
}
\vskip.4cm
Two-Loop Formfactors in Theories with Mass Gap\\
and $Z$-Boson Production
\vskip.4cm
}
\author{
A. Kotikov${}^{a,b}$, J.H. K\"uhn${}^c$ and O. Veretin${}^{b,d}$\\[0.5cm]
{\it ${}^a$Bogoliubov Laboratory of Theoretical Physics,}\\
{\it JINR, 141980 Dubna, Russia} \\
\\
{\it ${}^b$II Institute f\"ur Theoretische Physik,}\\
{\it Universit\"at Hamburg, 22761 Hamburg, Germany}\\
\\
{\it ${}^c$Institut f\"ur Theoretische Teilchenphysik,} \\
{\it Universit\"at Karlsruhe, 76128 Karlsruhe, Germany} \\
\\
{\it ${}^d$University of Petrozavodsk,} \\
{\it 185910 Petrozavodsk, Karelia, Russia} \\
%and\\
%{\it II Institute f\"ur Theoretische Physik,}\\
%{\it Universit\"at Hamburg, 22761 Hamburg, Germany}\\
}

\date{}
\maketitle



\abstract{ 
The non-factorizable two-loop corrections to
the formfactor both for a $U(1)\times U(1)$ 
and a $SU(2)\times U(1)$ gauge theory
with massive and massless gauge bosons respectively are evaluated at
arbitrary momentum transfer $q^2$.
%The two-loop formfactor both for a $U(1)\times U(1)$ 
%and a $SU(2)\times U(1)$ gauge theory
%with massive and massless gauge bosons respectively is evaluated at
%arbitrary momentum transfer $q^2$. 
The asymptotic behaviour for  $q^2\to\infty$ 
is compared to a recent calculation of Sudakov logarithms.  
The result is an important ingredient for the calculation
of radiative corrections to $Z$-boson production at hadron and lepton colliders. 
}

\vspace*{10mm}


%\setcounter{page}{0} \thispagestyle{empty}
\setcounter{equation}0

\section{Introduction}
Precise measurements of cross sections for the production of massive
and massless gauge bosons were one of the central topics of LEP
experiments. At the LHC similar reactions, namely the production of $W$
and $Z$-bosons, singly or in pairs, with or without additional quark
or gluon jets, will be crucial for precise studies of the electroweak
and strong interactions.
Single $W$- and $Z$-boson production will be used for the 
determination of parton
distributions and eventually even for luminosity
calibrations. A future linear collider, operating in the GIGA-$Z$ mode,
will measure the properties of the $Z$-resonance  with unprecedented
precision. All these measurements will rely on the theoretical knowledge
of radiative corrections to better than one percent accurracy, perhaps
even down to the level of several permille. 
QCD and electroweak radiative corrections, as well as their interplay,
thus will be crucial for the interpretation of these results.

QCD corrections to single $W$- and $Z$-production are identical to those
for the Drell--Yan process and have been evaluated in two-loop
approximation in \cite{Hamberg:1990np,Rijken:1994sh}, 
those for Higgs boson production
in \cite{Harlander:2003ai}. Electroweak corrections for the
on-shell process were computed some time ago (see e.g. \cite{Baur:2001ze}
and references therein).
The next step evidently requires to combine QCD and electroweak
effects, resulting in non-factorizable terms of order
$\alpha_{\rm weak}\alpha_s$. For the inclusive $Z$ decay rate these terms
have been calculated for final states with up-, down-, and bottom-quarks
\cite{Czarnecki:1996ei,Harlander:1997zb,Fleischer:1999iq}
and turned out to be relevant for the precise determination of the
strong coupling constant. However, these results cannot be directly
applied to the production process and to more differential
distributions. For the $Z$-boson such corrections for high $p_T$ distribution
have been obtained in \cite{Kuhn:2005az}.

In the present paper we describe conceptial developments and concrete
results which are important ingredients for the complete evaluation of
these non-factorizable terms of order $\alpha_{\rm weak}\alpha_s$.
In particular we consider those amplitudes which correspond to vertex
diagrams with a virtual gluon attached to one-loop electroweak
corrections. These are relevant for the ``mixed'' corrections of order
$\alpha_{\rm weak}\alpha_s$ to $Z$-boson production, and for hadronic $Z$
decay.
Essentially the same diagrams are also important ingredients for the
combination of photonic and weak corrections to $Z$ production in
electron-positron collisions and to leptonic $Z$ decays.

Our study identifies the infrared singular as well as the finite parts,
investigates the structure of these singularities and shows how they
can be combined with real radiation to arrive at a finite
result. The infrared finite remainder will be presented in analytical 
form in terms of generalized polylogarithms.

The form factor will also be investigated in the Sudakov limit $M^2/q^2\ll 1$. 
In the special case of an Abelian theory the result coincides
with the one of \cite{Penin} (see also \cite{Jantzen:2006jv}) 
and allows to contrast the logarithmic
approximation with the complete result. The calculational
method relies on an approach that has already been successfully employed 
in a number of cases \cite{FKV,Kniehl:2005bc}.
General considerations restrict the structure of the final
result to a sum of ``basis functions'' (in our case---generalized 
``harmonic'' polylogarithms up to fourth degree) 
with specific arguments and prefactors.
Calculating on one hand directly a large number of
terms in the low $q^2$ expansion with the technique of large mass
expansion, expanding the basis functions on the other hand, and
equating the results, the coefficients in front of the basis functions
can be determined. In a final step most of the basis functions are
transformed into Nielsen polylogarithms, leading to a fairly compact
result whose asymptotic behaviour can be analyzed in a straightforward
manner.

To facilitate the discussion, we present, in a first step, in section 2,
the results for a $U(1)\times U(1)$ theory with one of the gauge bosons
taken to be massive, the other one massless. The explicit analytical result
confirms the factorization of the infrared singularities and allows to
identify the infrared-finite remainder. 
%In the limit of large $q^2$
%(Sudakov limit) the result agrees with the one obtained in \cite{Penin}, and
%the relative size of subleading, power suppressed terms can be studied.
In section 3 the formalism will be extended to a massive nonabelian theory
and applied to the complete set of virtual corrections of order
$\alpha_{\rm weak}\alpha_s$, contributing to $Z$-boson production and decay.
The triple-boson coupling leads to additional diagrams with additional
generalized polylogarithms,  which cannot easily be transformed into
Nielson's polylogarithms. However, they can be evaluated
numerically with high precision \cite{hplog} and their asymptotic behaviour is
under control. The paper concludes with a brief summary. Much of the
formulae and calculational details will be collected in the Appendices.

\section{Abelian Theory}
For definiteness and simplicity we will, in a first step, consider the
form factor in a ficticious $U(1)\times U(1)$ theory with one massive and
one massless gauge boson and with coupling constants $\alpha$ and
$\alpha'$ respectively.

For the Abelian theory the form factor $F$ will be defined as matrix 
element of an external current
\begin{equation}
\gamma_\mu F(q,M) = \int dx e^{-ixq} \langle \psi'|J_\mu(x)|\psi \rangle.
\end{equation}
Here $\psi$ and $\psi'$ denote on-shell massless fermions of momenta $p$ and
$p'=p+q$, respectively,  $M$ the mass of the gauge  boson. 

In a perturbative expansion
\be
F(\alpha,\alpha',q,M,\varepsilon)=\sum_{m,n=0}^\infty
\left(\frac{\alpha}{4\pi}\right)^m \left(\frac{\alpha'}{4\pi}\right)^n
f^{(m,n)}(q,M,\varepsilon)
\ee
%
one needs to evaluate the expansion coefficients $f^{(m,n)}$. In Born and
one-loop approximation they are given by
%
\bea
f^{(0,0)}&=& 1 \,, \\
f^{(1,0)}&=&
   - \frac{7}{2} - \frac{2}{z}
   + \frac{2+3z}{z} \log(-z)
   + \frac{2(1+z)^2}{z^2} \Big( {\rm Li}_2(1+z) - \frac{\pi^2}{6}
   \Big)\,,\label{1loopMass}
\\
f^{(0,1)}&=& \left(\frac{\bar\mu^2}{-q^2}\right)^\varepsilon \left(
   -\frac{2}{\varepsilon^2} - \frac{3}{\varepsilon} 
   - 8 + \zeta_2 + \varepsilon \Big(
   -16 + \frac{3}{2} \zeta_2
   + \frac{14}{3} \zeta_3 \Big) \right) \,,
\label{1loopQCD}
\eea
where $z=q^2/M^2+i0$, 
%$\bar{\mu}$ is the renormalization scale in $\bar{MS}$ scheme 
$\zeta_n=\zeta(n)$ is the Riemann $\zeta$-function.
The infrared singularities are controlled by dimensional regularization 
in $d=4-2\epsilon$ dimension and $\bar\mu$ is related to the scale
parameter of dimensional regularisation by $\bar\mu^2 = \mu^2 4\pi
\exp(-\gamma)$.
%and the infrared singularities are controlled by
%dimensional regularization in $d=4-2\varepsilon$ dimensions.
In the euclidean region $q^2<0$, so that no
imaginary parts appear in the above formulae.



  The {\em two}-loop result for the massless case, $f^{(0,2)}$, can be
found e.g. in \cite{Gonsalves:1983nq,vanNeerven:1985xr}. 
The two-loop result for the fully massive case,
$f^{(2,0)}$, is only known in the large $q^2$ limit 
%\cite{Kuhn:1999nn}. 
\cite{Penin,Kuhn:1999nn} and \cite{Jantzen:2005xi} in abelian and non-abelian
theories respectively. The
evaluation of the mixed corrections is drastically simplified by the
fact that the infrared singularities factorize within infrared evolution
equation approach \cite{Fadin:1999bq,Kuhn:1999nn}, which gives in our case
\be
F(\alpha,\alpha',q^2/M^2,\varepsilon) = F_{\rm massless}(\alpha',q,\varepsilon)
                                      \tilde F(\alpha,\alpha',q^2/M^2)\,,
\label{infrared}
\ee
with $F_{\rm massless}=\sum (\alpha'/4\pi)^n f^{(0,n)}(q,\varepsilon)$ denoting the
formfactor for the massless theory and $\tilde F$ being free from infrared
singularities. 
The function $\tilde F$ can again be expressed as double series, and the
coefficients depend on the ratio $z=q^2/M^2$ only.
The terms $\tilde F^{(m,0)}=f^{(m,0)}$ coincide by definition with those
valid for the massive $U(1)$-theory. 
The evaluation of the nonfactorizable part of the two-loop contribution
\be
\tilde F^{(1,1)}\equiv \phi(q^2/M^2)
\ee
will be the central result of this section.
% We stress that $\phi(z)$ is
%the non-factorizable part of the two-loop correction to the formfactor of 
%order $\cal{O}(\alpha\alpha')$.

The Feynman diagrams
necessary for this computation have two thresholds: at $q^2=0$ and at $q^2=M^2$.
The analytical structure of vertex diagrams of this type has been explored in \cite{FKV}.
The coefficients of an expansion in $q^2/M^2$ (and $M^2/q^2$) can
always be expressed as combinations of so-called harmonic sums \cite{harmonicsums}
or more generally --- nested harmonic sums \cite{harmonicsums1}. These
sums correspond to (generalized) polylogarithms (\cite{Devoto}) \cite{Lewin}
of arguments $\pm q^2/M^2$ and their generalizations --- harmonic 
polylogarithms \cite{HPLdefinition} (see also \cite{MultiplePL}).
This structure suggests the following method for the
evaluation of Feynman integrals. First, using the method of
large mass expansion \cite{asymptotic}, one calculates a large number of
coefficients of the series in $q^2/M^2$. 
From the basis functions (polylogarithms) one then constructs an Ansatz (see Appendix A)
with unknown coefficients $r_i$.
Equating Ansatz and series one obtains a unique answer for parameters $r_i$.
This method has been applied earlier \cite{FKV98,FKV} to various scalar vertex 
masterintegrals.
(In a different context the method has also been applied in \cite{Harlander:2005rq}). 
Here it is applied to amplitudes deduced
from a a set of realistic Feynman diagrams representing a physical
process and leading to amplitudes with irreducible numerators and
shrunken lines.

A few comments on this procedure are in order. First, the main
problem is to write down the correct prefactors in the Ansatz.
Empirically one finds that the presence of a numerator or the absence of
a line may lead to the additional factors $M^2/q^2$ or $(M^2/q^2)^2$
in front of polylogarithms%
\footnote{
In a series representation such multiplications lead to shifts 
of the summation index in $c_n$. Indeed, if $z=q^2/M^2$ then, e.g.
$
  \frac{1}{z}\sum\nolimits_{n=1}c_n z^n = c_1 + \sum\nolimits_{n=1}c_{n+1}z^n
$
and so on.
}. % END of footnote
Therefore such factors should also be included in the Ansatz.
Second, only five functions could not be represented as Nielsen polylogarithms
with the argument $q^2/M^2$. These remaining functions belong to the class of harmonic
polylogarithms \cite{HPLdefinition} discussed in more detail in the
Appendix B.

%  For the completeness we the details of the Ansatz in the Appendix B.

Instead of expanding the amplitude in $q^2/M^2$
one could find the differential equation (see \cite{DEM}) for a diagram and 
again apply an Ansatz based on polylogarithms. This approach has  recently been used
for similar two-loop vertex diagrams in \cite{Bonciani1}.

  Altogether 16 one-particle-irreducible two-loop vertex diagrams
contribute to the formfactor. These diagrams can be obtained from the one-loop
one shown in Fig. \ref{ffdias}a by adding one gluon line.
The two-loop, one-particle reducible diagrams
which are obviously products of one-loop diagrams contribute to the term 
$\alpha\alpha' f^{(0,1)} f^{(1,0)}$ and are not repeated here. We also do not
display the contributions to the fermionic wave function renormalization, which receives
contributions from additional 6 diagrams.
For the generation of the input the program DIANA \cite{DIANA} has been used,
for the evaluation and expansion a program written in FORM \cite{FORM}. 
The evaluation of the Dirac traces has lead to about 700 different
integrals. For most of them the asymptotic expansion was performed up to 
order 45 which required in total several hours of CPU time on a Pentium IV
processor. 
For the remaining, most complicated cases (nonplanar diagram)
up to 60 expansion coefficients had to be computed. For this purpose
the parallel version of FORM \cite{Tentyukov:2004hz}, running on an
SGI machine with multiprocessor SMP architecture, was used.

  The function $\phi(z)$ can be cast into the
following form (here and below $z=q^2/M^2 + i0$)
\begin{eqnarray}
\label{res}
\phi(z) \z=\z 
      \frac{(1+z)^2}{z^2}  \Biggl(
          ( 6 \LOG^2 + 24 \zeta_2  - 24 \zeta_3 ) \log(1 + z)
         +( - 4 \LOG^2 - 6 \LOG - 20 \zeta_2 ) \log^2(1 + z)
\nonumber\\
\z\z
          + \frac{8}{3} \log^3(1 + z) \LOG
          + 8 \log^2(1 + z) \Li_2( - z)
          - 12 \log(1 + z) \Li_2( - z)
\nonumber\\
\z\z
          - 16 \log(1 + z) \Li_3( - z)
          + 16 \log(1 + z) \Si_{1,2}( - z)
          - 16 \Li_2( - z) \zeta_2
          - 4 \Li_2( - z) \LOG^2
\nonumber\\
\z\z
          - 8 \Li_2^2( - z)
          + 16 \Li_3( - z) \LOG
          - 24 \Li_4( - z)
          - 12 \Si_{1,2}( - z)
          + 16 \Si_{1,2}( - z) \LOG
          + 16 \Si_{1,3}( - z)
\nonumber\\
\z\z
          - 16 \Si_{2,2}( - z)
          + 24 h(z)
          - 48 H_3(z)
          + 8 H_2(z)
          + 32 H_4(z)
          \Biggr)
%
\nonumber\\
\z\z
       + \frac{1+3z+z^2}{z^2} \Biggl(
          - 32 \Li_2(z) \zeta_2
          - 8 \Li_2(z) \LOG^2
          - 8 \Li_2(z) \Li_2(z)
          + 32 \Li_3(z) \LOG
\nonumber\\
\z\z
          - 48 \Li_4(z)
          + 32 \Si_{2,2}(z)
          \Biggr)
%
       + \frac{1-z^2}{z^2} \Biggl(
            72 \log(1 - z) \zeta_2
          + 18 \log(1 - z) \LOG^2
\nonumber\\
\z\z
          + 36 \log(1 - z) \Li_2(z)
          + 36 \Li_2(z) \LOG
          + 72 \Si_{1,2}(z)
         \Biggr)
%
       + \frac{2+3z}{z}  ( 32 \zeta_2 + 12 \zeta_3 )
\nonumber\\
\z\z
       - \frac{34 + 51z}{z} \LOG
       + \frac{16 + 23z}{z} \LOG^2
       - \frac{2(1-z)(13 + 27z)}{z^2} \log(1-z)
       + \frac{4(3+4z)}{z^2} \Li_2( - z) \LOG
\nonumber\\
\z\z
       - \frac{2(1+z)(3+5z)}{z^2} \Bigl( \log(1 + z) \LOG + \Li_2( -z ) \Bigr)
       + \frac{4(11 + 9z)}{z}\Li_2(z)
\nonumber\\
\z\z
       - \frac{4(3 + 2z - 3z^2)}{z^2} \Li_3( - z)
       - \frac{4(9 + 4z - 6z^2)}{z^2} \Li_3(z)
       + \frac{2(8-z)}{z},
\end{eqnarray}
%
where $L=\log(-q^2/M^2)$, $\zeta_a=\zeta(a)$ is the Riemann $\zeta$-function,
$S_{a,b}(z)$ are Nielsen polylogarithms \cite{Devoto}.
The functions $h(z)$ and $H_j(z)$ are defined and discussed in Appendix B.


Note that the result (\ref{res}) depends only on one Vermaseren--Remiddi
function (we choose $H_2(z) \equiv H_{-1,0,0,1}(z)$) 
and a combination of Nielsen
polylogarithms (see Appendix A). This observation should be correct also for 
any results which can be representaed as a combination of massless Feynman
integrals and/or diagrams with one-massive particle thresholds.

%
\section{$Z$-Production}
%

For definitenes and simplicity, consider, in the next step,
$Z$-boson production in quark-antiquark annihilation.
To fix the notation, we recapitulate the one-loop results.
The weak corrections to the Born term
can be split into those involving the exchange of $W$- and $Z$-bosons, 
(Fig.\ref{ffdias}(a)) and those involving the
triple-boson coupling (Fig.\ref{ffdias}(b)).
The combination of photonic and QCD corrections follows
essentially from the two-loop  QED or QCD results and will not be addressed here.

  For a light quark the form factor can be decomposed as follows
\begin{eqnarray}
F(q^2)_\mu = \gamma_\mu\frac{1+\gamma_5}{2} F_R(q^2) 
        + \gamma_\mu\frac{1-\gamma_5}{2} F_L(q^2) \,.
\end{eqnarray}
At the Born level the expressions for the form factors $F_R$
and $F_L$ are given by 
\begin{eqnarray}
  F_R = i \frac{e}{s_W} g_R \,,\,\,\,
  F_L = i \frac{e}{s_W} g_L \,,
\end{eqnarray}
with $g_R=-Qs_W^2/c_W$ and $g_L=(I_3-Qs_W^2)/c_W$
being the right- and left- handed couplings of a quark to the $Z$-boson.
Here $I_3$ is the third component of the isospin of a quark,
$Q$ its electric charge 
and $s_W=\sin\theta_W$ and $c_W=\cos\theta_W$ denote sine and cosine
of the weak mixing angle, respectively.

Including radiative corrections and adopting a form similar 
to eq.\ref{infrared} the formfactors can be cast into the following form
\begin{eqnarray}
F_R &=&  i\frac{e}{s_W} 
    \Big( 1 + C_F\frac{\alpha_s}{4\pi} f^{(0,1)}\Big)
    \Bigg[ g_R + \frac{\alpha}{4\pi s_W^2} g_R^3 \rho_{\rm A}(q^2/m_Z^2,m_Z^2) 
  + C_F\frac{\alpha_s}{4\pi}\frac{\alpha}{4\pi s_W^2}
      g_R^3 \phi_{\rm A}(q^2/m_Z^2) \Bigg] \,,
\nonumber\\
%
F_L &=&  i\frac{e}{s_W} 
      \Big( 1 + C_F\frac{\alpha_s}{4\pi} f^{(0,1)}\Big)
      \Bigg[ g_L
\nonumber\\ 
         & +& \frac{\alpha}{4\pi s_W^2} \Bigg(
           g_L^3 \rho_{\rm A}(q^2/m_Z^2,m_Z^2) 
            + \frac{g_L}{2} \rho_{\rm A}(q^2/m_W^2,m_W^2) 
            + c_W \frac{I_3}{2} \rho_{\rm NA}(q^2/m_W^2,m_W^2) \Bigg)  
\nonumber\\
  &+& C_F\frac{\alpha_s}{4\pi}\frac{\alpha}{4\pi s_W^2}
      \Bigg( g_L^3 \phi_{A}(q^2/m_Z^2)
            + \frac{g_L}{2} \phi_{\rm A}(q^2/m_W^2) 
            + c_W\frac{I_3}{2} \phi_{\rm NA}(q^2/m_W^2)  \Bigg)\Bigg]  \,,
\end{eqnarray}
where first factors in the brackets in the above equations represent the
QCD corrections.
Terms given by $\rho_{\rm A}$ and $\rho_{\rm NA}$ account for
the one-loop electroweak corrections. The abelian part $\rho_{\rm A}$
is defined by the diagram of the abelian type (Fig. 1a)
and obviously closely related to $f^{(1,0)}$ defined in eq.\ref{1loopMass}. 
The unrenormalized result% 
 \footnote{We shall not discuss issues
 related to renormalization, since the non-factorizable part, which is the
 quantity of interest in this paper, is independent of the
 renormalization scheme.}
is given by
\begin{equation}
  \rho_{\rm A}(z,M^2) =
   \frac{1}{\varepsilon}
   - \ln(M^2/\bar\mu^2)
   - 4 - \frac{2}{z}
   + \frac{2+3z}{z} \log(-z)
   + \frac{2(1+z)^2}{z^2} \Big( {\rm Li}_2(1+z) - \frac{\pi^2}{6} \Big)\,.
\end{equation}
The nonabelian part $\rho_{\rm NA}(z,M_2)$ receives corrections from
both diagrams of Fig.1a and Fig.1b. It is given by
\begin{eqnarray}
  \rho_{\rm NA}(z,M^2) &=& - 2 \rho_{\rm A}(z,M^2) - 2\Lambda(z,M^2) \,, \\
  \Lambda(z,M^2) &=& - \frac{3}{\varepsilon}
    + 3 \ln(M^2/\bar\mu^2) - 2 + \frac{2}{z} 
    - \Big(1 + \frac{2}{z}\Big) \sqrt{1 - \frac{4}{z} } \, l
    - \Big( 1 + \frac{1}{2z}\Big) \frac{4}{z} l^2 
\end{eqnarray}
with
\begin{equation}
   l = \ln\left( \frac{\sqrt{1-4/z}+1+i0}{\sqrt{1-4/z}-1+i0} \right)\,.
\end{equation}
The function $\Lambda(z)$ can be taken from \cite{Grzadkowski:1986pm}
(see also \cite{Bohm:1986rj,Bardin:1999ak,Berends:1987ab} and references
therein for one-loop calculations in the Standard Model).
We do not include the terms from the renormalization of the coupling
and the $Z$-boson wave function%
\footnote{Hence the function $\Lambda(z)$ considered 
   in \cite{Grzadkowski:1986pm}
   differs by subtracting the term $3/\varepsilon-3\ln(M^2/\bar\mu^2)-1/2$.
   Furthermore, a typo in \cite{Grzadkowski:1986pm} has been corrected,
   flipping the signs of the terms proportional to $l$ and $l^2$}
which follow from textbook prescriptions and will not be 
considered in this work.

Evaluated for arbitrary $q^2\neq M_Z^2$, the above results are gauge dependent 
and are presented in Feynman gauge. 
For the offshell case they can be considered as
building blocks for a complete calculation.

The functions $\phi_{\rm A}(z)$ and $\phi_{\rm NA}(z)$, representing
the non-factorizable terms of ${\cal{O}}(\alpha\alpha_{s})$, are written in
a form completely analogous to the electroweak one-loop terms.  
The function $\phi_{\rm A}(z)$ has been given in the previous section.
The nonabelian part $\phi_{\rm NA}(z)$ involves new functions --- generalized 
polylogarithms. 
%
%
% PICTURE DIAGRAMS
\begin{figure}[h]
\centerline{
%\raisebox{30mm}{$\Delta r^{(2)}_{\rm bos}$}\hspace{-8mm}
%\raisebox{58mm}{$\times10^{-5}$}\hspace{-4mm}
%\begin{center}
\hbox{\hspace*{-0mm}\centerline{\epsfig{file=ffdias.eps, angle=0, width=80mm}}}
}
\caption{\label{ffdias}
  Diagrams, contributing to the vertex $Zq\bar{q}$ (a) and (b). 
The two-loop diagrams are obtained by attaching one virtual gluon in all possible ways.
%The case (b) represents nonabelian part. That gives contribution 
%$\phi_{\rm NA}(z)$ in the text. Diagram (a) with $W$ exchange also contributes
%to $\phi_{\rm NA}(z)$.
The case (b) represents the nonabelian part and contributes to 
$\phi_{NA}(z)$. Diagram (a) with $W$ exchange contributes 
both to $\phi_{A}(z)$ and $\phi_{NA}(z)$. The cross indicates the coupling 
to the external $Z$ boson
}
\end{figure}
Our result in Feynman gauge reads
\begin{eqnarray}
\phi_{\rm A}(z) &=& \phi(z) \,, \\ 
\label{resf2}
\phi_{\rm NA}(z) &=& - 2 \phi_{\rm A}(z)  \nonumber\\
     \z\z+
         4\frac{8-5z}{z}
       + 48\frac{1+2z}{z^2} H_{0,-r,-r,-1}(-z)
        - 12 H_{-r,-r,-1}(-z)
       + 8 H_{0,-r,-r}(-z)
\nonumber\\
\z\z
       + 4\frac{6+2z-3z^2}{z^2} H_{-r,-r}(-z)
       - 6\frac{(4-z)(4+3z)}{z^2} gH_{-r,-1}(-z)
\nonumber\\
\z\z
       - 2\frac{(4-z)(6+7z)}{z^2} gH_{-r}(-z)
       + \frac{16+23z}{z} \Bigl(
            8 \zeta_2
          + 2 \LOG^2
          \Bigr)
       - 4\frac{12-11z^2}{z^2} \Li_3(z)
\nonumber\\
\z\z
       + 2\frac{66+49z}{z} \Li_2(z)
       - 4\frac{(1-z)(13+34z)}{z^2} \log(1-z)
       - 16\frac{5+9z}{z} \LOG
\nonumber\\
\z\z
       + \frac{(1-z^2)}{z^2} \Biggl(
            96 \log(1 - z) \zeta_2
          + 24 \log(1 - z) \LOG^2
\nonumber\\
\z\z
\qquad\qquad\qquad\qquad\qquad
          + 48 \log(1 - z) \Li_2(z) 
          + 48 \Li_2(z) \LOG
          + 96 \Si_{1,2}(z)
          \Biggr)
\nonumber\\
\z\z
       + \frac{(1+4z+z^2)}{z^2} \Biggl(
          - 32 \Li_2(z) \zeta_2
          - 8 \Li_2(z) \LOG^2
          - 8 \Li_2^2(z)
\nonumber\\
\z\z
\qquad\qquad\qquad\qquad\qquad
          + 32 \Li_3(z) \LOG
          - 48 \Li_4(z)
          + 32 \Si_{2,2}(z)
          \Biggr)\,.
\end{eqnarray}
The function $\phi_{\rm NA}(z)$ receives contributions not only from digrams
of Fig. 1(b) but also from those of Fig. 1(a) with the exchange of $W$-boson.
The functions $H_{\dots}$ are considered in more detail in Appendix C.
For the special case $q^2=M^2$ one finds
%
\begin{eqnarray}
\phi_{\rm A}(1+i0) \z=\z
        14
          + 72 \zeta_2 l_2
          - 64 \zeta_2 l_2^2
          - \frac{16}{3} l_2^4
          + 22 \zeta_2
          - 28 \zeta_3
          + 16 \zeta_4
          - 128 {\rm Li}_4( {\textstyle \frac12} ) \nonumber\\
  \z\z     + i\pi   (
            85
          + 32 l_2
          + 24 l_2^2
          - \frac{32}{3} l_2^3
          + 14 \zeta_2
          - 120 \zeta_3
          )  \\
\z=\z -2.1073 - 19.0331 i\,,\\
\phi_{\rm NA}(1+i0) \z=\z
       - 16
          - 144 \zeta_2 l_2
          + 128 \zeta_2 l_2^2
          + \frac{32}{3} l_2^4
          + \frac{70}{3} \zeta_2
          + \frac{184}{3} \zeta_3
          - 236 \zeta_4  \nonumber\\
  \z\z
          + 26 \frac{\pi}{\sqrt{3}}
          + 256 {\rm Li}_4( {\textstyle \frac12} )
          - 84 \frac{1}{\sqrt{3}}{\rm Ls}_2( {\textstyle \frac{\pi}{3}} )
          - \frac{16}{3} \pi {\rm Ls}_2( {\textstyle \frac{\pi}{3}} ) 
          + 96 \Bigl( {\rm Ls}_2( {\textstyle \frac{\pi}{3}} ) \Bigr)^2  \nonumber\\
  \z\z    + i\pi   (
            54
          - 64 l_2
          - 48 l_2^2
          + \frac{64}{3} l_2^3
          - 28 \zeta_2
          + 48 \zeta_3
          ) \\
 \z=\z -7.5880 + 16.7194 i\,,
\end{eqnarray}
%
with $l_2=\log2$. Substituting the actual masses of the $W$- and
$Z$-bosons ($z=m_Z^2/m_W^2=1.2856$) we find:
\begin{eqnarray}
\phi_{\rm A}(1.2856+i0) \z=\z -1.3598 - 30.4095 i\,,\\
\phi_{\rm NA}(1.2856+i0)     \z=\z -10.1248 + 35.0336 i\,.
\end{eqnarray}

  In the limit $q^2\to\infty$ the function $\phi$, as given by Eq. (\ref{res})
coincides with
the result of \cite{Penin} where the power unsupressed logarithmic
and constant part have been evaluated. For the leading and the first
power suppressed  term we find 
\bea
\phi(z)& =&
(3 - 24\zeta_2 + 48\zeta_3)\log(-z)  - 2 + 40\zeta_2 - 84\zeta_3 + 14\zeta_4 \nonumber\\
  &+& \frac{1}{z} \bigg(
      ( - 26 + 8\zeta_2 ) \log^2(-z) 
      + ( - 120 - 16\zeta_2 + 128\zeta_3 ) \log(-z) \nonumber\\
  && \quad        - 188 - 8\zeta_2 - 8\zeta_3 + 116\zeta_4  \bigg) + \Oz{2}
\eea
In Fig. \ref{asympphi}  the exact result is compared with the Sudakov approximation
and with the approximation including the first power-suppressed term. 
For electroweak interactions the mass of the gauge boson can be taken to
be of order 100 GeV, the characteristic energy of order one to two TeV.
For  one TeV the relative error of the Sudakov approximation (the
logarithmic plus constant term) amounts to 15\%, at 2TeV it is reduced
to 2.5\%. 


  Assuming that this behaviour is characteristic for the general case, 
power suppressed terms start to become important for $M^2/s$ below 1/100.
Furthermore inclusion of the first subleading term only does not lead to a 
significant improvement of the approximation.



%
%
% PICTURE
\begin{figure}[h]
\centerline{
%\raisebox{30mm}{$\Delta r^{(2)}_{\rm bos}$}\hspace{-8mm}
%\raisebox{58mm}{$\times10^{-5}$}\hspace{-4mm}
%\begin{center}
\hbox{\hspace*{-0mm}\centerline{\epsfig{file=p1.ps, angle=-90, width=150mm}}}
}
%\hspace{95mm}\raisebox{-0mm}{\Large $z$}
%\hbox{\raisebox{-105mm}{\hspace*{-80mm}$q^2/M^2$}}\\
%\hspace{-190mm} \raisebox{-50mm}{$f$}
\caption{\label{asympphi}
  Non-factorizable two-loop correction to the abelian formfactor in the
euclidean regime ($z=q^2/m^2$).
The solid line represents the exact result, the dashed line the Sudakov 
approximation and the dash-dotted line includes the power suppressed terms.
}
\end{figure}




\section{Conclusions}
Using the technique of asymptotic expansions and the knowledge
of the general structure of integrals we evaluated analytically the 
two-loop formfactor in a $U(1)\times U(1)$ 
theory with one massive and one massless gauge boson
and constant terms have been evaluated. The Sudakov approximation gives 
a good description for $M^2/s$ below 1/100. For lower energies mass 
suppressed terms  must be included.

In the Sudakov limit full agreement is obtained with \cite{Penin},
where the logarithmic and constant terms had been evaluated obtained.
We furthermore perform the same caldulation for a 
 $SU(2)\times U(1)$ theory and derive the non-factorizable part of the
two-loop formfactor in the Standard Model. 
As an application we evaluate the mixed virtual $O(\alpha\alpha_s)$
radiative correction for Drell--Yan production of the $Z$-boson.


%Using the technique of asymptotic expansions and the knowledge
%of the general structure of integrals we evaluated analytically the 
%two-loop formfactor in a $U(1)\times U(1)$ 
%theory with one massive and one massless gauge boson.
%In the Sudakov limit full agreement is abtained with \cite{Penin},
%where the logarithmic and constant terms had been evaluated obtained.
%We furthermore perform the same caldulation for a 
% $SU(2)\times U(1)$ theory and derive the non-factorizable part of the
%two-loop formfactor in the Standard Model. 
%As an application we evaluate the mixed virtual $O(\alpha\alpha_s)$
%radiative correction for Drell--Yan production of the $Z$-boson.\\

{\it Acknowlegments.} We thank M.~Kalmykov for useful comments, discussions
%{\bf and information about reference \cite{Czakon}}
and communication 
and M.~Tentyukov for his help with DIANA. We acknowledge T.~Gehrmann
for information about the numerical program \verb+hplog+
and B.~Jantzen for the carefull reading the paper.
This work was supported by BMBF under grants
No. 05HT6VKA, 05HT4GUA4 and HGF grant No. NG-VH-008.
%At the beginning of the study A.K. was supported in part by an Alexander von
%Humboldt fellowship.
A.K. is supported in part by an Alexander von Humboldt Foundation
(a renewed academic stay in Germany) and RFBR Foundation through grant No. 05-02-17645-a.



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

\section{Appendix A}

\label{App:A}

\def\theequation{A\arabic{equation}}
\setcounter{equation}0

  In this Appendix we give the details of the Ansatz which is used to obtain
the analytical result for the amplitudes.

  Generally the result of the finite part of expansion for a given diagram $J$ looks like
\begin{equation}
  J = M^g \sum_{n=1}^{n_{\rm max}} z^n \Big( a_n + b_n \ln z + c_n \ln^2 z \Big) \,,
\label{dia}
\end{equation}
where $z=q^2/M^2$ and $g$ is the mass dimension of the diagram. In (\ref{dia}) we substitute
a linear combination
\begin{equation}
  r_1 \chi_1(z) + r_2 \chi_2(z) + \dots + r_N \chi_N(z)
\end{equation}
where $r_j$ are unknown rational numbers and $\chi_j(z)$ are basis functions.

  In our case the basis functions are chosen as follows. In the abelian case
it is enough to include all harmonic polylogarithms of the weigth one to four.
This includes
$$
  \ln(1\pm z)\,, {\rm Li}_2(\pm z)\,, {\rm Li}_3(\pm z)\,, {\rm Li}_4(\pm z)\,, 
     S_{1,3}(\pm z)\,, S_{2,2}(\pm z) \,,
$$
with $S_{a,b}$ being Nielsen polylogarithm. In addition one should add the following
harmonic polylogarithms
$$
  H_{-1,0,1}(-z)\,,  H_{-1,0,1,1}(-z)\,, H_{-1,0,0,1}(-z) \,, H_{0,-1,0,1}(-z) \,,
$$
which are considered in the Appendix B.
  In the nonabelian part (see Fig.1(b)) new functions appear which are listed
in formulae (C1)-(C6).

  We should also note that the above functions could be multiplied by additional
factors $1/z$ or $1/z^2$. All these combinations should be included in the Ansatz.
We also observed that functions proportional to $1/z^3$ and higher inverse powers of $z$
do not appear.

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

\section{Appendix B}

\label{App:B}

\def\theequation{B\arabic{equation}}
\setcounter{equation}0


%  In this Appendix 
Here we consider the asymptotic behaviour of the most
complicated basis functions in the limit $z=q^2/m^2\rightarrow\infty$. 
Most of basis functions can be expressed in terms
of Nielsen polylogarithms and then the standard transformations formulae
can be applied to go from argument $z$ to $1/z$ 
(see \cite{Devoto,Lewin}). Therefore we will consider here only
the five special cases, mentioned previously,
where complications arise.

  In our calculation the following functions appear in addition 
to usual Nielsen polylogarithms:
%
\begin{eqnarray}
h(z) \z=\z H_{-1,0,1}(z) \,,  \nonumber\\
H_1(z) \z=\z H_{-1,0,1,1}(z) \,, \nonumber\\
H_2(z) \z=\z H_{-1,0,0,1}(z) \,, \nonumber\\
H_3(z) \z=\z H_{-1,-1,0,1}(z) \,, \nonumber\\
H_4(z) \z=\z H_{0,-1,0,1}(z) \,,  \nonumber
\end{eqnarray}
%
where $H_{a,\dots,d}(z)$ are harmonic polylogarithms defined in \cite{HPLdefinition}.

  These functions correspond to the alternating Taylor series in $z$:
\begin{eqnarray}
h(z) \z=\z - \sum\limits_{n=1}^{\infty} \frac{S_{-2}(n-1)}{n}(-z)^n\,, \nonumber\\
H_1(z) \z=\z - \sum\limits_{n=1}^{\infty} \frac{S_{-2,1}(n-1)}{n}(-z)^n\,, \nonumber\\
H_2(z) \z=\z - \sum\limits_{n=1}^{\infty} \frac{S_{-3}(n-1)}{n}(-z)^n\,, \nonumber\\
H_3(z) \z=\z - \sum\limits_{n=1}^{\infty} 
     \frac{S_{-3}(n-1)+S_{-2,1}(n-1)-S_1(n-1)S_{-2}(n-1)}{n}(-z)^n\,, \nonumber\\
H_4(z) \z=\z - \sum\limits_{n=1}^{\infty} \frac{S_{-2}(n-1)}{n^2}(-z)^n\,, \nonumber
\end{eqnarray}
with finite harmonic sums $S_a(n)=\sum_{j=1}^n1/j^a$ 
and $S_{-a}(n)=\sum_{j=1}^n(-1)^j/j^a$ and
$S_{-2,1}(n) = \sum_{j=1}^n (-1)^j S_1(j)/j^2$.
It is interesting to note that the function $H_1$ cancels in the final 
result (\ref{res})
for the formfactor but is present in the particular integrals.

  Following \cite{FKV} it is not difficult to write down simple integral representations 
for the above series, e.g.
\begin{eqnarray}
\label{inth}
h(z) \z=\z \int\limits^{z}_{0} \frac{dx}{1+x}\,\Li_2(x)\,, \\
\label{intH1}
H_1(z) \z=\z \int\limits^{z}_{0} \frac{dx}{1+x}\,\Si_{1,2}(x)\,, \\
\label{intH2}
H_2(z) \z=\z \int\limits^{z}_{0} \frac{dx}{1+x}\,\Li_{3}(x)\,, \\
\label{intH3}
H_3(z) \z=\z  \log(1+z)h(z) - \int\limits^{z}_{0} \frac{dx}{1+x}\,\Li_{2}(x)\log(1+x)\,, \\
\label{intH4}
H_4(z) \z=\z \log(z)h(z) - \int\limits^{z}_{0} \frac{dx}{1+x}\,\Li_{2}(x)\log(x)\,. 
\end{eqnarray}
%
%Here $S_{1,2}(z)$ is Nielsen polylogarithm \cite{Devoto}.

Now the integrals can be expressed in terms of Nielsen polylogarithms 
of {\it nonlinear} arguments and only one harmonic polylogarithm
function $H_2$ (this choice being not unique, however). We have
\footnote{ Eq. (\ref{r1}) is taken from \cite{FKV}. Equations similar
to (\ref{r3})-(\ref{r2}) have been considered also in  
\cite{DavyKalmy,Czakon}. In particular, Eq. (\ref{r3}) coincides with
(A.29) of \cite{DavyKalmy}.}
%{\bf sled. formuly provereny chislenno}
%
\begin{eqnarray}
h(z) \z=\z \frac{1}{2} \Si_{1,2}(z^2) - \Si_{1,2}(z)-\Si_{1,2}(-z)
+\ln(1+z)\Li_{2}(z), \label{r1} \\
H_1(z) \z=\z \log(1+z) \Si_{1,2}(z) + \frac{1}{4} \Si_{1,3}(z^2)
             - \Si_{1,3}(-z) + \frac{1}{2} \Phi(z), \label{r3} \\
H_3(z) \z=\z \log(1+z) \Bigl( \frac{1}{2} \Si_{1,2}(z^2) - \Si_{1,2}(z)-\Si_{1,2}(-z) 
                    \Bigr) \nonumber \\
     \z\z  + \frac{1}{2}\log^2(1+z)\Li_{2}(z)
     + \frac{1}{4} \Si_{1,3}(z^2) - \Si_{1,3}(z) - \frac{1}{2} \Phi(z), \label{r4}\\
H_4(z) \z=\z \frac{1}{4} \Si_{2,2}(z^2) - \Si_{2,2}(z)-\Si_{2,2}(-z)
           + \log(1+z) \Li_3(z) - H_2(z), \label{r2} 
\end{eqnarray}
%
where
%
\begin{eqnarray}
\Phi(z) \z=\z - \frac{15}{8}\zeta_4 + \frac{1}{6} \log^3{s} \log{z} 
           + \frac{1}{2}\log^2{s} \Bigl(\Li_{2}(s) - \Li_{2}(-s)\Bigr) \nonumber \\
          \z\z -\log{s} \Bigl(\Li_{3}(s) - \Li_{3}(-s)\Bigr) +\Li_{4}(s) - \Li_{4}(-s),
\label{r5} \nonumber
\end{eqnarray}
%
with $s = (1-z)/(1+z)$.

  In order to find the asymptotic behaviour for $z\to\infty$ one needs to 
use the standard
formulae for polylogarithms and for the function $H_2(z)\equiv H_{-1,0,0,1}(z)$
the inversion formula (A.6) from \cite{DavyKalmy}. 
It is important to take care of imaginary parts, therefore
we approach the cut in $q^2$-plane from above,
which means that $z$ is replaced by $z+i0$.  Thus we obtain
% these are corrected formula after check with ``hplog''
\begin{eqnarray}
h(z+i0) &=& - \frac{1}{6}\log^3{z} + 2\zeta_2 \log{z}- \frac{3}{2}\zeta_3
      + \frac{1}{z} \bigg( - \frac{1}{2} \log^2{z} - \log{z} + 2\zeta_2 \bigg) \nonumber\\
      &+& i \pi \left\{ \frac{1}{2}\log^2{z} - \frac{1}{2}\zeta_2 
           + \frac{1}{z} \log{z} + \frac{1}{z} \right\}
      + {\cal O} \Bigl(\frac{1}{z^2}\Bigr),\\
H_1(z+i0) &=& \frac{1}{24}\log^4{z} - \frac{3}{2} \zeta_2 \log^2{z}
            + \zeta_3 \log{z} + \frac{57}{16} \zeta_4 \nonumber\\
       &+& \frac{1}{z} \bigg( \frac{1}{6} \log^3{z} + \frac{1}{2} \log^2{z}
           - 3\zeta_2\log{z} + \zeta_3 - 3\zeta_2 -1 \bigg) \nonumber\\
       &+& i \pi\left\{ -\frac{1}{6}\log^3{z} + \zeta_2 \log{z} - \frac{7}{8}\zeta_3
        + \frac{1}{z} \bigg(- \frac{1}{2}\log^2{z} - \log{z} + \zeta_2
        \bigg) \right\}
        + {\cal O} \Bigl(\frac{1}{z^2}\Bigr),
\label{19}\\
H_2(z+i0) &=& -\frac{1}{24}\log^4{z} + \zeta_2 \log^2{z} - \frac{5}{8}\zeta_4 \nonumber \\
     &+& \frac{1}{z} \bigg( - \frac{1}{6} \log^3{z} - \frac{1}{2} \log^2{z}
       + (2\zeta_2 - 1)\log{z} + 2\zeta_2 -2 \bigg) \nonumber\\
     &+& i \pi \left\{ \frac{1}{6} \log^3{z} - \frac{3}{4} \zeta_3 
       + \frac{1}{z} \bigg( \frac{1}{2}\log^2{z} + \log{z} + 1\bigg) \right\}
       + {\cal O} \Bigl(\frac{1}{z^2}\Bigr),
\label{18}\\
H_3(z+i0) &=& -\frac{1}{24}\log^4{z} + \zeta_2 \log^2{z}
          - \frac{3}{2}\zeta_3 \log{z} - \frac{3}{16}\zeta_4 \nonumber\\
     &+& \frac{1}{z} \bigg( - \frac{1}{6}\log^3{z} + (2\zeta_2+1)\log{z}
         - \frac{3}{2}\zeta_3 + 1\bigg) \nonumber\\
       &+& i \pi \left\{ \frac{1}{6}\log^3{z} - \frac{1}{2}\zeta_2 \log{z}
      + \frac{7}{8}\zeta_3
       + \frac{1}{z} \bigg( \frac{1}{2} \log^2{z} - \frac{1}{2}\zeta_2 - 1 \bigg) \right\}
      + {\cal O} \Bigl(\frac{1}{z^2}\Bigr),
\label{20}\\
H_4(z+i0) &=& -\frac{1}{24}\log^4{z} + \zeta_2 \log^2{z}
           - \frac{3}{2}\zeta_3 \log{z} + \frac{7}{8}\zeta_4 \nonumber \\
      &+& \frac{1}{z} \bigg( \frac{1}{2} \log^2{z} + 2\log{z} - 2\zeta_2 + 2\bigg) \nonumber\\
      & + & i \pi \left\{ \frac{1}{6}\log^3{z} - \frac{1}{2}\zeta_2 \log{z}
            +\frac{3}{2}\zeta_3
         + \frac{1}{z} \bigg( - \log{z} -2 \bigg)
         \right\} + {\cal O} \Bigl(\frac{1}{z^2}\Bigr).
\label{17}
\end{eqnarray}
Finally we used the program \verb+hplog+ \cite{hplog} to check numerically 
the asymptotic behaviour  of the $H$-functions.


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

\section{Appendix C}

\label{App:C}

\def\theequation{C\arabic{equation}}
\setcounter{equation}0




  In this Appendix we consider the $H$-functions contributing to
the nonabelian part of the formfactor. 
For the definitions and recursive constructions of these functions
we refer to \cite{Bonciani1}. However, for completeness we give
here explicitly the definitions of the functions which appear in our
calculation. The following six new functions arise in the evaluation of the 
 two-loop nonabelian formfactor:
%\begin{equation}
%H_{-r}(-z),\, H_{-r,-r}(-z),\, H_{-r,-1}(-z),\, H_{0,-r,-r}(-z),\,
%    H_{-r,-r,-1}(-z),\, H_{0,-r,-r,-1}(-z)
%\end{equation}
\begin{eqnarray}
H_{-r}(z) \z=\z \int\limits_0^z \frac{dt_1}{\sqrt{t_1(t_1+4)}} \,,\\
H_{-r,-r}(z) \z=\z \int\limits_0^z \frac{dt_2}{\sqrt{t_2(t_2+4)}} 
                   \int\limits_0^{t_2} \frac{dt_1}{\sqrt{t_1(t_1+4)}} \,,\\
H_{-r,-1}(z) \z=\z \int\limits_0^z \frac{dt_2}{\sqrt{t_2(t_2+4)}} 
                   \int\limits_0^{t_2} \frac{dt_1}{1+t_1} \,,\\
H_{-r,-r,-1}(z) \z=\z \int\limits_0^z \frac{dt_3}{\sqrt{t_3(t_3+4)}} 
                   \int\limits_0^{t_3} \frac{dt_2}{\sqrt{t_2(t_2+4)}}
                   \int\limits_0^{t_2} \frac{dt_1}{1+t_1} \,,\\
H_{0,-r,-r}(z) \z=\z \int\limits_0^{z} \frac{dt_3}{1+t_3} 
                   \int\limits_0^{t_3} \frac{dt_2}{\sqrt{t_2(t_2+4)}} 
                   \int\limits_0^{t_2} \frac{dt_1}{\sqrt{t_1(t_1+4)}} \,,\\
H_{0,-r,-r,-1}(z) \z=\z \int\limits_0^{z} \frac{dt_4}{t_4} 
                   \int\limits_0^{t_4} \frac{dt_3}{\sqrt{t_3(t_3+4)}} 
                   \int\limits_0^{t_3} \frac{dt_2}{\sqrt{t_2(t_2+4)}} 
                   \int\limits_0^{t_2} \frac{dt_1}{1+t_1} \,.
\end{eqnarray}
In the formula (\ref{resf2}) for the nonabelian part the functions with odd
number of indices ``$-r$'' appear always with the factor
\begin{equation}
\label{gfactor}
  g(-z) = \frac{1}{\sqrt{1-4/z}} \,.
\end{equation}
It is easy to check that $H_{-r}$ and $H_{-r,-1}$ cannot be expanded 
in the Taylor series of small arguments (they have a branche point at zero), 
but the combinations $gH_{-r}$ and $gH_{-r,-1}$ can. 

  The integral representations given above are not very suitable
for the analysis and numerics. The ultimate task would be
to relate them to the usual (harmonic) polilogarithms. 
In order to do this one should choose a ``right'' variable. From the previous 
expirience \cite{FKV} it is known that for the diagrams,
posessing a branch point at $q^2=4m^2$, the appropriate variable
is given by ($z=q^2/m^2$)
\begin{equation}
  y = \frac{1-\sqrt{z/(z-4)}}{1+\sqrt{z/(z-4)}} \,.
\end{equation}

  In terms of $y$ the $g$-factor (\ref{gfactor}) is expressed as
\begin{equation}
  g(-z) = \frac{1-y}{1+y} \,,
\end{equation}
and the required $H$-fucntions take form
\begin{eqnarray}
%
H_{-r}(-z) \z=\z  - \log y \,,\\
%
H_{-r,-r}(-z) \z=\z \frac{1}{2} \log^2 y \,,\\
%
H_{-r,-1}(-z) \z=\z \frac{1}{2} \log^2 y  + \frac{1}{3} \Li_{2}(-y^3) 
                - \Li_{2}(-y) - \frac{1}{3}\zeta_2 \,,\\
%
H_{0,-r,-r}(-z) \z=\z  - \frac{1}{6} \log^3 y + \log(1-y) \log^2 y 
             - 2 \Li_{3}(y) + 2 \log y \,\Li_{2}(y) +  2 \zeta_3\,,\\
%
H_{-r,-r,-1}(-z) \z=\z - \frac{1}{6} \log^3 y +  \frac{1}{3}  \zeta_2 \log y 
              + \frac{2}{3}\zeta_3 + \Li_{3}(-y) -  \frac{1}{9} \Li_{3}(-y^3) 
\,,\\ 
H_{0,-r,-r,-1}(-z) \z=\z  \frac{1}{24} \log^4 y -  \frac{1}{6}  \zeta_2 
\log^2 y  - \frac{2}{3}\zeta_3 \log y    + \frac{89}{108}\zeta_4
 - \Li_{4}(-y) + \frac{1}{27} \Li_{4}(-y^3)  \nonumber \\
\z+\z 2 \Si_{1,3}(1-y) -  \frac{2}{3} \zeta_2 \Li_{2}(1-y)
  + 2\ln(1-y) \biggl( 
\frac{2}{3}\zeta_3 + \Li_{3}(-y) -  \frac{1}{9} \Li_{3}(-y^3) \biggr)
 \nonumber \\
\z+\z
 2\Li_{2}(y) \biggl( \Li_{2}(-y) -  \frac{1}{3} \Li_{2}(-y^3) \biggr)
+ 2 N_1(1)- 2N_1(y)
\,,
\end{eqnarray}
where
\begin{eqnarray}
N_1(y)= \int\limits^y_0 \frac{dt}{t} \Li_{2}(t) \ln(1-t+t^2), \qquad
N_1(1)= - \frac{11}{54} \zeta_4\,.
\end{eqnarray}

  As it is seen from the above fomulae, the $H$-functions with index ``$-r$''
can be rewritten in terms of harmonic polylogarithms but of {\it nonlinear}
argument $y$.

  In the limit when $z\to +\infty+i0$ we obtain
\begin{eqnarray}
%
\z\z g(-z)H_{-r}(-z) = 
     \lz{} + \frac{1}{z} \left( 2\lz{} - 2 \right) 
    + i\pi\left\{ - 1 - \frac{2}{z} \right\} + \Oz{2}  \,, \\
%
\z\z g(-z)H_{-r,-1}(-z) = 
     \frac{1}{2}\lz{2} - \frac{10}{3}\zeta_2
     + \left( \lz{2} - 2\lz{} - \frac{20}{3}\zeta_2 -1 \right) 
\nonumber\\
\z\z\qquad\qquad
     + i\pi\left\{ - \lz{} - \frac{1}{z}\left( 2\lz{} - 2 \right)\right\} + \Oz{2} \,,\\
%
\z\z H_{-r,-r}(-z) = 
     \frac{1}{2}\lz{2} - 3\zeta_2 - \frac{2\lz{}}{z}
     + i\pi\left\{ - \lz{} + \frac{2}{z} \right\} + \Oz{2} \,,\\
%
\z\z H_{0,-r,-r}(-z) = 
     \frac{1}{6}\lz{3} - 3\zeta_2\lz{} + 2\zeta_3 
      + \frac{1}{z} \left( 2\lz{} + 2 \right)  
\nonumber\\
\z\z\qquad\qquad
     + i\pi\left\{ - \frac{1}{2}\lz{2} + \zeta_2 - \frac{2}{z} \right\} + \Oz{2} \,,\\
%
\z\z H_{-r,-r,-1}(-z) = 
     \frac{1}{6}\lz{3} - \frac{10}{3}\zeta_2\lz{} + \frac{2}{3}\zeta_3 
      + \frac{1}{z} \left( - \lz{2} + \frac{20}{3}\zeta_2 + 1 \right)  
\nonumber\\
\z\z\qquad\qquad
     + i\pi\left\{ - \frac{1}{2}\lz{2} + \frac{4}{3}\zeta_2 
                   + \frac{2}{z}\lz{} \right\} + \Oz{2} \,,\\
%
%H_{0,-r,-r,-1}(-z) \z=\z  \frac{1}{24} \log^4 y -  \frac{1}{6}  \zeta_2 
%\log^2 y  - \frac{2}{3}\zeta_3 \log y    + \frac{125}{108}\zeta_4
%+2 N_1(1) \nonumber \\
%\z+\z y  \biggl( \frac{1}{3} \log^3 y - \log^2 y  
%+2 \Bigl(1-\frac{1}{3}\zeta_3 \Bigr) \log y  -1  + \frac{2}{3}\zeta_2
%- \frac{4}{3}\zeta_3 \biggr) + \Oz{2}  \,,\\
%
\z\z H_{0,-r,-r,-1}(-z) =
\frac{1}{24} \log^4{z} - \frac{5}{3} \zeta_2 \log^2{z} + \frac{2}{3} \zeta_3 \log{z} 
   + 7\zeta_4  \nonumber\\
\z\z\qquad\qquad
  + \frac{1}{z}\bigg( 
     \log^2{z} 
   + ( 2 + \frac{2}{3} \zeta_2 - \frac{2}{3} \zeta_3 )\log{z} 
    + 1 - \frac{20}{3} \zeta_2 
   \bigg) \nonumber\\
\z\z\qquad\qquad
   + i\pi\bigg\{
     - \frac{1}{6} \log^3{z}
     + \frac{4}{3} \zeta_2 \log{z}
     - \frac{2}{3} \zeta_3  \nonumber\\
\z\z\qquad\qquad 
     + \frac{1}{z} \bigg(
        - 2\log{z} -2 - \frac{2}{3}\zeta_2 + \frac{2}{3}\zeta_3 - 2\log{z} \bigg)
       \bigg\} + \Oz{2} \,.
\end{eqnarray}

  And finaly we give the values of $H$-functions at the particular 
point $z=1$:
\begin{eqnarray}
g(-1)H_{-r}(-1) \z=\z - \frac{\pi}{3\sqrt{3}} \,, \\
g(-1)H_{-r,-1}(-1)  \z=\z  \frac{2}{3} \frac{\Ls_2( {\textstyle \frac{\pi}{3}} )}{\sqrt{3}} \,, \\
H_{-r,-r}(-1) \z=\z  - \frac{1}{3} \zeta_2 \,, \\
H_{0,-r,-r}(-1) \z=\z \frac{4}{3} \zeta_3 
                      - \frac{2}{3} \pi \Ls_2( {\textstyle \frac{\pi}{3}} ) \,, \\
H_{-r,-r,-1}(-1) \z=\z \frac{1}{9} \zeta_3 \,, \\
H_{0,-r,-r,-1}(-1) \z=\z - \frac{7}{12} \zeta_4 
               + \frac{2}{3} \Bigl( \Ls_2( {\textstyle \frac{\pi}{3}} ) \Bigr)^2 \,,
\end{eqnarray}
%
where $\zeta_n$ is the Riemann $\zeta$-function and 
$\Ls_n(x)$ is the log-sine integral defined as
%
\begin{equation}
   \Ls_n(x) = -\int\limits_0^x \log^{n-1}\Bigl(2\sin\frac{t}{2}\Bigr) \, dt \,.
\end{equation}
In particular the constant $\Ls_2({\textstyle \frac{\pi}{3}})$,
sometimes denoted as Clausen's integral ${\rm Cl}_2({\textstyle \frac{\pi}{3}})$
(see, e.g., \cite{Lewin}), is given by
\begin{equation}
\Ls_2({\textstyle \frac{\pi}{3}})=1.014941606409653625\dots
\end{equation}

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

\begin{thebibliography}{99}

%\cite{Hamberg:1990np}
\bibitem{Hamberg:1990np}
  R.~Hamberg, W.~L.~van Neerven and T.~Matsuura,
  %``A Complete Calculation Of The Order Alpha-S**2 Correction To The Drell-Yan
  %K Factor,''
  Nucl.\ Phys.\ B {\bf 359} (1991) 343
  [Erratum-ibid.\ B {\bf 644} (2002) 403].
  %%CITATION = NUPHA,B359,343;%%

%\cite{Rijken:1994sh}
\bibitem{Rijken:1994sh}
  P.~J.~Rijken and W.~L.~van Neerven,
  %``Order alpha-s**2 contributions to the Drell-Yan cross-section at fixed
  %target energies,''
  Phys.\ Rev.\ D {\bf 51} (1995) 44.
%  [arXiv:hep-ph/9408366].
  %%CITATION = HEP-PH 9408366;%%

%\cite{Harlander:2003ai}
\bibitem{Harlander:2003ai}
R.~V.~Harlander and W.~B.~Kilgore,
%``Higgs boson production in bottom quark fusion at  next-to-next-to-leading
%order,''
Phys.\ Rev.\ D {\bf 68} (2003) 013001.
%[arXiv:hep-ph/0304035].
%%CITATION = HEP-PH 0304035;%%

%\cite{Baur:2001ze}
\bibitem{Baur:2001ze}
  U.~Baur, O.~Brein, W.~Hollik, C.~Schappacher and D.~Wackeroth,
  %``Electroweak radiative corrections to neutral-current Drell-Yan  processes
  %at hadron colliders,''
  Phys.\ Rev.\ D {\bf 65} (2002) 033007.
 % [arXiv:hep-ph/0108274].
  %%CITATION = HEP-PH 0108274;%%

%\cite{Czarnecki:1996ei}
\bibitem{Czarnecki:1996ei}
A.~Czarnecki and J.~H.~Kuhn,
%``Nonfactorizable QCD and electroweak corrections to the hadronic Z boson
%decay rate,''
Phys.\ Rev.\ Lett.\  {\bf 77} (1996) 3955.
%[arXiv:hep-ph/9608366].
%%CITATION = HEP-PH 9608366;%%

%\cite{Harlander:1997zb}
\bibitem{Harlander:1997zb}
R.~Harlander, T.~Seidensticker and M.~Steinhauser,
%``Complete corrections of O(alpha alpha(s)) to the decay of the Z boson  into
%bottom quarks,''
Phys.\ Lett.\ B {\bf 426} (1998) 125.
%[arXiv:hep-ph/9712228].
%%CITATION = HEP-PH 9712228;%%

%\cite{Fleischer:1999iq}
\bibitem{Fleischer:1999iq}
J.~Fleischer, F.~Jegerlehner, M.~Tentyukov and O.~Veretin,
%``Non factorizable O(alpha alpha(s)) corrections to the process  Z $\to$ b
%anti-b,''
Phys.\ Lett.\ B {\bf 459} (1999) 625.
%[arXiv:hep-ph/9904256].
%%CITATION = HEP-PH 9904256;%%

%\cite{Kuhn:2005az}
\bibitem{Kuhn:2005az}
  J.~H.~Kuhn, A.~Kulesza, S.~Pozzorini and M.~Schulze,
  %``One-loop weak corrections to hadronic production of Z bosons at large
  %transverse momenta,''
  Nucl.\ Phys.\ B {\bf 727} (2005) 368.
  %[arXiv:hep-ph/0507178].
  %%CITATION = HEP-PH 0507178;%%

\bibitem{Penin}
B.~Feucht, J.~H.~Kuhn, A.~A.~Penin and V.~A.~Smirnov,
%``Two-loop Sudakov form factor in a theory with mass gap,''
%arXiv:hep-ph/0404082.
Phys.\ Rev.\ Lett.\ {\bf 93} (2004) 101802.
%%CITATION = HEP-PH 0404082;%%

%\cite{Jantzen:2006jv}
\bibitem{Jantzen:2006jv}
  B.~Jantzen and V.~A.~Smirnov,
  %``The two-loop vector form factor in the Sudakov limit,''
  Eur.\ Phys.\ J.\  C {\bf 47} (2006) 671.
  %[arXiv:hep-ph/0603133].
  %%CITATION = EPHJA,C47,671;%%

\bibitem{FKV}
J.~Fleischer, A.~V.~Kotikov and O.~L.~Veretin,
%``Analytic two-loop results for selfenergy- and vertex-type diagrams  with one
%non-zero mass,''
Nucl.\ Phys.\ B {\bf 547} (1999) 343.
%[arXiv:hep-ph/9808242].
%%CITATION = HEP-PH 9808242;%%

%\cite{Kniehl:2005bc}
\bibitem{Kniehl:2005bc}
  B.~A.~Kniehl, A.~V.~Kotikov, A.~Onishchenko and O.~Veretin,
  %``Two-loop sunset diagrams with three massive lines,''
  Nucl.\ Phys.\ B {\bf 738} (2006) 306.
  %[arXiv:hep-ph/0510235].
  %%CITATION = HEP-PH 0510235;%%

\bibitem{hplog}
T.~Gehrmann and E.~Remiddi,
%``Numerical evaluation of two-dimensional harmonic polylogarithms,''
Comput.\ Phys.\ Commun.\  {\bf 144} (2002) 200.
%[arXiv:hep-ph/0111255].
%%CITATION = HEP-PH 0111255;%%

\bibitem{Gonsalves:1983nq}
R.~J.~Gonsalves,
%``Dimensionally Regularized Two Loop On-Shell Quark Form-Factor,''
Phys.\ Rev.\ D {\bf 28} (1983) 1542.
%%CITATION = PHRVA,D28,1542;%%

%\cite{vanNeerven:1985xr}
\bibitem{vanNeerven:1985xr}
  W.~L.~van Neerven,
  %``Dimensional Regularization Of Mass And Infrared Singularities In Two Loop
  %On-Shell Vertex Functions,''
  Nucl.\ Phys.\ B {\bf 268} (1986) 453.
  %%CITATION = NUPHA,B268,453;%%

%\cite{Kuhn:1999nn}
\bibitem{Kuhn:1999nn}
  J.~H.~Kuhn, A.~A.~Penin and V.~A.~Smirnov,
  %``Summing up subleading Sudakov logarithms,''
  Eur.\ Phys.\ J.\ C {\bf 17} (2000) 97.
  %[arXiv:hep-ph/9912503].
  %%CITATION = HEP-PH 9912503;%%

%\cite{Jantzen:2005xi}
\bibitem{Jantzen:2005xi}
  B.~Jantzen, J.~H.~Kuhn, A.~A.~Penin and V.~A.~Smirnov,
  %``Two-loop electroweak logarithms,''
  Phys.\ Rev.\  D {\bf 72} (2005) 051301
  [Erratum-ibid.\  D {\bf 74} (2006) 019901];\\
  %[arXiv:hep-ph/0504111].
  %%CITATION = PHRVA,D72,051301;%%
%\cite{Jantzen:2005az}
%\bibitem{Jantzen:2005az}
  B.~Jantzen, J.~H.~Kuhn, A.~A.~Penin and V.~A.~Smirnov,
  %``Two-loop electroweak logarithms in four-fermion processes at high
  %energy,''
  Nucl.\ Phys.\  B {\bf 731} (2005) 188
  [Erratum-ibid.\  B {\bf 752} (2006) 327].
  %[arXiv:hep-ph/0509157].


%\cite{Fadin:1999bq}
\bibitem{Fadin:1999bq}
  V.~S.~Fadin, L.~N.~Lipatov, A.~D.~Martin and M.~Melles,
  %``Resummation of double logarithms in electroweak high energy processes,''
  Phys.\ Rev.\ D {\bf 61} (2000) 094002.
  %[arXiv:hep-ph/9910338].
  %%CITATION = HEP-PH 9910338;%%

\bibitem{harmonicsums}
A.~Gonzalez-Arroyo, C.~Lopez and F.~J.~Yndurain,
%``Second Order Contributions To The Structure Functions In Deep Inelastic
%Scattering. I. Theoretical Calculations,''
Nucl.\ Phys.\ B {\bf 153} (1979) 161;\\
%%CITATION = NUPHA,B153,161;%%
%\cite{Gonzalez-Arroyo:1979he}
A.~Gonzalez-Arroyo and C.~Lopez,
%``Second Order Contributions To The Structure Functions In Deep Inelastic
%Scattering. 3. The Singlet Case,''
Nucl.\ Phys.\ B {\bf 166} (1980) 429;\\
%%CITATION = NUPHA,B166,429;%%
D. I. Kazakov and A.~V.~Kotikov, Nucl.\ Phys.\ B {\bf 307} (1988) 721;
%%CITATION = NUPHA,B307,721;%%
Theor.\ Math.\ Phys.\  {\bf 73} (1988) 1264
[Teor.\ Mat.\ Fiz.\  {\bf 73} (1987) 348].
%%CITATION = TMFZA,73,348;%%

\bibitem{harmonicsums1}
%\cite{Vermaseren:1998uu}
J.~A.~M.~Vermaseren,
%``Harmonic sums, Mellin transforms and integrals,''
Int.\ J.\ Mod.\ Phys.\ A {\bf 14} (1999) 2037;\\
%[arXiv:hep-ph/9806280].
%%CITATION = HEP-PH 9806280;%%
%\cite{Blumlein:1998if}
J.~Blumlein and S.~Kurth,
%``Harmonic sums and Mellin transforms up to two-loop order,''
Phys.\ Rev.\ D {\bf 60} (1999) 014018.
%[arXiv:hep-ph/9810241].
%%CITATION = HEP-PH 9810241;%%

\bibitem{Devoto}
%\cite{Devoto:1983tc}
A.~Devoto and D.~W.~Duke,
%``Table Of Integrals And Formulae For Feynman Diagram Calculations,''
Riv.\ Nuovo Cim.\  {\bf 7N6} (1984) 1.
%%CITATION = RNCIB,7N6,1;%%

\bibitem{Lewin}
L.~Lewin, {\it Polylogarithms and associated functions}
(North-Holland, Amsterdam, 1981).

\bibitem{HPLdefinition}
E.~Remiddi and J.~A.~M.~Vermaseren,
%``Harmonic polylogarithms,''
Int.\ J.\ Mod.\ Phys.\ A {\bf 15} (2000) 725.
%[arXiv:hep-ph/9905237].
%%CITATION = HEP-PH 9905237;%%

\bibitem{MultiplePL}
A.B.~Goncharov, Math. Res. Lett. {\bf 5} (1998) 497.

\bibitem{asymptotic}
F.V.~Tkachov, Preprint INR P-0332, Moscow (1983); P-0358, Moscow 1984;\\
K.G.~Chetyrkin,
Teor. Math. Phys. {\bf 75} (1988) 26; ibid {\bf 76} (1988) 207;
%%CITATION = TMPHA,75,346;%%
%%CITATION = TMPHA,76,809;%%
Preprint, MPI-PAE/PTh-13/91, Munich (1991);\\
~~V.A.~Smirnov,
Comm. Math. Phys.{\bf 134} (1990) 109;
%%CITATION = CMPHA,134,109;%%
{\it Renormalization and asymptotic expansions}
(Birkh\"auser, Basel, 1991);
{\it Applied asymptotic expansions in momenta and masses},
Berlin, Germany: Springer (2002), (Springer tracts in modern physics. 177).

\bibitem{FKV98}
J.~Fleischer, A.~V.~Kotikov and O.~L.~Veretin,
Phys.\ Lett.\ B {\bf 417} (1998) 163.
%%CITATION = HEP-PH 9707492;%%

%\cite{Harlander:2005rq}
\bibitem{Harlander:2005rq}
  R.~Harlander and P.~Kant,
  %``Higgs production and decay: Analytic results at next-to-leading order
  %QCD,''
  JHEP {\bf 0512} (2005) 015.
  %[arXiv:hep-ph/0509189].
  %%CITATION = HEP-PH 0509189;%%

\bibitem{DEM}
  A.V.~Kotikov,
  Phys.\ Lett.\ B {\bf 254} (1991) 158;
%%CITATION = PHLTA,B254,158;%%
  ibid. {\bf 259} (1991) 314;
%%CITATION = PHLTA,B259,314;%%
  ibid. {\bf 267} (1991) 123;\\
%%CITATION = PHLTA,B267,123;%%
 E. Remiddi,
  Nuovo Cim.\ A {\bf 110} (1997) 1435.
%  [arXiv:hep-ph/9711188].
%%CITATION = HEP-TH 9711188;%%

\bibitem{Bonciani1}
U.~Aglietti and R.~Bonciani,
%``Master integrals with one massive propagator for the two-loop  electroweak
%form factor,''
Nucl.\ Phys.\ B {\bf 668} (2003) 3;\\
%[arXiv:hep-ph/0304028].
%%CITATION = HEP-PH 0304028;%%
%\cite{Aglietti:2004nj}
U.~Aglietti, R.~Bonciani, G.~Degrassi and A.~Vicini,
%``Two-loop light fermion contribution to Higgs production and decays,''
%arXiv:hep-ph/0404071.
 Phys.\ Lett.\ B {\bf 595} (2004) 432;
%%CITATION = HEP-PH 0404071;%%
 ibid. {\bf 600} (2004) 57;
%%CITATION = HEP-PH 0407162;%%
JHEP {\bf 0701} (2007) 021.
%%CITATION = JHEPA,0701,021;%%

\bibitem{DIANA}
M.~Tentyukov and J.~Fleischer,
%``A Feynman diagram analyser DIANA,''
Comput.\ Phys.\ Commun.\  {\bf 132} (2000) 124.
%[arXiv:hep-ph/9904258].
%%CITATION = HEP-PH 9904258;%%

\bibitem{FORM}
J.A.M. Vermaseren, {\it Symbolic Manipulation with FORM},
Amsterdam, Computer Algebra, Netherland, 1991.

%\cite{Tentyukov:2004hz}
\bibitem{Tentyukov:2004hz}
  M.~Tentyukov, D.~Fliegner, M.~Frank, A.~Onischenko, A.~Retey, H.~M.~Staudenmaier and J.~A.~M.~Vermaseren,
  %``ParFORM: Parallel Version of the Symbolic Manipulation Program FORM,''
  arXiv:cs.sc/0407066.
  %%CITATION = CS-SC 0407066;%%

\bibitem{DavyKalmy}
A.~I.~Davydychev and M.~Y.~Kalmykov,
%``Massive Feynman diagrams and inverse binomial sums,''
%arXiv:hep-th/0303162.
Nucl.\ Phys.\ B {\bf 699} (2004) 3.
%%CITATION = HEP-TH 0303162;%%

%\cite{Grzadkowski:1986pm}
\bibitem{Grzadkowski:1986pm}
  B.~Grzadkowski, J.~H.~Kuhn, P.~Krawczyk and R.~G.~Stuart,
  %``ELECTROWEAK CORRECTIONS ON THE TOPONIUM RESONANCE,''
  Nucl.\ Phys.\ B {\bf 281} (1987) 18.
  %%CITATION = NUPHA,B281,18;%%

%\cite{Bohm:1986rj}
\bibitem{Bohm:1986rj}
  M.~Bohm, H.~Spiesberger and W.~Hollik,
  %``ON THE ONE LOOP RENORMALIZATION OF THE ELECTROWEAK STANDARD MODEL AND ITS
  %APPLICATION TO LEPTONIC PROCESSES,''
  Fortsch.\ Phys.\  {\bf 34} (1986) 687.
  %%CITATION = FPYKA,34,687;%%

%\cite{Bardin:1999ak}
\bibitem{Bardin:1999ak}
  D.~Y.~Bardin and G.~Passarino,
  {\it The standard model in the making: Precision study of the electroweak
  interactions},
%\href{http://www.slac.stanford.edu/spires/find/hep/www?irn=4320115}{SPIRES entry}
Oxford, UK: Clarendon (1999) 685 p, (International series of monographs on physics. 104).

%\cite{Berends:1987ab}
\bibitem{Berends:1987ab}
  F.~A.~Berends, W.~L.~van Neerven and G.~J.~H.~Burgers,
  %``Higher Order Radiative Corrections At Lep Energies,''
  Nucl.\ Phys.\ B {\bf 297} (1988) 429
  [Erratum-ibid.\ B {\bf 304} (1988) 921].
  %%CITATION = NUPHA,B297,429;%%



%\cite{Czakon:2004wm}
\bibitem{Czakon}
  M.~Czakon, J.~Gluza and T.~Riemann,
  %``Master integrals for massive two-loop Bhabha scattering in QED,''
  Phys.\ Rev.\  D {\bf 71} (2005) 073009;
 % [arXiv:hep-ph/0412164].
  %%CITATION = PHRVA,D71,073009;%%
%\cite{Czakon:2005jd}
%\bibitem{Czakon:2005jd}
%  M.~Czakon, J.~Gluza and T.~Riemann,
  %``Harmonic polylogarithms for massive Bhabha scattering,''
  Nucl.\ Instrum.\ Meth.\  A {\bf 559} (2006) 265.
%  [arXiv:hep-ph/0508212].
  %%CITATION = NUIMA,A559,265;%%


\end{thebibliography}

\end{document}

