%
\subsection{QCD corrections to the photon polarization
  function\label{sec::polfunc}}
%
\vspace{2ex}
{\bf Notation and methods}\\[2ex]
%
This section is concerned with the computation of QCD corrections
to the photon propagator. We define the polarization function in the
following way:
\begin{eqnarray}
\left(-g_{\mu\nu}q^2+q_\mu q_\nu\right)\,\Pi(q^2)
&=&i\int \dd^4 x\,e^{iqx}\langle 0 |Tj_\mu(x) j_\nu(0)|0 \rangle
\,,
\end{eqnarray}
where $j^\mu=\bar{\psi}\gamma^\mu\psi$ is the diagonal vector current of
two quarks with mass $m$.  The main motivation is the simple connection
of $\Pi(q^2)$ to the physical quantity $R(s)$ which is defined as the
normalized total cross section to the production of heavy quarks:
\begin{eqnarray}
R(s)&\equiv&\frac{\sigma\left(e^+e^-\to\mbox{hadrons}\right)}
            {\sigma\left(e^+e^-\to\mu^+\mu^-\right)}
    \,\,=\,\,12\pi\mbox{Im}\Pi(q^2=s+i\epsilon)
\,.
\end{eqnarray}
The advantage of making a detour to $\Pi(q^2)$ and not considering
$R(s)$ from the very beginning is the implicit summation of real
radiation and infra-red singularities. The fact that the number of loops
to be evaluated for $\Pi(q^2)$ is larger by one is made up for by having
to deal with two-point functions only, for which a huge amount of
technology is available.

It is convenient to separately consider the following contributions:
\begin{eqnarray}
\Pi(q^2) &=& \Pi^{(0)}(q^2) 
         + \frac{\alpha_s(\mu^2)}{\pi} C_F \Pi^{(1)}(q^2)
         + \left(\frac{\alpha_s(\mu^2)}{\pi}\right)^2\Pi^{(2)}(q^2)
         + \ldots\,\,,
\\
\Pi^{(2)} &=&
                C_F^2       \Pi_A^{(2)}
              + C_A C_F     \Pi_{\it NA}^{(2)}
              + C_F T   n_l \Pi_l^{(2)}
              + C_F T       \Pi_F^{(2)},
\label{eqpi2}
\end{eqnarray}
and similarly for $R(s)$.  The colour factors $C_F=(N_c^2-1)/(2N_c)$
and $C_A=N_c$ correspond to the Casimir operators of the fundamental
and adjoint representations of $SU(3)$, respectively.  The case of QCD
corresponds to $N_c=3$, the trace normalization of the fundamental
representation is $T=1/2$, and the number of light (massless) quark
flavours is denoted by $n_l$.  In Eq.~(\ref{eqpi2}), $\Pi_A^{(2)}$ is
the Abelian contribution (corresponding to quenched QED) and 
$\Pi_{\it NA}^{(2)}$ is
the non-Abelian part specific for QCD.  There are two fermionic
contributions arising from diagrams with two closed fermion lines,
so-called double-bubble diagrams: For $\Pi_l^{(2)}$ the quark in the
inner loop is massless, the one in the outer loop massive, whereas for
$\Pi_F^{(2)}$ both fermions carry mass $m$.  The case where the external
current couples to massless quarks and these via gluons to massive ones
will not be addressed here.  Its imaginary part was considered in
\cite{HoaJezKueTeu94}.  The Feynman diagrams contributing to one-, two-
and three-loop order are shown in Fig.~\ref{figpol}.

%
\begin{figure}
  \begin{center}
    \leavevmode
    \epsfxsize=10.cm
    \epsffile[90 430 510 730]{figpol.ps}
    \hfill
    \parbox{\captionwidth}{
    \caption[]{\label{figpol}\sloppy
      Diagrams contributing to the one-, two- and three-loop polarization
      function. Solid lines represent quarks, loopy ones are gluons.
      }}
  \end{center}
\end{figure}
%

The two-loop corrections of ${\cal O}(\alpha_s)$ were computed in
analytic form in the context of QED quite some time ago~\cite{KalSab55}.
In subsequent works the calculation was redone and more convenient
representations were found~\cite{BarRem73,Kni90,BroFleTar93}.  It is, at
least with the currently available technology, out of range to compute
the three-loop diagrams in complete analytic form for arbitrary $q$ and
$m$.  Only for a subclass of diagrams, namely the ones with two closed
fermion lines where one of them is massless, was it possible to compute
the corresponding contribution to $R(s)$
analytically~\cite{HoaJezKueTeu94,HoaKueTeu95,Teudiss,Hoadiss}.
Nevertheless, there are essentially two approaches which provide very
good approximations for $\Pi^{(2)}(q^2)$, both of them relying heavily
on the use of computer algebra.

One method is to consider the polarization function $\Pi(q^2)$ in the
limit $q^2\gg m^2$. Application of the large-momentum procedure (see
Section~\ref{subasymp}) leads to an expansion in $m^2/q^2$ with the
coefficients still depending on logarithms $\ln(-q^2/m^2)$.  The idea is
to evaluate as many terms as possible (the limitations essentially
coming from the CPU time and the size of the intermediate expressions)
and thus to approximate the true function $\Pi(q^2)$ even for rather
small values of $q^2$.

Whereas the above method results in analytical expressions, an
alternative way is to construct a semi-numerical result for $\Pi(q^2)$
which will, however, be valid for all values of $q$ and $m$. In addition
to the high energy expansion discussed above, this method also requires
the expansion of $\Pi(q^2)$ in the limit $q^2\ll m^2$.  Together
with some information about the (two-particle) threshold $q^2=4m^2$, a
suitable conformal mapping and the use of Pad\'e approximation, one
obtains an approximation for $\Pi(q^2)$ over the full kinematical region
\cite{CheKueSte96,CheHarSte98}.

Referring to the literature for physical applications (e.g.,
\cite{ckk96,chkst98}), in the following we will concentrate on the
technical aspect of the calculation in the limits $q^2\ll m^2$ and
$q^2\gg m^2$. While the latter case directly corresponds to the first
method above, the former one is needed for the second method.

As can be seen from Fig.~\ref{figpol} only a small number of diagrams
contribute at three-loop level. This makes the use of a diagram generating
program like {\tt QGRAF} unnecessary.  The main challenge instead is to
compute as many terms in the expansions as possible.  Apart from fast
computers this requires an efficient algebraic language and optimized
programs to deal with a large number of terms. The choice for this
problem was {\tt MATAD} and {\tt MINCER}, both written in {\tt FORM}
(see Section~\ref{sec::matmin}).

%
\vspace{4ex}
{\bf Large mass limit}\\[2ex]
%
In the limit $q^2\ll m^2$ application of the hard-mass procedure (see
Section~\ref{subasymp}) shows that it suffices to keep the naive Taylor expansion
of the diagrams in their external momentum. Thus one stays with the
calculation of massive tadpole integrals. Let us consider the moments $C_n$ of
the polarization function, defined through the Taylor series
\begin{eqnarray}
\Pi^{(2)}(q^2) &=& 
              \frac{3}{16\pi^2}
              \sum_{n>0} C_{n}^{(2)} \left(\frac{q^2}{4m^2}\right)^n
\,,
\label{eqpolfunpi}
\end{eqnarray}
where $m$ is the on-shell mass.
Although this series does not develop an imaginary part, one may gain
information on the rate from it by exploiting the analyticity of
$\Pi(q^2)$. We do not want to discuss this approach here and refer the
interested reader to~\cite{BaiBro95,CheKueSte96}.  At three-loop level
the evaluation of the coefficients up to $C_8$ was performed in
\cite{CheKueSte96}. The calculation required disk space of the order of
several GB for {\tt FORM} (so-called ``formswap'') to store the
intermediate expressions.  The total CPU time on a 256~MHz DEC-Alpha
workstation with 128~MB main memory was roughly two weeks.  With the
input for the ladder-type diagram already shown in
Section~\ref{sec::matmin}, let us now have a look at the output of {\tt
  MATAD}. The expansion up to terms of order $(q^2/m^2)^4$ reads:
\begin{verbatim}
   ladder =
       + ep^-3 * (  - 8/9*Q.Q )

       + ep^-2 * (  - 1/693*Q.Q^5*M^-8 + 8/945*Q.Q^4*M^-6 - 2/35*Q.Q^3*M^-4 + 
         8/15*Q.Q^2*M^-2 - 134/27*Q.Q )

       + ep^-1 * (  - 20*M^2 - 4042699/108056025*Q.Q^5*M^-8 + 273004/1488375*
         Q.Q^4*M^-6 - 63029/66150*Q.Q^3*M^-4 + 10658/2025*Q.Q^2*M^-2 - 473/81*
         Q.Q - 4/3*Q.Q*z2 )

       + 96*z3*M^2 - 482/3*M^2 - 19508/225*Q.Q^5*z3*M^-8 - 1/462*Q.Q^5*M^-8*z2
          + 129728592122581/1248047088750*Q.Q^5*M^-8 + 42496/675*Q.Q^4*z3*M^-6
          + 4/315*Q.Q^4*M^-6*z2 - 559517166977/7501410000*Q.Q^4*M^-6 - 44*
         Q.Q^3*z3*M^-4 - 3/35*Q.Q^3*M^-4*z2 + 6119442979/125023500*Q.Q^3*M^-4
          + 328/27*Q.Q^2*z3*M^-2 + 4/5*Q.Q^2*M^-2*z2 - 880963/121500*Q.Q^2*
         M^-2 + 11909/486*Q.Q - 392/9*Q.Q*z3 - 67/9*Q.Q*z2;
\end{verbatim}
Note that this expression still contains an overall factor $q^2$ (the
results above are in Euclidean space, indicated by the capital $Q$,
i.e.\ $Q^2 = -q^2$).  After dividing by $q^2$ there are terms of ${\cal
  O}(m^2/q^2)$ which cancel in the sum of all three-loop diagrams.  The
constant terms disappear after requiring the QED-like on-shell condition
$\Pi(0)=0$, and finally the structure of Eq.~(\ref{eqpolfunpi}) is
obtained.

After all diagrams up to three loops are computed and added up (taking
into account the correct colour factors), the parameters $\alpha_s$ and
$m$ are renormalized. To give an impression of the structure of
the final result, we list the first two and the eighth term of the
Abelian part, $\Pi_A^{(2)}$:
\begin{eqnarray}
   C^{(2)}_{A,1} &=&
          - {8687\over 864}
          - {32\over 5}\,\zeta_2\,\ln 2
          + 4\,\zeta_2
          + {22781\over 1728}\,\zeta_3\,,
\nonumber\\
   C^{(2)}_{A,2} &=&
          - {223404289\over 1866240}
          - {192\over 35}\,\zeta_2\,\ln 2
          + {24\over 7}\,\zeta_2
          + {4857587\over 46080}\,\zeta_3\,,
\label{eq::pfcn}
\\
&\cdots&
\nonumber\\
   C^{(2)}_{A,8} &=&
          - {190302182417255312898886115648452691\over 
             63063833203636585050931200000}
\nonumber\\
&&\mbox{}
          - {786432\over 230945}\,\zeta_2\,\ln 2
          + {98304\over 46189}\,\zeta_2
          + {31209476560803609727258477\over 12432176686773043200}\,\zeta_3
\nonumber
\,,
\end{eqnarray}
where $\zeta_n \equiv \zeta(n)$, with Riemann's $\zeta$-function. For
the specific values above it is $\zeta_2 = \pi^2/6$ and $\zeta_3 =
1.20206\ldots$.  The left-out coefficients and the ones for the other
colour structures of (\ref{eqpi2}) can be found in~\cite{CheKueSte96}.

%
\vspace{4ex}
{\bf Large momentum limit}\\[2ex]
%
The calculation for $q^2\gg m^2$ is more involved in the sense that
here the asymptotic expansion does not reduce to a naive Taylor
expansion.  Instead, the large-momentum procedure to be applied in this
case generates a large number of subgraphs for each of the diagrams of
Fig.~\ref{figpol}. Whereas for the small-momentum expansion a computer
only had to be used for the computation of the diagrams, in this
case even the generation of the subdiagrams is a non-trivial task and
program packages like {\tt EXP} or {\tt LMP} (see
Section~\ref{subsubexplmp}) are indispensable.
For the three-loop case the 18 initial diagrams produce 240 subgraphs
and their manual selection would be very tedious and error-prone if not
impossible.  The corresponding numbers in the two- (one-) loop case are
2 (1) initial and 14 (3) subdiagrams which would still allow a manual
treatment.  These numbers also show that there is a big step in
complexity when going from two to three loops.

\begin{figure}[t]
  \begin{center}
  \leavevmode
   \epsfxsize=8em
   \raisebox{-2.7em}{\epsffile[150 260 420 450]{L12c.ps}}\hspace{0em}
   \raisebox{0em}{\Large $\star$}\hspace{1em}
   \epsfxsize=8em
   \raisebox{-2.7em}{\epsffile[150 260 420 450]{LAlmp.ps}}
   \parbox{\captionwidth}{
   \caption[]{\label{figsub}
     Sixth subdiagram as generated by {\tt LMP} for the large-momentum
     procedure of the ladder diagram (cf.\ Fig.~\ref{figladder}).  }}
 \end{center}
\end{figure}


For the calculation the package {\tt LMP} (see
Section~\ref{subsubexplmp}) was used.  As an example let us again
consider the ladder diagram of Fig.~\ref{figladder}.  The input file,
looking similar to the one of the low-energy case above, produces 27
subdiagrams.  The one of Fig.~\ref{figsub} (the sixth subgraph generated
by {\tt LMP}), for example, corresponds to the following code:
\begin{verbatim}
*--#[ d3l3_6 :
        ((-1)
        *Dg(nu1,nu2,+aexp,-p2)
        *Dg(nu3,nu4,-bexp,+p2)
        *S(nu2,+p1mexp,nu4,-bexp,+qmexp,mu2,-p21m,nu3,-p2mexp,nu1,
           -p11m,mu1,-aexp,+qmexp)
        *1);
#define DIANUM "6"
#define TOPOLOGY "arb"
#define INT1 "inpl1"
#define MASS1 "0"
#define INT2 "topL1"
#define MASS2 "M"
#define INT3 "topL1"
#define MASS3 "M"
*--#] d3l3_6 :
\end{verbatim}
After Taylor expansion, the integral factorizes into a massless and two
massive one-loop integrals. The topology of each of these one-loop
integrals is encoded in the pre-processor variables {\tt INT1}, {\tt
  INT2}, and {\tt INT3}, where {\tt inp..} denotes a massless, {\tt
  top..} a massive topology (see also Section~\ref{subsubexplmp}). The
corresponding mass is given by {\tt MASS1}, {\tt MASS2}, and {\tt MASS3}.

Again the first two and the highest available (in this case the seventh)
term in the expansion for $\Pi_A^{(2)}$ will be displayed:
\begin{eqnarray}
   \bar{\Pi}^{(2)}_{A} &=& {3 \over 16 \pi^2}\, \bigg\{
%
       - {143\over 72}
          - {37\over 6}\,\zeta_3
          + 10\,\zeta_5
          + {1\over 8}\,\logqmums
%
\nonumber\\&&\mbox{}
       + {\bar{m}^2\over q^2} \, \bigg[
            {1667\over 24}
          - {5\over 3}\,\zeta_3
          - {70\over 3}\,\zeta_5
          - {51\over 2}\,\logqmums
          + 9\,\logqmums^2
          \bigg] 
\nonumber\\[2ex]&&\mbox{}
+ \,\,\cdots
\nonumber\\[1ex]&&\mbox{}
       + \left({\bar{m}^2\over q^2}\right)^{6} \, \bigg[
          - {420607059143\over 19440000}
          - {13059229\over 2700}\,\zeta_3
          - 1120\,\zeta_4
          + 4000\,\zeta_5
          + {560\over 3}\,B_4
\nonumber\\&&\mbox{\hspace{1.0cm}}
          + \left( {133099291\over 972000}
          + {16842\over 5}\,\zeta_3 \right)\,\logqmms
          + {54076013\over 7200}\,\logqmms^2
          + {8575579\over 2700}\,\logqmms^3
\nonumber\\&&\mbox{\hspace{1.0cm}}
          + \left( {13274779\over 450}
          - {142256\over 15}\,\logqmms
          - 8992\,\logqmms^2 \right)\,\logqmums
\nonumber\\&&\mbox{\hspace{1.0cm}}
          + \left( - 10854
          + 7560\,\logqmms \right)\,\logqmums^2
          \bigg]
%
\bigg\} + \ldots\,,
\label{pi2abar}
\end{eqnarray}
with $\logqmms = \ln(-q^2/\bar m^2)$ and $\logqmums = \ln(-q^2/\mu^2)$,
$\zeta_n$ as in (\ref{eq::pfcn}) and the additional values
$\zeta_4 = \pi^4/90$, $\zeta_5 = 1.03693\ldots$.  $B_4$ was
analytically calculated in \cite{Bro92} and evaluates numerically to
$B_4 = -1.76280\ldots$.  $\bar m$ is the \msbar\ mass which can easily be
transformed to the on-shell scheme~\cite{GraBroGraSch90}. For the
left-out terms in (\ref{pi2abar}) we refer to
\cite{CheHarKueSte96,CheHarKueSte97}.

In Fig.~\ref{figrAv}~(a) the imaginary part of the Abelian contribution
is plotted as a function of $x=2m/\sqrt{s}$ including successively
higher orders in $x$. Up to $x\approx 0.7$ very quick convergence is
observed. For larger values of $x$, however, the inclusion of higher
terms does not significantly improve the approximation.  This becomes
manifest in an even more drastic way in Fig.~\ref{figrAv}~(b), showing
the same curves as functions of $v=\sqrt{1-4m^2/s}$.  This choice of
variable enlarges the threshold region and thus demonstrates the
breakdown of the high-energy expansion close to threshold.
%
\begin{figure}[ht]
\begin{center}
\begin{tabular}{cc}
    \leavevmode
    \epsfxsize=5.5cm
    \epsffile[110 265 465 560]{ravx.ps}
&
    \epsfxsize=5.5cm
    \epsffile[110 265 465 560]{ravv.ps}
\\
$(a)$ & $(b)$
\end{tabular}
\parbox{\captionwidth}{
    \caption[]{\label{figrAv}\sloppy
      The Abelian contribution $R_A^{(2)}$ as functions of $(a)$ $x =
      2m/\sqrt{s}$ and $(b)$ $v=\sqrt{1-4m^2/s}$.  Wide dots: no mass
      terms; dashed lines: including mass terms $(m^2/s)^n$ up to $n=5$;
      solid line: including mass terms up to $(m^2/s)^6$; narrow dots:
      semi-analytical result obtained via Pad\'e
      approximation~\cite{CheKueSte96}.  }}
\end{center}
\end{figure}

%
% end of polfunc.tex
%
