%
\section{Introduction}
%----------------------------------------------------------------------
\subsection{Quantum field theory in particle physics}
%----------------------------------------------------------------------
%
The huge amount of data collected in the recent years mainly at the
Large Electron Positron Collider (LEP) at CERN (Geneva) and the
Stanford Linear Collider (SLC)
as well as at the proton--anti-proton machine TEVATRON 
at FERMILAB (Chicago) has
resulted in an impressive experimental precision for many parameters in
particle physics. For example, the mass and width of the $Z$ boson are
measured with an accuracy which is of the same order of magnitude as
for Fermi's constant~\cite{Kar98,Hol98}.
Therefore it is essential to also improve the
accuracy of predictions following from theoretical models. It is the
purpose of this article to give an overview of some of the technical
tools to perform the corresponding calculations.

At present, the most favoured theoretical model in particle physics is
the so-called Standard Model. It incorporates all known quarks and
leptons and describes their interactions by the gauge group
SU(3)$\times$SU(2)$\times$U(1). Here, the SU(3) generates the strong
interactions among the quarks by the exchange of eight massless gluons,
while the SU(2)$\times$U(1), usually called the electroweak sector,
comprises electromagnetic and weak interactions, carried by the
massless photon and the massive charged $W^\pm$ and neutral $Z$ boson,
respectively.

Up to now all theoretical predictions deduced from the Standard Model
have found perfect confirmation in experiment. For example, the
existence of the top quark was predicted after the discovery
of the $\Upsilon$-resonances in 1977 \cite{Herb77} merely out of the
need for an SU(2) doublet partner for the bottom quark. At least equally
impressive, however, was the accurate prediction of its mass, $M_t$, by
comparing precision measurements of certain observables at LEP and SLC
to their theoretical values. The top quark mass dependence in these
observables only appears as a higher order effect through the so-called
radiative corrections, typical for quantum field theories. So when the
first evidence of the discovery of the top quark was announced by the
CDF and D0 collaborations~\cite{TOP95}, the agreement of directly and
indirectly measured values for the top mass served as a very convincing
argument.  Today statistics have improved considerably and the original
signal has been impressively confirmed~\cite{CDFD098}.

A similar situation exists with the Higgs boson, at the present time
being the only particle of the Standard Model not yet discovered. Its
existence, however, is by no means guaranteed as was the top quark after
$b$ discovery.  Furthermore, predictions for the Higgs mass $M_H$ suffer
from much larger uncertainties because the dependence of the relevant
radiative corrections on $M_H$ is only logarithmic~\cite{Vel77} while it
is quadratic on $M_t$.  Nevertheless, it is possible to set the limits
to $M_H=76^{+85}_{-47}$~GeV with an upper bound of $M_H< 262$~GeV at
$95\%$ confidence level~\cite{Teu98}.  The current lower limit for $M_H$
from the direct search at LEP reads $M_H\agt 89.8$~GeV~\cite{Mh_LEP}.

From the theoretical point of view the Standard Model does not look like
a ``final'' theory. First of all, it contains a lot of free parameters
and any attempt to predict their values needs a generalization of the
model.  Furthermore, there is no way to find unification of the strong,
weak and electromagnetic coupling constant within the Standard Model.  A
very promising and theoretically attractive extension is the Minimal
Supersymmetric Standard Model (MSSM) which not only allows for gauge
unification but even solves the problem of quadratic divergences and
many others. Phenomenological calculations, however, are even more
difficult here since the particle content is approximately doubled.
Moreover, some important technical tools are no longer applicable.

Let us now turn to the question how to extract physical information like
cross sections and decay rates from the Standard Model, leading to
predictions like the values for $M_t$ and $M_H$ mentioned above. At the
moment the vast majority of the calculations are based on perturbation
theory, represented by the extremely depictive Feynman diagrams.
Together with the Feynman rules, which constitute the translation rules
from the graphical to a mathematical notation, they are in complete
one-to-one correspondence to the individual terms of the perturbative
series. As long as one is interested only in perturbative results, one
may formulate the theory in terms of the Feynman rules alone, i.e., by
providing the list of possible vertices, describing particle
interactions, and propagators. Any physical process may then be
evaluated by connecting the vertices and propagators in all possible
ways to give the desired initial and final states.  The specific order
in perturbation theory is reflected in the number of vertices of the
corresponding Feynman diagram. Given a fixed initial state and
increasing the number of vertices leads to the appearance either of
closed loops or of additional particles in the final state. In the
former case, this implies an integration over the corresponding loop
momentum, while the latter case results in a more complicated phase
space integral when calculating the rate. Usually, diagrams with more
than one closed loop are called multi-loop, those with more than two
particles (``legs'') in the final state are called multi-leg diagrams.
The technology, in particular the automation of the computations in the
multi-loop sector is much further developed than for multi-leg diagrams.
Furthermore, to a certain extent it is possible to relate multi-leg
diagrams to multi-loop ones via the optical theorem. Thus our main focus
in this article lies in the calculation of loop diagrams.

%----------------------------------------------------------------------
\subsection{Motivation for automatic computation of Feynman diagrams}
%----------------------------------------------------------------------
%
The first realization of the idea to pass purely algebraic or, more
precisely, symbolic operations to a computer indeed was driven by
particle physics, when in 1967 M.~Veltman developed the program {\tt
  SCHOONSCHIP} \cite{VelSS}, mainly to control the evaluation of fermion
traces, i.e.\ traces of gamma matrices. Today, {\tt SCHOONSCHIP} to a
large extent is superceded by other systems which have, however, taken
over many of its basic ideas.  Some important examples from the early
days that also have their roots in particle physics are {\tt
  MACSYMA}~\cite{macsyma}, {\tt REDUCE}~\cite{reduce} and others. The by
now highly developed algebraic system {\tt Mathematica}~\cite{Wolfram}
actually is a derivative of the program {\tt SMP}, created in 1980 by
the particle physicist S.~Wolfram.  While some of these programs claim
to be fairly general, the most direct descendant of {\tt SCHOONSCHIP} is
{\tt FORM}~\cite{form}, a system mainly tailored to high energy physics.
In Section~\ref{secalgprg} a few representatives of the most important
algebraic systems will be introduced.  The main concern of this review
are not these algebraic programs, but rather their applications, i.e.,
packages based on these systems.  It should be clear, however, that they
are very important prerequisites for the automation of Feynman diagram
calculations.

The higher the order of perturbation theory under consideration, the
larger is the number of contributing Feynman diagrams and the more
complicated is their individual evaluation. This makes it desirable and
even unavoidable to develop efficient algorithms which can be
implemented on a computer in a simple way.  Therefore, large efforts in
the field of Feynman diagram evaluation are devoted to find algorithms
for the specific types of operations arising in a typical calculation.
One may characterize these operations by analyzing the steps necessary
for a perturbative field theoretic calculation: given a set of Feynman
rules, first the contributing Feynman diagrams for the process under
consideration have to be generated.  This, being mainly a question of
combinatorics, is clearly a well suited problem for computerization, and
meanwhile very effective algorithms have been implemented with emphasis
on slightly diverse purposes. The next step is the evaluation of the
corresponding mathematical expressions for the diagrams which can be
divided into an algebraic and an analytic part. The algebraic part
consists of operations like contracting Lorentz indices and calculating
fermion traces.  Nowadays many algebraic systems provide functions
optimized for this kind of manipulations.

The analytic part of a Feynman diagram calculation is concerned with the
Feynman integrals, and usually this is the point where computers fail to
be applicable without human intervention.  However, Feynman integrals
establish a very special class of integrals, and for huge subclasses
algorithms exist that {\em algebraically} reduce each of them to basic
sets of meanwhile tabulated expressions.  Two of the most important
techniques of this kind are the so-called tensor
reduction-~\cite{PasVel79} and the integration-by-parts
algorithm~\cite{CheTka81}.  The first one reduces any Lorentz structure
in the integrand to invariants and has been worked out explicitly at the
one-loop level for the general case and up to two loops for propagator
diagrams. Using the integration-by-parts algorithm, on the other hand,
recurrence relations can be derived which express diagrams with
non-trivial topologies through simpler ones at the cost of increasing
their number and the degree of the denominators.  The practical
realization of these algorithms would be unfeasible in realistic
calculations without the use of powerful computer systems.

From the considerations of the previous paragraphs one can see that as
long as one is interested in the {\it exact} evaluation of Feynman
diagrams, the number of problems that can completely be passed to
computers is rather limited which is mainly due to the relatively small
number of exactly solvable Feynman integrals. There are, however, many
processes where some kind of approximate result may be equally helpful.
This becomes clear by recalling that, on the one hand, working at fixed
order in perturbation theory is an approximation anyway, and, on the
other hand, experimental results are not free of uncertainties
either.
Therefore, it is enormously important to have approximation procedures
for Feynman diagrams. For instance, the numerical evaluation of
complicated Feynman integrals, which is a typical task for a computer,
is a possibility to obtain results with finite accuracy.
This approach is currently favoured for the automation of tree- and one-loop
calculations involving lots of different mass parameters like in the
Standard Model or even the MSSM. For recent developments in the
numerical evaluation of Feynman diagrams beyond the one-loop level
we refer to~\cite{numcalc}.

A different kind of approximation may be obtained by expanding the
integrals with respect to ratios of different scale parameters. Here the
so-called asymptotic expansions of Feynman diagrams are becoming more
and more popular, since they allow any multi-scale diagram to be reduced
to single-scale ones, at least in principle. While the evaluation of the
former may be a hard job even at one-loop level, the latter ones can be
computed up to three loops using the integration-by-parts algorithm
mentioned before.  The required efforts for a manual application of
asymptotic expansions increase steeply with the number of loops.
Meanwhile, however, the most important variants have been implemented
and have also been applied to the calculation of three-loop
contributions to several important physical observables.

Performing huge calculations with computers inevitably leads to seemingly
unrelated difficulties of purely organisational nature: Computers
nowadays are not yet as stable as one would like them to be. In this way
they force the user to cover on possible breakdowns, keeping the loss
minimal. In the case where the long run-time is mainly due to a huge
number of diagrams, each of them taking only a few minutes to be
evaluated, the solution to this problem is rather straightforward.
Saving each diagram on disk after computation naturally guarantees that
the loss will be of only a few minutes (apart from more severe problems
like disk crashes or the like).  For more complicated diagrams, however,
it is helpful if the algebra program itself is structured in a way that
intermediate results will not be lost after system crashes.

All the algorithms mentioned above allow the calculation of Feynman
diagrams to be automated to more or less high degree, where one will
certainly choose different combinations of them for different purposes.
It is quite difficult to give an exhaustive definition of automation.
For instance, often the computer system either applies only for one
special problem (or a very restricted class of problems), nevertheless
saving the user a great deal of effort and preventing otherwise unavoidable
errors. In other cases the system is rather flexible but requires a
certain amount of human intervention for each specific problem.  The
final goal of automation where a typical screen snapshot would look like
\begin{verbatim}
  model: SM 
  loops: 3 
  initial-state: ep,em 
  final-state: t,tb 
  parameters: default 
  observable: sigma_tot 
  energy: 500 GeV

    sigma = 0.8 pb

\end{verbatim}
is at the moment still out of range.
Nevertheless, many of the recently performed calculations in high energy
physics would never have been possible without a high degree of
automation. Let us, at the end of this introduction, list a few of these
applications that will be addressed in more detail in
Section~\ref{secapplications}.

%----------------------------------------------------------------------
\subsection{Examples for physical applications}
%----------------------------------------------------------------------
%
Most of the multi-loop and multi-leg computations have been performed
for the QCD sector of the Standard Model. On the one hand they are
numerically most important simply because the coupling constant is
relatively large. On the other hand they are simpler to evaluate due to
the fact that very often only one dimensionful scale appears like the
external momentum or a large internal mass.

The ``classical'' example which is currently known to order $\alpha_s^3$
are QCD corrections to the hadronic R ratio,
$R(s)=\sigma(e^+e^-\to\mbox{hadrons})/\sigma(e^+e^-\to\mu^+\mu^-)$.
Already in 1979 ${\cal O}(\alpha_s^2)$ corrections have been
computed~\cite{CheKatTka79DinSap79CelGon80} in the massless limit.  A
large part of the calculations was carried out without using computers
and lots of efforts had to be spent on nowadays trivial things like
fermion traces.  It took more than ten years until the next order was
completed~\cite{GorKatLar91SurSam91,Che97R}.  Even then the manipulation
of the diagrams was largely done ``by hand''. As a consequence,
important checks like gauge parameter independence could not be
performed.  Only recently the ${\cal O}(\alpha_s^3)$ result was checked
by a completely independent calculation which used the powerful methods
of automatic generation and computation of diagrams~\cite{Che97R}.

The full mass dependence of the ${\cal O}(\alpha_s^2)$ corrections to
$R(s)$ has become available~\cite{CheKueSte96} only very recently, when
high moments of the polarization function were evaluated.  In this
calculation the problem was not in the number of diagrams which is less
than 20, but in the huge expressions occurring in intermediate steps.  (For
some diagrams up to four giga-bytes of disk space were required.)  For
such computations it is indispensible that enormously powerful algebra
systems are available and that their use is completely automated to
avoid any possible human error in intermediate steps.

The number of diagrams to be considered for the ${\cal O}(\alpha_s^3)$
corrections to $R(s)$ at $m=0$ is of the order of one hundred and thus
it is still possible to organize the calculation by hand.  The situation
is different for the decay of the Higgs boson into gluons, for example.
Among other reasons this is a very important process because the inverse
reaction of gluon fusion is the dominant production mechanism of Higgs
bosons at the future Large Hadron Collider (LHC). Since the leading
order QCD corrections to the loop-induced process $H\to gg$ amount to
roughly $68\%$~\cite{InaKubOka83,DjoSpiZer91}, it was necessary to
evaluate the next-to-leading order contribution. However, the number of
three-loop diagrams to be computed in this case is of the order one
thousand.  This in combination with the large expressions in the
intermediate calculational steps is the reason why one has to rely on a
high degree of automation as will be described in
Section~\ref{sec:applic:hgg}.

An even larger number of diagrams had to be dealt with when the
renormalization group functions $\beta$ and $\gamma_m$, governing the
running of the strong coupling constant and the quark masses,
respectively, were computed at four-loop level.  Roughly $50\,000$
diagrams for $\beta$~\cite{RitVerLar97}
and $2\,000$ for 
$\gamma_m$~\cite{Che97VerLarRit97}
contributed, and it is
quite clear that such calculations would never have been possible
without the intensive use of extremely efficient computer algebra
systems in combination with administrative software concerned with
book-keeping.

So far, none of the examples above really had to rely on one of the
approximation procedures like numerical integration or asymptotic
expansion mentioned above. This becomes relevant for the decay of the
$Z$ boson into quarks at ${\cal O}(\alpha\alpha_s)$, for example. Here,
the bottom quark channel is of special interest since the top quark
appears as a virtual particle.  With the present technology, a
calculation of the full $M_t$-dependence is out of reach because up to
four different mass scales are involved.  However, asymptotic expansions
provide a very promising method to get a result which is almost
equivalent to the full answer. In the approach which was used to tackle
the ${\cal O}(\alpha\alpha_s)$ terms to this process~\cite{HarSeiSte97},
the 69 contributing diagrams are split into 234 subdiagrams.  The manual
application of the method of asymptotic expansions was therefore
impossible and the existence of program packages that apply this
procedure automatically was important.

Meanwhile, experimental accuracy allows the electroweak sector of the
Standard Model to be tested even at the two-loop level, and several
groups have tackled the corresponding calculations.  As was already
outlined, the strategy here is to reduce all integrals to scalar ones.
The latter depend on a lot of different mass scales, such that their
evaluation must be done numerically.  In Section~\ref{sec:applic:deltar}
we will describe the computation~\cite{BauWei97} of two-loop corrections
induced by the Higgs boson to $\Delta r$ which enters the relation
between $G_F, M_Z, M_W$ and $\alpha$.

It is not only since the increase of the center-of-mass energy at LEP
above the $W$ boson production threshold that processes involving four
(or more) particles in the final state became very topical.  In such
reactions even the contribution from lowest order perturbation theory may
pose a serious problem as quite a lot of diagrams are involved and a
highly non-trivial phase space integration has to be performed. In
Section~\ref{sec:applic:strongW} the scattering of vector bosons 
is discussed in the background of a potentially large coupling among the
bosons at high energies~\cite{BHKPYZ98}.

Other applications that will not be discussed in more detail but shall
further substantiate the success of computer algebra in Feynman diagram
calculations are, for example, QCD corrections to
$\Delta\rho$~\cite{Avd95,CheKueSte951} and $\Delta
r$~\cite{CheKueSte952}, moments of QCD structure
functions~\cite{LarTkaVer91,LarRitVer94,LarNogRitVer96}, top mass
effects in the decay of intermediate~\cite{LarRitVer95CheKwi96} and
heavy~\cite{HarSte97} Higgs bosons as well as in $e^+e^-$ collisions
\cite{HarSte98}.  Involved one-loop calculations concerning radiative
corrections to the gauge boson scattering were considered
in~\cite{DenDitHah97}.  The muon anomalous magnetic
moment~\cite{CzaKraMar96} and the neutron anomalous electric
moment~\cite{CzaKra97} are important examples of two-loop calculations
in the electroweak sector.  In~\cite{FriKniKreRie96} the main difficulty
was the evaluation of two-loop vertex diagrams in order to obtain
corrections to the decay of a heavy Higgs boson.  Investigations of $b$
decay~\cite{CzaMel97,CzaMel98b} and threshold production of heavy
quarks~\cite{CzaMel98,BenSmi98,HoaTeu98,BenSigSmi98} are examples for
the developments of new algorithms and their implementation.  A large
list of results obtained with the help of programs that automatically
compute tree-level processes involving many particles can be found
in~\cite{comphep}.

Many more important calculations for physical observables in the
multi-loop and multi-leg sector have been performed --- mainly within
the last few years, but let us close the list at this point.

The outline of this review is as follows: Section~\ref{sec:tools}
introduces the most important technical tools that allow the automation
of Feynman diagram calculations. This includes on the one hand some of
the required algorithms, and on the other hand the most frequently used
algebraic programming languages to implement these algorithms.  Concrete
examples for such implementations will be described in
Section~\ref{secimpl}. This splits into a surveying part whose main
purpose is to collect a list of the most important programs, and a more
specific part concerned with the actual realization in the light of a
representative set of selected packages.  Finally,
Section~\ref{secapplications} gives a list of physical applications that
have been considered in the literature and that would not have been
feasible without the help of a certain degree of automation. Some of
these examples actually were the driving force for some of the
afore-mentioned packages to be developed. Others are real applications
in the sense that existing programs could successfully be used to tackle
as yet unsolved problems even by people that were not involved in the
development of the programs.

%
% end of introduction.tex
%
