%Title: THREE-LOOP HEAVY QUARK VACUUM POLARIZATION
%Authors: K.G. Chetyrkin, J.H. Kuehn, M. Steinhauser
%Published: Proceedings of the 1996 Zeuthen Workshop on <EM>Elementary Particle Theory, QCD and QED in Higher Orders, Nucl. Phys.</EM><B>B</B> (Proc.Suppl.)<B>51C</B>(1996)66-70.   
%%%%%%%%%% espcrc2.tex %%%%%%%%%%
\documentstyle[twoside,fleqn,espcrc2,epsf]{article}

% put your own definitions here:
%   \newcommand{\cZ}{\cal{Z}}
%   \newtheorem{def}{Definition}[section]
%   ...
\newcommand{\ttbs}{\char'134}
\newcommand{\AmS}{{\protect\the\textfont2
  A\kern-.1667em\lower.5ex\hbox{M}\kern-.125emS}}


\voffset -3.35cm
\hoffset +0.5cm
\renewcommand{\topfraction}{0.95}
\renewcommand{\bottomfraction}{0.95}
\renewcommand{\textfraction}{0.02}



% add words to TeX's hyphenation exception list
\hyphenation{author another created financial paper re-commend-ed}

% declarations for front matter
\title{THREE-LOOP HEAVY QUARK VACUUM POLARIZATION}


\author{K.G.~Chetyrkin\address{Max-Planck-Institut f\"ur Physik, 
                               Werner-Heisenberg-Institut,
                               D-80805 Munich, Germany\\
                   and\\
         Institute for Nuclear Research, Russian Academy of Sciences,
         Moscow 117312, Russia}
   J.H.~K\"uhn, 
   M.~Steinhauser\address{Institut f\"ur Theoretische Teilchenphysik,
                         Universit\"at Karlsruhe, 
                         D-76128 Karlsruhe, Germany}
        }


\begin{document}

\begin{abstract}
The real and imaginary part of the vacuum polarization function
$\Pi(q^2)$ induced by a massive quark is calculated in perturbative
QCD up to order $\alpha_s^2$.
We combine the information from small and large momentum
region and from the threshold using conformal mapping
and Pad\'e approximation. This leads us to formulae for
$\Pi(q^2)$ valid for arbitrary $m^2/q^2$. Taking subsequently
the imaginary part we get the ${\cal O}(\alpha_s^2)$ to
$R  \equiv \sigma(e^+ e^- \to \mbox{hadrons})/
           \sigma(e^+ e^- \to \mu^+ \mu^-)$.
This extends the calculation by
K\"all\'en and Sabry from two to three loops.
\end{abstract}

% typeset front matter (including abstract)
\maketitle

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

\section{Introduction}

The measurement of the total cross section for
electron positron annihilation into hadrons allows for
a unique test of pertubative QCD. The decay rate
$\Gamma(Z \to \mbox{hadrons})$ provides one of the
most precise determinations of the strong coupling 
constant $\alpha_s$. In the high energy limit the quark masses
can often be neglected.
In this approximation QCD corrections to 
$R  \equiv \sigma(e^+ e^- \to \mbox{hadrons})/
           \sigma(e^+ e^- \to \mu^+ \mu^-)$
have been calculated up to order $\alpha_s^3$
\cite{CheKatTka79DinSap79CelGon80,GorKatLar91SurSam91}.
For precision measurements the dominant mass corrections must be
included through an expansion in $m^2/s$. Terms up to order 
$\alpha_s^3 m^2/s$ 
\cite{CheKue90}
and $\alpha_s^2 m^4/s^2$ 
\cite{CheKue94}
are available at present,
providing an acceptable approximation from the high energy
region down to intermediate energy values.
For a number of measurements, however, the information on the 
complete mass dependence is desirable. This
includes charm and bottom meson production above the resonance 
region, say
above $4.5$~GeV and $12$~GeV, respectively, and, of course,
top quark production at a future electron positron collider.

To order $\alpha_s$ this calculation was performed by
K\"all\'en and Sabry in the context of QED a long time ago
\cite{KaeSab55}.
With measurements of ever increasing precision, predictions
in order $\alpha_s^2$ are needed for a reliable
comparison between theory and experiment. Furthermore, 
when one tries to apply the ${\cal O}(\alpha)$ result 
to QCD, with its running coupling
constant, the choice of scale becomes important.
In fact, the distinction between the two intrinsically different 
scales, the relative momentum versus the center of mass
energy, is crucial for a stable numerical prediction.
This information can be obtained from a full calculation 
to order $\alpha_s^2$ only. 
Such a calculation then allows to predict the cross section 
in the complete energy region where perturbative QCD can be applied
--- from close to threshold up to high energies.

In this contribution results for the cross section are presented in 
order $\alpha_s^2$.
They are obtained from the vacuum polarization $\Pi(q^2)$ which is
calculated up to three loops.
The imaginary part of the ``fermionic contribution'' ---  derived
from diagrams with a massless quark 
loop inserted in the gluon propagator --- has been calculated in
\cite{HoaKueTeu95}.
All integrals could be performed to the end and the result was expressed
in terms of polylogarithms. 
In this paper the calculation is extended to the full set of 
diagrams relevant for QCD.
Instead of trying to perform the integrals analytically, we use information
of $\Pi(q^2)$ from the large $q^2$ behaviour, the expansion around
$q^2=0$ and from threshold.
Only results without renormalization group improvement 
and resummation of the Coulomb singularities from higher orders  
are presented. Resummation of leading higher order terms, phenomenological
applications and a more detailed discussion of our methods will be
presented elsewhere.

\section {Outline of the Calculation}

The different behaviour at threshold makes it necessary to decompose
$\Pi$ according to its colour structure. It is convenient to
write:
\begin{eqnarray}
\Pi(q^2) &=& \Pi^{(0)}(q^2) 
         + \frac{\alpha_s(\mu^2)}{\pi}C_F\Pi^{(1)}(q^2)
\nonumber\\
&&
         + \left(\frac{\alpha_s(\mu^2)}{\pi}\right)^2\Pi^{(2)}(q^2)
         + \cdots,
\\
\Pi^{(2)} &=& 
                C_F^2       \Pi_{\mbox{\scriptsize\it A}}^{(2)}
              + C_A C_F     \Pi_{\mbox{\scriptsize\it NA}}^{(2)}
              + C_F T n_l \Pi_{\mbox{\scriptsize\it l}}^{(2)}
\nonumber\\
&&
              + C_F T     \Pi_{\mbox{\scriptsize\it F}}^{(2)}.
\end{eqnarray}
The same notation is adopted to the
physical observable $R(s)$ which is related to $\Pi(q^2)$ by
\begin{eqnarray}
R(s)   &=&  12\pi\, \mbox{Im}\Pi(q^2=s+i\epsilon).
\end{eqnarray}

The contributions from diagrams with $n_l$ light 
or one massive 
internal fermion loop are denoted
by $C_F T n_l\Pi_l^{(2)}$ and $C_F T \Pi_F^{(2)}$, respectively.
The purely gluonic corrections
are proportional to $C_F^2$ or $C_A C_F$ where the former are the only
contributions in an abelian theory and the latter are characteristic for
the nonabelian aspects of QCD.

All steps described below have been performed seperately for
the first three contributions to $\Pi^{(2)}$. 
In fact, new information is only
obtained for $\Pi_{\mbox{\scriptsize\it A}}^{(2)}$ 
and $\Pi_{\mbox{\scriptsize\it NA}}^{(2)}$ since
$\mbox{Im}\Pi_{\mbox{\scriptsize\it l}}^{(2)}$ 
is already known analytically \cite{HoaKueTeu95}.
The contribution from a four particle cut 
with threshold at $4m$ is given
in terms of a two dimensional integral
\cite{HoaKueTeu95}
which can be solved 
easily numerically,
so $\Pi_F^{(2)}$ will not be treated in this paper.
Also the contributions to 
$\Pi(q^2)$ and $R(s)$ which originate from 
diagrams where the electromagnetic current couples to a massless
quark and the massive quark is produced through a virtual
gluon 
will not be discussed here. 
They have been calculated in 
\cite{HoaKueTeu94}.
In the following the behaviour of $\Pi(q^2)$ in the three different 
kinematical regions and the approximation method is discussed.

\vspace{.5em}
\noindent
{\bf Analysis of the high $q^2$ behaviour}
\vspace{.5em}

The high energy behaviour of $\Pi$ provides important 
constraints on the complete answer.
In the limit of small $m^2/q^2$ the constant term and the one
proportional to $m^2/q^2$ (modulated by powers of $\ln \mu^2/q^2$) have been
calculated a long time ago
\cite{GorKatLar86}.
For the imaginary part even the $m^4/q^4$ terms are available 
\cite{CheKue94}.
This provides an important test of the numerical results 
presented below.

\vspace{.5em}
\noindent
{\bf Threshold behaviour}
\vspace{.5em}

General arguments based on the influence of Coulomb exchange close to 
threshold, combined with the information on the perturbative QCD
potential and the running of $\alpha_s$ dictate the singularities
and the structure of the leading cuts close to threshold, that
is for small $v=\sqrt{1-4m^2/s}$.
The $C_F^2$ term 
is directly related to the QED result with internal photon lines only.
The leading $1/v$
singularity  and the constant term  of $R_A$
can be predicted from the nonrelativistic
Greens function for the Coulomb potential
and the ${\cal O}(\alpha_s)$ calculation. 
The next-to-leading
term is determined by the combination of one loop
results again with the Coulomb singularities
\cite{BarGatKoeKun75,VolSmi94,BaiBro95}. One finds
\begin{eqnarray}
R_{\mbox{\scriptsize\it A}}^{(2)} &=& 
3\left(\frac{\pi^4}{8v} - 3\pi^2 + \ldots\right).
\label{Ra}
\end{eqnarray}

The contributions $\sim C_A C_F$ and $\sim C_F T n_l$ can be treated
in parallel. 
For these colour structures the perturbative QCD potential
\cite{Fis77}
\begin{eqnarray}
 V_{\mbox{\scriptsize QCD}}(\vec{q}\,^2) &=& 
         -4\pi C_F\frac{\alpha_V(\vec{q}\,^2)}{\vec{q}\,^2},
\\
\alpha_V(\vec{q}\,^2)  &=&   \alpha_s(\mu^2)\Bigg[
      1 + \frac{\alpha_s(\mu^2)}{4\pi}\bigg(
          \left(\frac{11}{3}C_A
\right.
\nonumber\\
&&
\left.
\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
-\frac{4}{3}T n_l\right)
          \left(-\ln\frac{\vec{q}\,^2}{\mu^2}+\frac{5}{3}\right)
          -\frac{8}{3}C_A            \bigg)         
\label{alphav}
\Bigg]
\end{eqnarray}
will become important. 
The leading $C_AC_F$ and $C_F T n_l$ term in $R$ is proportional
to $\ln v$ and is responsible for the evolution of the
coupling constant close to threshold. Also the constant term can
be predicted by the observation, that the leading term in
order $\alpha_s$ is induced by the potential.
The ${\cal O}(\alpha_s)$ result
\begin{equation}
R=3\frac{v(3-v^2)}{2}
\left(1
       + C_F\frac{\pi^2(1+v^2)}{2v}\frac{\alpha_s}{\pi}+\ldots
     \right)
\end{equation}
is employed to predict the logarithmic and constant $C_FC_A$ and
$C_FTn_l$ terms of ${\cal O}(\alpha_s^2)$ by replacing $\alpha_s$
by $\alpha_V(4\vec{p}\,^2=v^2 s)$ as given in Eq.(\ref{alphav}).
This implies the following threshold behaviour:
\begin{eqnarray}
R_{\it NA}^{(2)}&=&3\frac{\pi^2}{3}(3-v^2)(1+v^2)
\nonumber\\
&&
                         \left(
                              -\frac{11}{16}\ln\frac{v^2 s}{\mu^2}
                              +\frac{31}{48}
                     +\ldots    \right),
\label{Rna}
\\
R_l^{(2)}&=&3\frac{\pi^2}{3}(3-v^2)(1+v^2)       
\nonumber\\
&&
                         \left(   
                                  \frac{1}{4}\ln\frac{v^2 s}{\mu^2}
                                 -\frac{5}{12}
                     +\ldots    \right).
\label{Rnl}
\end{eqnarray}

This ansatz can be verified for the $C_F T n_l$ term since in this case 
the result is known in analytical form
\cite{HoaKueTeu95}. 
Extending the ansatz from the behaviour of the imaginary part
close to the branching point into the complex plane
allows to predict the leading term of
$\Pi(q^2)$ $\sim \ln v$ and $\sim \ln^2 v$.  

\vspace{.5em}
\noindent
{\bf Behaviour at $q^2=0$}
\vspace{.5em}

Important information is contained in the Taylor series of $\Pi(q^2)$
around zero. The calculation of the first seven nontrivial terms is 
based on the evaluation of three-loop tadpole integrals with 
the help of the algebraic program MATAD written in FORM
\cite{VerFORM}
which
performs the traces,
calculates the derivatives with respect to the external momenta. 
It reduces the large number of different
integrals to one master integral and a few simple ones
through an elaborate set of
recursion relations based on the integration-by-parts method
\cite{CheTka81,Bro92}.
The result can be written in the form:
\begin{eqnarray}
\Pi^{(2)} &=& 
              \frac{3}{16\pi^2}
              \sum_{n>0} C_{n}^{(2)} \left(\frac{q^2}{4m^2}\right)^n,
\end{eqnarray}
where the first seven moments are listed in
\cite{CheKueSte96POL}.

\vspace{.5em}
\noindent
{\bf Conformal mapping and Pad\'e approximation}
\vspace{.5em}

The vacuum polarization function 
$\Pi^{(2)}$ is analytic in the complex plane
cut from $q^2=4m^2$ to $+\infty$. The Taylor series in $q^2$ is
thus convergent in the domain $|q^2|<4m^2$ only. The
conformal mapping 
which corresponds to the variable 
transformation
\begin{eqnarray}
\omega = \frac{1-\sqrt{1-q^2/4m^2}}{1+\sqrt{1-q^2/4m^2}},
\,
\frac{q^2}{4m^2} = \frac{4\omega}{(1+\omega)^2},
\label{omega}
\end{eqnarray}
transforms the cut complex $q^2$ plane into the
interior of the
unit circle. The special points
$q^2=0,4m^2,-\infty$ correspond to $\omega=0,1,-1$, respectively.

The upper (lower) part of the cut is mapped onto the upper (lower) 
perimeter of the circle.
The Taylor series in $\omega$ thus converges in the interior of the
unit circle. To obtain predictions for
$\Pi(q^2)$ at the boundary it has been suggested
\cite{FleTar94,BroFleTar93}
to use
the Pad\'e approximation 
which converges towards $\Pi(q^2)$ even on the perimeter.
To improve the accuracy 
the singular threshold behaviour  
and the large $q^2$ behaviour
is incorporated into simple analytical functions
which are removed 
from $\Pi^{(2)}$ before the Pad\'e approximation is 
performed.
The quality of this
procedure can be tested by comparing the prediction with
the known result for $\mbox{Im}\Pi_{\mbox{\scriptsize\it l}}^{(2)}$.

The logarithmic singularities at threshold and large $q^2$
are removed by subtraction, the $1/v$ singularity, which is present
for the $C_F^2$ terms only, by multiplication with $v$ as 
suggested in
\cite{BaiBro95}.
The imaginary part of the remainder which is actually 
approximated by the Pad\'e method is thus smooth in the
entire circle, numerically small and vanishes at 
$\omega=1$ and $\omega=-1$. 

\begin{figure}[ht]
 \leavevmode
 \begin{center}
 \begin{tabular}{c}
   \epsfxsize=7.0cm
   \epsffile[110 330 460 520]{plots/pimorcf22.ps}  
   \\
   \epsfxsize=7.0cm
   \epsffile[110 330 460 520]{plots/pimorcacf2.ps} 
   \\
   \epsfxsize=7.0cm
   \epsffile[110 330 460 520]{plots/pimornf2.ps}
 \end{tabular}
 \caption{\label{fighigh}Complete results (full line) are
                         compared to the threshold
                         approximations and the high energy
                         approximations including the $m^2/s$
                         (dash-dotted) and the $m^4/s^2$ (dashed) terms
                         ($x=2m/\protect\sqrt{s}$).
         }
 \end{center}
\end{figure}

\begin{figure}[ht]
 \leavevmode
 \begin{center}
 \begin{tabular}{c}
   \epsfxsize=7.0cm
   \epsffile[110 330 460 520]{plots/pilocf23.ps}  
   \\
   \epsfxsize=7.0cm
   \epsffile[110 330 460 520]{plots/pilocacf3.ps} 
   \\
   \epsfxsize=7.0cm
   \epsffile[110 330 460 520]{plots/pilonf3.ps} 
 \end{tabular}
 \caption{\label{figthr}The 
                        threshold behaviour of the remainder $\delta R$ 
                        for three different Pad\'e approximants is shown.
                        (The singular
                        and constant parts around threshold are subtracted.)
        }
 \end{center}
\end{figure}

\section{Results}

After performing the Pad\'e approximation for the smooth remainder
with $\omega$ as natural variable, the transformation (\ref{omega})
is inverted and the full vacuum polarization function reconstructed
by reintroducing the threshold and high energy terms. This
procedure provides real and imaginary parts of $\Pi^{(2)}$.
Subsequently only the absorbtive part of $\Pi^{(2)}$ (multiplied by
$12\pi$) will be presented.

In Fig.~\ref{fighigh} the complete results are shown 
for $\mu^2=m^2$ with 
$R_{\mbox{\scriptsize\it A}}^{(2)}$, $R_{\mbox{\scriptsize\it NA}}^{(2)}$
and $R_{\mbox{\scriptsize\it l}}^{(2)}$
displayed separately. 
The solid line represents the full correction. The threshold
approximation is given  by the dotted curve. In the
high energy region besides the corrections containing the $m^2/s$ terms
(dash-dotted line) also the quartic (dashed line) approximations are shown.
It should be stressed that they are not incorporated into the 
construction of $R^{(2)}$ but they are evidently very well 
reproduced by the method presented here.

Different Pad\'e approximations of the same degree 
and approximants with a
reduced number of parameters give rise to practically indentical
predictions, which could hardly be distinguished in 
Fig.~\ref{fighigh}. Minor variations are observed close
to threshold, {\it after} subtracting the singular and constant parts.
The remainder $\delta R$ for up to ten different Pad\'e 
approximants is shown in 
Fig.~\ref{figthr}. 
There is a perfect agreement for $R^{(2)}_l$. It is hard to detect the
exact result which is represented by the dotted line.
$R_{\it NA}^{(2)}$ seems to converge to the solid line 
($[4/4], [5/3]$ and $[3/5]$)
when more moments from small $q^2$ are included. The dashed lines are from
the $[3/3], [4/2], [2/4]$ and $[3/4]$, the dotted
ones from lower order Pad\'e approximants.
The dash-dotted curve is the [4/3] Pad\'e approximant and has a pole very 
close to $\omega=1 (1.07\ldots)$.
For the abelian part a classification of the different results
can be seen: the dashed lines are $[4/2]$ and $[2/4]$,
the solid ones $[3/2], [2/3], [5/3], [3/5]$ and $[5/4]$ 
Pad\'e approximants.

To summarize: Real and imaginary part of the vacuum polarization
function $\Pi(q^2)$ from a massive quark have been calculated up to
three loops for QCD and QED. This result extends the classic
calculations  of K\"all\'en and Sabri \cite{KaeSab55} to 
next-to-leading order. The imaginary part
can be used to predict the cross section for production of
massive quarks for arbitrary $m^2/s$, wherever perturbative QCD can be 
justified --- from above the quarkonium resonance region up to high 
energies.

\vspace{5ex}                  
\noindent
{\bf Acknowledgments}

\noindent
     Work supported by BMFT under Contract 056KA93P6, 
     DFG under Contract Ku502/6-1 and INTAS under Contract INTAS-93-0744
     and the Graduiertenkolleg Elementarteilchenphysik, Karlsruhe.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\sloppy
\raggedright
\def\app#1#2#3{{\it Act. Phys. Pol. }{\bf B #1} (#2) #3}
\def\apa#1#2#3{{\it Act. Phys. Austr.}{\bf #1} (#2) #3}
\def\lhc{Proc. LHC Workshop, CERN 90-10}
\def\npb#1#2#3{{\it Nucl. Phys. }{\bf B #1} (#2) #3}
\def\plb#1#2#3{{\it Phys. Lett. }{\bf B #1} (#2) #3}
\def\prd#1#2#3{{\it Phys. Rev. }{\bf D #1} (#2) #3}
\def\pR#1#2#3{{\it Phys. Rev. }{\bf #1} (#2) #3}
\def\prl#1#2#3{{\it Phys. Rev. Lett. }{\bf #1} (#2) #3}
\def\prc#1#2#3{{\it Phys. Reports }{\bf #1} (#2) #3}
\def\cpc#1#2#3{{\it Comp. Phys. Commun. }{\bf #1} (#2) #3}
\def\nim#1#2#3{{\it Nucl. Inst. Meth. }{\bf #1} (#2) #3}
\def\pr#1#2#3{{\it Phys. Reports }{\bf #1} (#2) #3}
\def\sovnp#1#2#3{{\it Sov. J. Nucl. Phys. }{\bf #1} (#2) #3}
\def\jl#1#2#3{{\it JETP Lett. }{\bf #1} (#2) #3}
\def\jet#1#2#3{{\it JETP Lett. }{\bf #1} (#2) #3}
\def\zpc#1#2#3{{\it Z. Phys. }{\bf C #1} (#2) #3}
\def\ptp#1#2#3{{\it Prog.~Theor.~Phys.~}{\bf #1} (#2) #3}
\def\nca#1#2#3{{\it Nouvo~Cim.~}{\bf #1A} (#2) #3}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{thebibliography}{99}

\bibitem{CheKatTka79DinSap79CelGon80}
K.G. Chetyrkin, A.L. Kataev and F.V. Tkachov,
\prl{85}{1979}{227};\\
M. Dine and J. Sapirstein,
\prl{43}{1979}{668};\\
W. Celmaster and R.J. Gonsalves,
\prl{44}{1980}{560}.

\bibitem{GorKatLar91SurSam91}
S.G. Gorishny, A.L. Kataev and S.A. Larin,
\plb{259}{1991}{144};\\
L.R. Surguladze and M.A. Samuel,
\prl{66}{1991}{560};
erratum ibid, 2416.

\bibitem{CheKue90}
K.G. Chetyrkin and J.H. K\"uhn,
\plb{248}{1990}{359}.

\bibitem{CheKue94}
K.G. Chetyrkin and J.H. K\"uhn,
\npb{432}{1994}{337}.

\bibitem{KaeSab55}
G. K\"all\'en and A. Sabry, K. Dan. Vidensk. Selsk. Mat.-Fys. Medd. 
{\bf 29} (1955) No. 17,\\ 
see also 
J. Schwinger, Particles, {\it Sources and Fields}, Vol.II,
   (Addison-Wesley, New York, 1973).

\bibitem{HoaKueTeu95}
A.H. Hoang, J.H. K\"uhn and T. Teubner,
\npb{452}{1995}{173}.

\bibitem{HoaKueTeu94}
A.H. Hoang, M. Je\.zabek, J.H. K\"uhn and T. Teubner,
\plb{338}{1994}{330}.

\bibitem{GorKatLar86}
S.G. Gorishny, A.L. Kataev and S.A. Larin,
\nca{92}{1986}{119}.

\bibitem{Fis77}
W. Fischler, \npb{129}{1977}{157};\\
A. Billoire, \plb{92}{1980}{343}.

\bibitem{BarGatKoeKun75}
R. Barbieri, R. Gatto, R. K\"ogerler and Z. Kunszt, \plb{57}{1975}{455}.

\bibitem{VolSmi94}
B. H. Smith and M.B. Voloshin,
\plb{324}{1994}{117}; erratum ibid, {\bf B 333} (1994) 564.

\bibitem{BaiBro95}
P.A. Baikov and D.J. Broadhurst, Report Nos. OUT-4102-54, INP-95-13/377 and
hep-ph/9504398, to appear in {\it New Computing Techniques in Physics
Research IV}, (World Scientific, in press).

\bibitem{VerFORM}
J.A.M. Vermaseren, {\it Symbolic Manipulation with FORM},
(Computer Algebra Netherlands, Amsterdam, 1991).

\bibitem{CheTka81}
F.V. Tkachov, \plb{100}{1981}{65};\\
K.G. Chetyrkin and F.V. Tkachov, \npb{192}{1981}{159}.

\bibitem{Bro92}
D.J. Broadhurst, \zpc{54}{1992}{54}.

\bibitem{CheKueSte96POL}
K.G. Chetyrkin, J.H. K\"uhn and M. Steinhauser,
\plb{371}{1996}{93};\\
K.G. Chetyrkin, J.H. K\"uhn and M. Steinhauser, in preparation.

\bibitem{FleTar94}
J. Fleischer and O.V. Tarasov, \zpc{64}{1994}{413}.

\bibitem{BroFleTar93}
D.J. Broadhurst, J. Fleischer and O.V. Tarasov, \zpc{60}{1993}{287}.

\end{thebibliography}

\end{document}







