%\documentclass[preprint,prl,aps,showpacs]{revtex4}
%\documentclass[twocolumn,prl,aps,showpacs,floatfix,endfloats]{revtex4}
%\documentclass[twocolumn,prl,aps,showpacs,floatfix]{revtex4}

\documentclass[twocolumn,prl,aps,showpacs]{revtex4}

%\usepackage{axodraw}
\usepackage{graphicx}%

\usepackage{dcolumn}
\usepackage{amsmath}
\usepackage{subfig}


\newcommand{\ed}{\end{document}}

\newcommand{\wtD}{\widetilde{D}}
\newcommand{\wtR}{\widetilde{R}}
\newcommand{\tdr}{\tilde{r}}
\newcommand{\tdd}{\tilde{d}}



\newcommand{\ice}[1]{\relax}

\newcommand{\re}[1]{(\ref{#1})}

\newcommand{\unl}{\underline}

\newcommand{\as}{a_s}
\newcommand{\api}{\frac{\alpha_s}{\pi}}
\newcommand{\als}{\alpha_s}
\newcommand{\beq}{\begin{equation}}

\newcommand{\ba}{\begin{array}}

\renewcommand{\arraystretch}{1.7}

\newcommand{\ea}{\end{array}}


\newcommand{\eeq}{\end{equation}}

\newcommand{\bea}{\begin{eqnarray}}

\newcommand{\eea}{\end{eqnarray}}
\newcommand{\g}{\gamma}

\newcommand{\asthree}{{}^{2\times{\cal O}(\as^3)}_{{\rm no}\;{\cal O}(\as^3)}}
\newcommand{\asfour}{{}^{2\times{\cal O}(\as^4)}_{{\rm no}\;{\cal O}(\as^4)}}


%belowfrom my SS-PLB paper of 97

\newcommand{\ds}{\displaystyle}
%\newcommand{\as}{\alpha_s}
\newcommand{\gm}{\gamma_m}
%\newcommand{\G}{\Gamma}
\newcommand{\ep}{\epsilon}

\newcommand{\EQN}{\label}
\newcommand{\ovl}{\overline}


%below is from SS.tex

\newcommand{\BreakO}{ \nonumber \\ &{}&   }

\newcommand{\Break}{ \right. \nonumber \\ &{}& \left. }
\newcommand{\BreakI}{ \right. \nonumber \\ &{}& \left. }
\newcommand{\BreakII}{ \right. \right.  \nonumber \\ &{}& \left.\left. }
\newcommand{\BreakIII}{ \right. \right. \right.  \nonumber \\ &{}& \left.\left.\left. }
\newcommand{\BreakIV}{
\right. \right. \right.\right.  \nonumber \\ &{}& \left. \left.\left.\left.
                     }


\newcommand{\lsb}{\left[}
\newcommand{\rsb}{\right]}

%%% for qcd constraints
\newcommand{\nn}{\nonumber}
\newcommand{\ra}{\rightarrow}

\newcommand{\cA}{{\cal A}}
\newcommand{\cB}{{\cal B}}
\newcommand{\cD}{{\cal D}}
\newcommand{\cF}{{\cal F}}
\newcommand{\cH}{{\cal H}}
\newcommand{\cK}{{\cal K}}
\newcommand{\cL}{{\cal L}}
\newcommand{\cM}{{\cal M}}
\newcommand{\cN}{{\cal N}}
\newcommand{\cO}{{\cal O}}
\newcommand{\cP}{{\cal P}}
\newcommand{\cT}{{\cal T}}
\newcommand{\cW}{{\cal W}}

\newcommand{\CF}{C_{\rm F}}
\newcommand{\Imm}{\mbox{\rm Im}}
\newcommand{\QCD}{\mbox{\rm {\tiny QCD}}}
\newcommand{\Tr}{\mbox{\rm Tr}}
\newcommand{\tr}{\mbox{\rm tr}}
\newcommand{\Li}{\mbox{\rm Li}}
\newcommand{\MeV}{\mbox{\rm MeV}}
\newcommand{\GeV}{\mbox{\rm GeV}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\newcommand{\dd}{{\rm d}}
\newcommand{\nnb}{\nonumber}


\begin{document}    

\mbox{\phantom{QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ}SFB/CPP-12-05, TTP11-31}

\title{
%%%%%%%%%%
Complete QCD Corrections to   Hadronic $Z$--Decays 
in Order $\alpha_s^4$
 }     

\author{P.~A.~Baikov}
\affiliation{
Skobeltsyn Institute of Nuclear Physics, Moscow State University,
Moscow 119991, Russia
        }
\author{K.~G.~Chetyrkin}
%\thanks{{\small Permanent address:
%Institute for Nuclear Research, Russian Academy of Sciences,
% Moscow 117312, Russia}}

\author{J.~H.~K\"uhn}
\affiliation{Institut f\"ur Theoretische Teilchenphysik,
  Universit\"at Karlsruhe, D-76128 Karlsruhe, Germany}


\author{J.Rittinger}
\affiliation{Institut f\"ur Theoretische Teilchenphysik,
  Universit\"at Karlsruhe, D-76128 Karlsruhe, Germany}

\begin{abstract}
%\vspace*{.7cm}
\noindent

Corrections of order $\alpha_s^4$ for
the axial singlet contributions for the decay rate of the $Z$-boson
into hadrons are evaluated in  the limit of the heavy  top quark  mass.
%
Combined with recently finished ${\cal O}(\alpha_s^4)$ calculations
of the non-signlet corrections, the new results directly  lead us
to  the first {\em complete}  ${\cal O}(\alpha_s^4)$ prediction for
the total hadronic decay rate of the Z-boson. 

The  new ${\cal O}(\alpha_s^4)$  term in Z-decay rate   lead to a
significant stabilization of the perturbative series, to a reduction
of the theory uncertainty in the strong coupling constant $\alpha_s$,
as extracted from these measurements, and to a small shift of the
central value.



\end{abstract}


\pacs{12.38.-t 12.38.Bx  12.20.-m }

\maketitle

%\section{Introduction}   

%\restylefloat{figure}

The precise determination of the $Z$-boson decay rate into hadrons at
LEP \cite{Alcaraz:2007ri} has led to one of the most precise determinations of
the strong coupling constant $\alpha_s$. From the experimental side, in
view of the fully inclusive nature of this measurement, the result is
fairly robust, in particular since it is insensitive to simulations of
the hadronic final state. Hence the error is essentially dominated by
the statistical uncertainty. From the theory side the advantage of the
measurement is its high energy, and as a result, the irrelevance of
nonperturbative and power-law suppressed terms. The smallness of
$\alpha_s$ at high energies then leads to a rapid decrease of higher
order corrections in the perturbative series and, correspondingly, to
a significant reduction of the theory error. 

A variety of methods has
been suggested to estimate the remaining uncertainty in the theory
prediction. Using the  last calculated term is probably the most
conservative approach, varying the renormalization scaBle $\mu$ within an
energy range characteristic for the problem 
(e.g. $M_Z/3 < \mu < 3\,M_Z$) is frequently
used, albeit with considerable ambiguity in the actual choice of the
region of the $\mu$-variation. In order to reduce the theoretical uncertainty in the
extraction of $\alpha_s$ to a level significantly smaller than the
experimental one, the knowledge of the corrections of 
${\cal O}(\alpha_s^4)$ is necessary. At the same time this calculation
opens the window for a considerable improvement in the
$\alpha_s$-determination at $Giga\,Z$, the project of a high-luminosity linear
collider operating at the $Z$-resonance (see e.g. \cite{Winter:2001av}). The dominant
part of the $\alpha_s^4$-corrections, the ''non-singlet''-piece, has
been evaluated in \cite{Baikov:2008jh}. This has lead  
to a slight shift of  the central value upwards 
from $0.1185\pm 0.0026$ to $0.1190\pm 0.0026$ \cite{Baikov:2008jh}
and a reduction of the theory error far below the error of 0.0026
from experiment. However, as noted already in \cite{Baikov:2008jh}, for a
complete evaluation of the decay rate in ${\cal O}(\alpha_s^4)$ an
additional set of corrections, namely those for the ``singlet''
contributions, is required. For the axial current correlator these start
at ${\cal O}(\alpha_s^2)$ \cite{Kniehl:1989qu,Kniehl:1989bb}, for the vector correlator at 
${\cal O}(\alpha_s^3)$. Both of them are presently known to third order in
$\alpha_s$ only 
\cite{Gorishnii:1990vf,Surguladze:1990tg,Chetyrkin:1993jm,Chetyrkin:1993ug,Larin:1993ju}. 
Hence, for a completely consistent ${\cal O}(\alpha_s^4)$
extraction of the strong coupling the extension of these results by one
order in $\alpha_s$ is required.

Before describing this calculation in detail, let us briefly recall 
the basic structure of QCD corrections to the correlator of the
electromagnetic and the neutral current, respectively, their similarities 
and their main differences. After splitting off inessential kinematic
factors, the absorptive part of the current-current correlator of the
electromagnetic current 
%$j^{\rm  em}_\mu=\sum_f Q_f\bar\psi_f\gamma_\mu\psi_f$ 
is expressed by the familiar $R$-ratio 
\begin{equation}
R^{\rm  em}=3\Big[\sum_f q_f^2 r_{\rm NS}^V + (\sum_f q_f)^2 r_{\rm S}^V\Big]~, \label {Rem}
\end{equation}
where $r_{\rm NS}^V$ and $r_{\rm S}^V$ stand for the (numerically dominant) non-singlet
and the singlet part respectively. The corresponding decomposition
for the correlator of the neutral current 
%$j^{\rm  nc}_\mu=j^V_\mu + j^A_\mu= \sum_f v_f\bar\psi_f\gamma_\mu\psi_f
% + \sum_f a_f\bar\psi_f\gamma_\mu\gamma_5 \psi_f$ 
involves the following four terms
\begin{equation}
R^{\rm  nc}=3\Big[\sum_f\!v_f^2 r_{\rm NS}^V + \big(\sum_f\!v_f\big)^2 r_{\rm S}^V
+\sum_f\!a_f^2 r_{\rm NS}^A + r_{\rm S;t,b}^A\Big]~,
\end{equation}
with $v_f\equiv 2I_f-4q_fs_W^2 $, $a_f \equiv 2I_f$ and $s_W$ defined as effective weak mixing angle. 
Here all but the top quark are assumed to be massless.  
%

(Mass corrections to both vector and axial vector correlator due to
other massive quarks are dominated by the bottom quark and can be
classified by orders in $m_b^2/M_Z^2$ and $\alpha_s$.
Up to  ${\cal O}(\als^2 m_b^2/M_Z^2)$ and
${\cal O}(\als^2 m_b^4/M_Z^4)$ they can be found in 
\cite{Chetyrkin:1996ia}, as well terms of  order 
$\als^2 m_b^2/M_Z^2$ (const + $\log\   m_b^2/M_Z^2$)
and  $\als^2 m_b^2/M_t^2$ (const + $\log\  m_b^2/M_Z^2$)
that arise  from  axial vector singlet contributions. Terms  of order 
$\als^3 m_b^4/M_Z^4$ and $\als^4 m_b^2/M_Z^2$ can be found in  
\cite{Chetyrkin:2000zk} and
\cite{Baikov:2004ku} respectively.
Corrections of order $\als^2 m_Z^2/m_t^2$ and  $\als^3 m_Z^2/m_t^2$ from 
singlet and non-singlet terms are known from \cite{Kniehl:1989qu,Kniehl:1989bb,Chetyrkin:1993tt} and \cite{Larin:1994va} respectively. 
These are important for the
actual $\alpha_s$-determination, but will not be discussed further in
the present paper.)





From the prefactors of the
non-singlet contributions in electromagnetic, vector and axial
correlator it is evident that different quark flavours contribute
incoherently, hence additive to the rate. Thus their contribution is significantly
enhanced in comparison with the singlet terms where amplitudes from
different flavours interfere destructively, with prefactors $(\sum_f
q_f)^2$ and $(\sum_f v_f)^2$ for the electromagnetic and neutral
current respectively.

Non-singlet contributions are present at the parton level
and the QCD corrections are known in second \cite{Chetyrkin:1979bj}, third
\cite{Gorishnii:1990vf,Surguladze:1990tg} and fourth \cite{Baikov:2008jh} order in $\alpha_s$. In terms of
Feynman diagrams, non-singlet contributions are characterized by the
fact that one quark loop connects the two external currents (Fig.\ref{PropNS}).
In the absorptive part of this fermion loop no top quark is present due to kinematic reasons,
whence the non-singlet functions are identical 
$r_{\rm NS}^V=r_{\rm NS}^A\equiv r_{\rm NS}$. 

In the case of singlet contributions of the vector current the two currents couple to two
different quark loops (Fig.\ref{PropSV}) requiring a three-gluon intermediate
state. Correspondingly the leading term is of ${\cal O}(\alpha_s^3)$
and has been obtained already long time ago \cite{Gorishnii:1990vf,Surguladze:1990tg}. The NLO
corrections to this   result are of ${\cal O}(\alpha_s^4)$. They serve
to soften the strong scale dependence of the ${\cal O}(\alpha_s^3)$  result, stabilize the theory prediction and will be the
subject of this paper.

\begin{figure}[b!]
%\begin{figure}[H]
%\centering
\subfloat[]{
 \includegraphics[width=.25\linewidth]{PropNS.eps}
\label{PropNS}
}
\subfloat[]{
 \includegraphics[width=.3\linewidth]{PropSV.eps}
\label{PropSV}
}
\subfloat[]{
 \includegraphics[width=.3\linewidth]{PropSA.eps}
\label{PropSA}
}
\caption{Different contributions to $r$-ratios: (a) non-singlet, (b) vector singlet and (c) axial vector singlet.}
\end{figure}

The situation is different in the case of the singlet axial vector 
current correlator. The axial couplings of the two members of an isospin 
doublet are opposite equal. Hence their singlet contributions vanishes, if the
corresponding quark masses are equal. This approximation is valid for
the two lightest quark doublets. The only remaining contribution
originates from the combination of bottom and top quarks with their 
specific mass hierarchy  $m_b^2<<M_Z^2<<m_t^2$ (Fig.\ref{PropSA}). In this case
the contribution starts at ${\cal O}(\alpha_s^2)$ and is further
enhanced by the ''large'' logarithm $\log (m_t^2/M_Z^2)$ 
\cite{Kniehl:1989qu,Kniehl:1989bb}. Corrections of ${\cal O}(\alpha_s^3)$
have been calculated in \cite{Chetyrkin:1993jm,Chetyrkin:1993ug,Larin:1993ju}, those of  ${\cal O}(\alpha_s^4)$
will be the subject of this paper.

The evaluation of the NLO terms of $r_V^{\rm S}$ requires the
calculation of the absorptive parts of five-loop diagrams with
massless propagators which, with the help of some complicated
combinatorics based on the $R^*$-operation \cite{Chetyrkin:1984xa}, can be boiled down
to the calculation of  four-loop propagator
diagrams. The  latter have been computed  
via reduction to 28  master integrals, based on evaluating  
sufficiently many terms of the $1/D$ expansion \cite{Baikov:2005nv} of
the corresponding coefficient functions \cite{Baikov:1996rk}. 
This direct procedure required huge computing resources and
was performed using a parallel version \cite{Tentyukov:2004hz} of FORM
\cite{Vermaseren:2000nd}. The master integrals are reliably  known from
\cite{Baikov:2010hf,Smirnov:2010hd,Lee:2011jt}. 
The details of the calculation, the results in analytic form and their 
relation to the Gross-Llewellyn Smith sum rule will be given in  \cite{Baikov:2012}.



The evaluation of the NNLO terms of $R_{\rm S;t,b}^A$ involves again
absorptive parts of five-loop diagrams with massless propagators,
however, in addition also absorptive parts of four-loop diagrams combined
with one-loop massive tadpoles, etc. down to one-loop massless
diagrams together with four-loop massive tadpoles. 
The methods employed in our calculations  will be described in more detail in
\cite{Baikov:2012}, together with the results expressed in terms of
group theoretical coefficients and transcendental numbers.


The result is valid
in the limit $M_Z^2 \ll 4\,M_t^2$, an excellent approximation as evident
from the lower orders.
The relative importance of the various terms
is best seen from the results for the various $r$-ratios introduced
above, expressed in numerical form
\begin{align}
 r_{\rm NS}=&1 + \as + 1.40923\,\as^2 - 12.7671\,\as^3 - 79.9806\,\as^4~, \nnb \\
 r_{\rm S}^V=&-0.41318\,\as^3 - 4.9841\,\as^4~, \nnb \\
 r_{\rm S:t,b}^A =\, &  (-3.08333+l_t)\,\as^2  \\ &+ (-15.9877+3.72222\,l_t+1.91667\,l_t^2)\,\as^3 \nnb \\
		 &~\hspace{-1cm} +( 49.162-17.6822\,l_t+14.7153\,l_t^2 + 3.67361\,l_t^3)\,\as^4  ~, \nnb
\end{align}
with $\as=\alpha_s(M_Z)/\pi$ and $l_t=\ln(M_Z^2/M_t^2)$. Since all
three $r$-ratios are seperately scale invariant, the corresponding results
for a generic value of $\alpha_s(\mu)$ can easily be reconstructed. 
Using the pole mass $M_t=172$ GeV the axial singlet contribution is given
in numerical form by
\begin{align}
 r_{\rm S;t,b}^A = -4.35248\,\as^2 - 17.6245\,\as^3 + 87.7962\,\as^4~.
\end{align}
Collecting now all QCD terms, the decay rate of the $Z$-boson into hadrons 
can be cast into the following form
\begin{align}
 \Gamma_Z=\Gamma_0\,R^{\rm nc}=\tfrac{G_F\,M_Z^3}{24\pi\sqrt{2}}\,R^{\rm nc}~.
\end{align}
Here all electroweak corrections are assumed to be collected in the
prefactor $\Gamma_0$, and the forementioned mass corrections are
ignored as well as electroweak and mixed QCD-electroweak corrections 
\cite{Czarnecki:1996ei,Harlander:1997zb,Czarnecki:1996ac}. Thus the $R$-ratio is now known up to $\mathcal{O}(\as^4)$
\begin{align}
 R^{\rm nc}=&20.1945 + 20.1945\,\as \nnb \\ &+ (28.4587-13.0575+0)\,\as^2 \nnb \\ &+(-257.825-52.8736-2.12068)\,\as^3 \nnb \\ &+(-1615.17+263.388-25.5814)\,\as^4~,
\end{align}
with $s_W^2=0.231$. The three terms in the brackets display separately non-singlet, axial singlet and vector singlet contributions.

Let us now evaluate the impact of the newly calculated terms on the
$\alpha_s$-determination from $Z$-decays. Following our approach for
the non-singlet terms (where a shift $\delta\alpha_s=0.0005$ had been
obtained \cite{Baikov:2008jh}, consistent with an analysis \cite{Quast}  based on
results of the electroweak working group \cite{Alcaraz:2007ri} and confirmed by
the G-fitter collaboration \cite{Goebel:2010ux}), we consider the
quantity $R^{\rm nc}$ as ``pseudo-observable''. With a starting value
$R^{\rm nc}=20.9612$, if evaluated for $\alpha_s=0.1190$ and without the
$\alpha_s^4$ singlet terms, a shift $\delta\alpha_s=-0.00008$
is obtained after including the newly calculated
contributions.

\begin{figure}[b!]
\centering
\subfloat[]{
 \hspace{.45cm}\includegraphics[width=.95\linewidth,height=5.4cm]{rns.eps}
\label{rns}
} \\
\subfloat[]{
 \includegraphics[width=\linewidth,height=5.4cm]{rsv.eps}
\label{rsv}
} \\
\subfloat[]{
 \hspace{.2cm}\includegraphics[width=.98\linewidth,height=5.4cm]{rsa.eps}
\label{rsa}
}
\caption{Scale dependence of (a) non-singlet $r_{\rm NS}$, (b) vector singlet $r_S^V$ and (c) axial vector singlet $r^A_{\rm S;t,b}$. Dotted, dash-dotted, dashed and solid curves refer to $\mathcal{O}(\alpha_s)$ up to $\mathcal{O}(\alpha_s^4)$ predictions.
$\alpha_s (M_Z) =0.1190$ and $n_l =5 $ is adopted in all these curves.
}
\label{Plots}
\end{figure}

As discussed in \cite{Baikov:2008jh}, the non-singlet $\alpha_s^4$
term leads to a considerable stabilization of the theory prediction,
and, correspondingly, to a reduction of the theory error. A similar
statement holds true for the singlet contribution. To illustrate this
aspect, the dependence on the renormalization scale $\mu$ is shown in
Fig. \ref{Plots} for $r_{\rm NS}$, $r_{\rm S}^V$ and $r^A_{\rm S;t,b}$. The relative
variation is significantly reduced in all three cases. In particular
for the vector singlet case we observe a shift of the result by about
a factor 1.45 (for $\mu=M_Z$) and a considerable flattening of the
result. Using for example the Principle of Minimal Sensitivity (PMS)
\cite{Stevenson:1981vj} as a guidance for the proper choice of scale,
$\mu=0.3\,M_Z$ seems to be favoured, leading to an amplification of
the LO result by a factor 1.68 (if the latter is evaluated for
$\mu=M_Z$, as done traditionally). On the other hand LO and NLO
predictions, if evaluated for $\mu=0.3\,M_Z$, are fairly close which
is another argument in favor of this choice. 

Let us assume that the remaining theory uncertainties from $r_{\rm NS}$, $r^V_{\rm S}$ and
$r^A_{\rm S;t,b}$ can be estimated by varying $\mu$
between $M_Z/3$ and $3\,M_Z$ and using the maximal variation as twice
the uncertainty $\delta r$. This leads to $\delta\Gamma_{\rm NS}=0.0971$
MeV, $\delta\Gamma^V_{\rm S}=0.0026$ MeV and $\delta\Gamma^A_{\rm S}=0.046$
MeV. Even adding these terms linearly, they are far below the
experimental error of $\delta\Gamma_{exp}=2.3$ MeV \cite{Nakamura:2010zzi}.
% and, in total lead to a combined theory error of: $\delta\alpha_s=0.1189\pm0.0026$. 
In combination with the quadratic and quartic mass terms, which are
known to $\mathcal{O}(\alpha_s^4)$ and $\mathcal{O}(\alpha_s^3)$
respectively, this analysis completes the QCD corrections to the $Z$
decay rate.

Let us also comment on the impact of the $\alpha_s^4$ singlet result
on the measurement of $R^{\rm em}$ at low energies, i.e. in the region
accessible at BESS or at  B-factories, say between 3 GeV and 10
GeV. Considering the large luminosities collected at these machines, a
precise $\alpha_s$ determination from $R^{\rm em}$ seems possible
\cite{Chetyrkin:1996tz}. In the 
low energy region only $r_{\rm S}^V$ and $r_{\rm NS}^V$ contribute. Since
$\sum_{f=u,d,s}\,q_f=0$, the singlet contribution vanishes in the
three flavour case. If we consider the region above charm and below
bottom threshold, say at 10 GeV, only $u,d,s$ and $c$ quarks
contribute, the relative weight of the $r_S^V$ in eq. (\ref{Rem}) is
given by $(\sum q_f)^2/(\sum q_f^2) = 2/5$, and thus is fairly
suppressed. At energy of 10 GeV, in the absence of open bottom quark
contribution, it seems appropriate to analyze the results in an
effective four flavour theory with 
\[
r^V_S = -0.41318 \,\as^3(\mu)  - (5.1757 + 2.5824 \ln \mu^2/s)\,\as^4(\mu)
{}.
\]
As shown in Fig. \ref{Plot10gev}, it is evident that the
scale dependence is softened in NLO. Again a scale $\mu$ around
$0.3\,\sqrt{s}$ is suggested by PMS. Note, that in this case LO and
NLO results differ by nearly a factor two; nevertheless the
contribution to $R^{\rm em}$ remains  small.

\begin{figure}[h!]
\centering
 \includegraphics[width=\linewidth]{rsv10gev.eps}
\label{rsv10}
\caption{Scale dependence of the vector singlet $r_S^V$ around 10 GeV. Dashed and solid curves refer to $\mathcal{O}(\alpha_s^3)$ and $\mathcal{O}(\alpha_s^4)$ predictions.
$n_l =4 $ and $\alpha_s (10 \mbox{GeV}) = 0.1806$ as obtained with the use of  package RunDec \protect\cite{Chetyrkin:2000yt}
have  been assumed.
}
\label{Plot10gev}
\end{figure}

In conclusion we want to mention that all our calculations have been
performed on a SGI ALTIX 24-node IB-interconnected cluster of 8-cores
Xeon computers and on the HP XC4000 supercomputer of the federal state
Baden-W\"urttemberg
%(claster of Dual Core AMD Opteron 2.6 GH),
using  parallel  MPI-based \cite{Tentyukov:2004hz} as well as thread-based \cite{Tentyukov:2007mu} versions  of FORM
\cite{Vermaseren:2000nd}.  For evaluation of color factors we have used the FORM program { COLOR}
\cite{vanRitbergen:1998pn}. The diagrams have been generated with QGRAF \cite{Nogueira:1991ex}.


This work was supported by the Deutsche Forschungsgemeinschaft in the
Sonderforschungsbereich/Transregio SFB/TR-9 ``Computational Particle
Physics'' and  by RFBR grants   11-02-01196 and  10-02-00525.

\begin{thebibliography}{99}
%
% GZfull4_ref.tex -- generated by sortref-2.3.6  
% ((C) R. Harlander, http://www.robert-harlander.de/software/)
% on Wed Jan 25 18:15:07 CET 2012
%

%1

\ice{
\bibitem{Alcaraz:2007ri}
  J.~Alcaraz {\it et al.}  [LEP Collaborations and ALEPH Collaboration and
                  DELPHI Collaboration an],
  %``Precision Electroweak Measurements and Constraints on the Standard Model,''
  arXiv:0712.0929 [hep-ex].
  %%CITATION = ARXIV:0712.0929;%%

}

\bibitem{Alcaraz:2007ri}
  J.~Alcaraz {\it et al.} [ LEP and ALEPH and DELPHI and L3 and OPAL and LEP Electroweak Working Group Collaborations ],
  [arXiv:0712.0929 [hep-ex]].
  %%CITATION = ARXIV:0712.0929;%%


%2
\bibitem{Winter:2001av}
  M.~Winter, 
  ECFA-DESY Linear Collider Note LC-PHSM-2001-016.

%3
\bibitem{Baikov:2008jh}
  P.~A.~Baikov, K.~G.~Chetyrkin and J.~H.~K\"uhn,
  Phys.\ Rev.\ Lett.\  {\bf 101}, 012002 (2008)
  [arXiv:0801.1821 [hep-ph]].
  %%CITATION = PRLTA,101,012002;%%

%4
\bibitem{Kniehl:1989qu}
  B.~A.~Kniehl and J.~H.~K\"uhn,
less  Nucl.\ Phys.\  B {\bf 329}, 547 (1990).
  %%CITATION = NUPHA,B329,547;%%

%5
\bibitem{Kniehl:1989bb}
  B.~A.~Kniehl and J.~H.~K\"uhn,
  Phys.\ Lett.\  B {\bf 224}, 229 (1989).
  %%CITATION = PHLTA,B224,229;%%

%6
\bibitem{Gorishnii:1990vf}
  S.~G.~Gorishnii, A.~L.~Kataev and S.~A.~Larin,
  Phys.\ Lett.\  B {\bf 259}, 144 (1991).
  %%CITATION = PHLTA,B259,144;%%

%7
\bibitem{Surguladze:1990tg}
  L.~R.~Surguladze and M.~A.~Samuel,
  Phys.\ Rev.\ Lett.\  {\bf 66}, 560 (1991)
  [Erratum-ibid.\  {\bf 66}, 2416 (1991)].
  %%CITATION = PRLTA,66,560;%%

%8
\bibitem{Chetyrkin:1993jm}
  K.~G.~Chetyrkin and J.~H.~K\"uhn,
  Phys.\ Lett.\  B {\bf 308}, 127 (1993).
  %%CITATION = PHLTA,B308,127;%%

%9
\bibitem{Chetyrkin:1993ug}
  K.~G.~Chetyrkin and O.~V.~Tarasov,
  Phys.\ Lett.\  B {\bf 327}, 114 (1994)
  [arXiv:hep-ph/9312323].
  %%CITATION = PHLTA,B327,114;%%

%10
\bibitem{Larin:1993ju}
  S.~A.~Larin, T.~van Ritbergen, J.~A.~M.~Vermaseren,
  Phys.\ Lett.\  {\bf B320}, 159-164 (1994).
  [hep-ph/9310378].

%11
\bibitem{Chetyrkin:1996ia}
  K.~G.~Chetyrkin, J.~H.~K\"uhn and A.~Kwiatkowski,
  Phys.\ Rept.\  {\bf 277} (1996) 189 and and references therein.
  %%CITATION = PRPLC,277,189;%%

%12
\bibitem{Chetyrkin:2000zk}
  K.~G.~Chetyrkin, R.~V.~Harlander and J.~H.~K\"uhn,
  Nucl.\ Phys.\  B {\bf 586} (2000) 56
  [Erratum-ibid.\  B {\bf 634} (2002) 413]
  [arXiv:hep-ph/0005139].
  %%CITATION = NUPHA,B586,56;%%

%13
\bibitem{Baikov:2004ku}
  P.~A.~Baikov, K.~G.~Chetyrkin, J.~H.~K\"uhn,
  Nucl.\ Phys.\ Proc.\ Suppl.\  {\bf 135}, 243-246 (2004).

%14
\bibitem{Chetyrkin:1993tt}
  K.~G.~Chetyrkin,
  Phys.\ Lett.\  B {\bf 307} (1993) 169.
  %%CITATION = PHLTA,B307,169;%%

%15
\bibitem{Larin:1994va}
  S.~A.~Larin, T.~van Ritbergen and J.~A.~M.~Vermaseren,
  Nucl.\ Phys.\  B {\bf 438} (1995) 278
  [arXiv:hep-ph/9411260].
  %%CITATION = NUPHA,B438,278;%%

%16
\bibitem{Chetyrkin:1979bj}
  K.~G.~Chetyrkin, A.~L.~Kataev and F.~V.~Tkachov,
  Phys.\ Lett.\  B {\bf 85}, 277 (1979).
  %%CITATION = PHLTA,B85,277;%%

%17
\bibitem{Chetyrkin:1984xa}
  K.~G.~Chetyrkin and V.~A.~Smirnov,
  Phys.\ Lett.\  B {\bf 144} (1984) 419.
  %%CITATION = PHLTA,B144,419;%%

%18
\bibitem{Baikov:2005nv}
  P.~A.~Baikov,
  Phys.\ Lett.\  B {\bf 634} (2006) 325
  [arXiv:hep-ph/0507053].
  %%CITATION = PHLTA,B634,325;%%

%19
\bibitem{Baikov:1996rk}
  P.~A.~Baikov,
  Phys.\ Lett.\  B {\bf 385} (1996) 404
  [arXiv:hep-ph/9603267].
  %%CITATION = PHLTA,B385,404;%%

%20
\bibitem{Tentyukov:2004hz}
  M.~Tentyukov, D.~Fliegner, M.~Frank, A.~Onischenko, A.~Retey, H.~M.~Staudenmaier and J.~A.~M.~Vermaseren,
  arXiv:cs/0407066.
  %%CITATION = CS/0407066;%%

%21
\bibitem{Vermaseren:2000nd}
  J.~A.~M.~Vermaseren,
  arXiv:math-ph/0010025.
  %%CITATION = MATH-PH/0010025;%%

%22
\bibitem{Baikov:2010hf}
  P.~A.~Baikov and K.~G.~Chetyrkin,
  Nucl.\ Phys.\  B {\bf 837} (2010) 186
  [arXiv:1004.1153 [hep-ph]].
  %%CITATION = NUPHA,B837,186;%%

%23
\bibitem{Smirnov:2010hd}
  A.~V.~Smirnov and M.~Tentyukov,
  Nucl.\ Phys.\  B {\bf 837} (2010) 40
  [arXiv:1004.1149 [hep-ph]].
  %%CITATION = NUPHA,B837,40;%%

%24
\bibitem{Lee:2011jt}
  R.~N.~Lee, A.~V.~Smirnov and V.~A.~Smirnov,
  Nucl.\ Phys.\  B {\bf 856} (2012) 95
  [arXiv:1108.0732 [hep-th]].
  %%CITATION = NUPHA,B856,95;%%
%

%25
\bibitem{Baikov:2012}
  P.~A.~Baikov, K.~G.~Chetyrkin, J.~H.~K\"uhn and J.~Rittinger,
  in preparation.
  %%CITATION = PHLTA,B315,179;%%

%26
\bibitem{Czarnecki:1996ei}
  A.~Czarnecki and J.~H.~K\"uhn,
  Phys.\ Rev.\ Lett.\  {\bf 77}, 3955 (1996)
  [arXiv:hep-ph/9608366].
  %%CITATION = PRLTA,77,3955;%%

%27
\bibitem{Harlander:1997zb}
  R.~Harlander, T.~Seidensticker and M.~Steinhauser,
  Phys.\ Lett.\  B {\bf 426}, 125 (1998)
  [arXiv:hep-ph/9712228].
  %%CITATION = PHLTA,B426,125;%%

%28
\bibitem{Czarnecki:1996ac}
  A.~Czarnecki and K.~Melnikov,
  Phys.\ Rev.\  D {\bf 56}, 1638 (1997)
  [arXiv:hep-ph/9609302].
  %%CITATION = PHRVA,D56,1638;%%

%29
\bibitem{Quast}
  G.~Quast,
  private communication.

%30
\bibitem{Goebel:2010ux}
  M.~Goebel  [on behalf of the Gfitter Group],
  PoS {\bf ICHEP2010}, 570 (2010)
  [arXiv:1012.1331 [hep-ph]].
  %%CITATION = POSCI,ICHEP2010,570;%%

%31

\bibitem{Stevenson:1981vj}
  P.~M.~Stevenson,
  %``Optimized Perturbation Theory,''
  Phys.\ Rev.\  D {\bf 23} (1981) 2916.
  %%CITATION = PHRVA,D23,2916;%%


%32
%\cite{Nakamura:2010zzi}
\bibitem{Nakamura:2010zzi}
  K.~Nakamura {\it et al.}  [Particle Data Group],
  %``Review of particle physics,''
  J.\ Phys.\ G {\bf 37} (2010) 075021.
  %%CITATION = JPHGB,G37,075021;%%



%33
\bibitem{Chetyrkin:2000yt}
  K.~G.~Chetyrkin, J.~H.~K\"uhn and M.~Steinhauser,
  Comput.\ Phys.\ Commun.\  {\bf 133} (2000) 43
  [arXiv:hep-ph/0004189].
  %%CITATION = CPHCB,133,43;%%

%34
\bibitem{Chetyrkin:1996tz}
  K.~G.~Chetyrkin, J.~H.~K\"uhn and T.~Teubner,
  Phys.\ Rev.\  D {\bf 56} (1997) 3011
  [arXiv:hep-ph/9609411].
  %%CITATION = PHRVA,D56,3011;%%

%35
\bibitem{Tentyukov:2007mu}
  M.~Tentyukov and J.~A.~M.~Vermaseren,
  Comput.\ Phys.\ Commun.\  {\bf 181} (2010) 1419
  [arXiv:hep-ph/0702279].
  %%CITATION = CPHCB,181,1419;%%

%36
\bibitem{vanRitbergen:1998pn}
  T.~van Ritbergen, A.~N.~Schellekens and J.~A.~M.~Vermaseren,
  Int.\ J.\ Mod.\ Phys.\  A {\bf 14} (1999) 41
  [arXiv:hep-ph/9802376].
  %%CITATION = IMPAE,A14,41;%%

%37
\bibitem{Nogueira:1991ex}
  P.~Nogueira,
  J.\ Comput.\ Phys.\  {\bf 105} (1993) 279.
  %%CITATION = JCTPA,105,279;%%


%\input{GZfull4_ref}
%\input{bib}

\end{thebibliography}


\end{document}

