%\documentclass[11pt]{cernrep}
%\usepackage{graphicx}
%\usepackage{here}
%\usepackage{cite}

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

\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}

\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\nn}{\nonumber}

\def\eqn#1{Eq.~(\ref{#1})}
\def\eqns#1#2{Eqs.~(\ref{#1}) and~(\ref{#2})}
\def\eqnss#1#2{Eqs.~(\ref{#1})-(\ref{#2})}
\def\fig#1{Fig.~{\ref{#1}}}
\def\sec#1{Section~{\ref{#1}}}
\def\app#1{Appendix~\ref{#1}}
\def\tab#1{Table~\ref{#1}}

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


\def\ket#1{|{#1}\rangle}
\def\bra#1{\langle{#1}|}
\def\braket#1#2{\langle #1 |#2 \rangle}
\def\ord{{\cal O} }
\def\e{\epsilon}
\def\eps{\epsilon}
\def\cm{{\cal M}}
\def\cM{{\cal M}}
\def\cmb{\overline{{\cal M}}}
\def\cmqq{{\cal M}_{q\bar q}}
\def\cmbqq{\overline{{\cal M}}_{q\bar q}}
\def\cmqqg{{\cal M}_{q\bar qg}}
\def\cmbqqg{\overline{{\cal M}}_{q\bar qg}}
\def\cmqqgg{ {\cal M}_{q\bar qgg}}
\def\cmqqqq{ {\cal M}_{q\bar qq\bar q}}
\def\asopi{\left(\frac{\alpha_s}{2\pi}\right)}
\def\asmuopi{\left(\frac{\alpha_s(\mu^2)}{2\pi}\right)}
\newcommand\sss{\scriptscriptstyle}
\newcommand\als{\alpha_{\sss S}} 
\newcommand\gs{g_{\sss S}}
\def\tgs{\tilde\gs}
\def\CF{C_{\sss F}}
\def\CA{C_{\sss A}}
\def\NF{N_{\sss F}}
\def\TF{T_{\sss R}}
\def\ca#1{{\cal A}^{( #1)}}
\def\caqq#1{{\cal A}^{( #1)}_{q\bar q}}
\def\caqqg#1{{\cal A}^{( #1)}_{q\bar qg}}
\def\caqqgg#1{{\cal A}^{( #1)}_{q\bar qgg}}
\def\caqqqq#1{{\cal A}^{( #1)}_{q\bar qq\bar q}}
\def\sab{s_{12}}
\def\sac{s_{13}}
\def\sbc{s_{23}}
\def\yab{y_{12}}
\def\yac{y_{13}}
\def\ybc{y_{23}}
\def\sij{s_{ij}}
\def\sabc{s_{123}}
\def\Ione{{\mbox{\boldmath $I$}^{(1)}(\e)}}
\def\dIone{{\mbox{\boldmath $\Delta I$}^{(1)}(\e)}}
\def\Ionep{{\mbox{\boldmath $I$}^{(1)}(2\e)}}
\def\Htwo{{\mbox{\boldmath $H$}^{(2)}(\e)}}
\def\MSbar{$\overline{{\rm MS}}$}


%\begin{document}
%\bibliographystyle{unsrt}

%\title{THE REPORT OF THE WORKING GROUP ON HIGHER ORDERS}
%\author{T. Binoth$^1$,   V. Del Duca$^2$, 
%T. Gehrmann$^3$, A. Gehrmann-De Ridder$^3$,
%E.W.N. Glover$^4$, J.-Ph. Guillet$^5$ and G. Heinrich$^4$}
%\institute{$^1$Department of Physics and Astronomy, 
%           The University of Edinburgh,
%	   EH9 3JZ Edinburgh, Scotland\\
%$^2$Istituto Nazionale di Fisica Nucleare, Sez. di Torino,
%via P. Giuria, 1 - 10125 Torino, Italy\\  
%$^3$TH Division, CERN, CH-1211, Geneve 23,  Switzerland\\ 
%$^4$Institute for Particle Physics Phenomenology, 
%University of Durham, Durham, DH1 3LE, U.K.\\
%$^5$Laboratoire d'Annecy-Le-Vieux de Physique Theorique LAPTH, Chemin de
%Bellevue, B.P. 110, F-74941 Annecy-Le-Vieux, France}

\section{HIGHER ORDERS\protect\footnote{Section coordinator:
    E.W.N.~Glover}$^{,}$~\protect\footnote{Contributing authors: 
    T.~Binoth, V.~Del~Duca, T.~Gehrmann, A.~Gehrmann-De~Ridder, 
    E.W.N.~Glover, J.-Ph.~Guillet and G.~Heinrich}}

%\maketitle
%\begin{abstract}
In this report, we summarize the issues discussed and worked out during
the 'Working Group on Higher Orders'. The two current frontiers of higher
order QCD calculations at colliders are the description of multiparton
final states at next-to-leading order and the extension of calculations
for precision observables beyond this order.  Considerable progress has
been made on both issues, with the highlights being the first calculation
of the Yukawa-model one-loop six-point amplitude and the two-loop
corrections to the $e^+e^- \to 3$~jets matrix element. Developments
towards the construction of NNLO parton level Monte Carlo programs are
described for the example of $\gamma^* \to q\bar q$. Finally, the first
applications of two-loop matrix elements to the improved description
of the high energy limit of QCD are reported.
%\end{abstract}

%\tableofcontents

\subsection{Introduction}
Present and future collider experiments will confront  us with large sets of
data on multi--particle final states. This is particularly true for the hadron
colliders Tevatron  and LHC, which will be the machines operating at the 
highest attainable energies in the near future.  Hence the comparison of jet
observables to theoretical  predictions will become increasingly important. 

Since the theoretical predictions based on leading order (LO) calculations  
are typically plagued by large scale uncertainties,  it becomes necessary to
calculate the next-to-leading order (NLO) corrections  in order to make
meaningful predictions which match the experimental precision.  Indeed, the
$N$-jet cross section  in hadronic collisions  is proportional to $\alpha_s^N$
at leading order, which means that  theoretical uncertainties  are actually
amplified for growing $N$. 

For processes with relatively few jets, such as dijet production,
next-to-next-to-leading order (NNLO) perturbative predictions will be needed to
reduce the theoretical uncertainties and enable useful physics to be extracted
from the copious high-precision data.

In this report, we address issues connected with theoretical progress in
calculating higher order corrections both at NLO (in the context of $2 \to 4$
scattering processes) and at NNLO (in the context of $2\to 2$ processes).   At
present, making numerical predictions for these types of processes lies
well beyond our capabilities.   However, there has been very rapid progress in
the last two years and it is very likely that the technical stumbling blocks
will be removed.   

Our report is structured as follows.   First,  we address the issue of NLO
corrections to multi-particle final states. The main technical problems
associated with dimensionally regulated pentagon integrals were solved some time
ago~\cite{Bern:1994kr} and the next-to-leading order matrix
elements for $2\rightarrow 3$ processes have become available in the recent 
years\,\cite{Bern:1993em,Bern:1994kr,Kunszt:1994tq,Signer:1995nk,
Bern:1997ka,Bern:1998sd,Campbell:1997tv,Glover:1997eh,DelDuca:1999pa,
Beenakker:2001rj,Reina:2001bc,DelDuca:2001fn}.
However, the step to
$2\rightarrow 4$ or even higher processes at NLO has not been made yet. The
reason lies in the fact that the computation of the  corresponding amplitudes
is highly nontrivial.   Although the calculation techniques  for amplitudes
with an arbitrary number of external  legs  are
available\,\cite{Binoth:1999sp}, it turns out that a brute force approach is
not viable. In order to  avoid intractably large expressions in the calculation
of  six-point (or higher) amplitudes, it is indispensable to understand better
recombination and cancellation mechanisms at intermediate steps of the
calculation.   This issue is addressed in Sec.~\ref{sec:ho;nlo}, in the
context of the Yukawa model, where all external legs are  massless scalars
attached to a massless fermion loop.

In Sec.~\ref{sec:ho;nnlo} we consider the NNLO corrections to $2\to 2$
scattering processes.   The rationale for going beyond the next-to-leading
order is reviewed in Sec.~\ref{sec:ho;nnlorationale} while the  various
building blocks necessary for such a calculation are discussed in
Secs.~\ref{sec:ho;nnlopdf}-- \ref{sec:ho;two-loop}.   While the individual
components are in relatively good shape - the infrared limits are well studied,
many two-loop matrix elements exist and the NNLO evolution of parton
distributions is almost under control -- a systematic procedure for combining
them to give numerical predictions is not established.   Therefore in
Sec.~\ref{sec:ho;exercise} we examine the infrared singular structure of the
various pieces for $\gamma^* \to 2$~jets which is one of the simplest
non-trivial processes at NNLO.

Sec.~\ref{sec:ho;highenergy} contains a summary of the present status of the
analytic structure of QCD amplitudes in the limit of forward and backward
scattering. In these high energy limits, the scattering process  is  dominated
by the exchange of a particle in the $t$- (or $u$-)channel respectively. This
may be the gluon or the quark.   In both cases, the amplitude reggeises and the
large logarithms can be resummed to next-to-leading-logarithmic accuracy by
simple forms involving the reggeised particle exchange together with modified
vertex functions (or impact factors).  We evaluate the gluon and quark Regge
trajectories to two-loop accuracy and show that they are strikingly similar:
the gluon Regge trajectory can be obtained from the quark trajectory by mapping
$\CF \to \CA$.
 
A brief summary and outlook is given in Sec.~\ref{sec:ho;conc}.


\subsection{NLO}
\label{sec:ho;nlo}

For NLO amplitudes with $N\ge 6$ external particles,
standard calculational methods are not adequate because of
the complexity of intermediate expressions.

In the last few years, methods either directly based on string perturbation 
theory\,\cite{Bern:1992aq,Bern:1991ux,Bern:1993wt}  or on a world line
formulation of field theory  amplitudes\,\cite{Strassler:1992zr} have been used
to derive a number of ``master formulae'' for one-loop $N$--point amplitudes.
Those are generating functionals which  yield, for any $N$,  a closed parameter
integral expression for the amplitude. The resulting integral representations
are related to standard Feynman parameter integrals in a well-understood 
way\,\cite{Bern:1992an}. Nevertheless,  due to their superior organisation they
often allow one to exploit at the integral level properties of an amplitude
which normally would be seen only at later stages in a Feynman graph
calculation\,\cite{Schubert:2001he,Frizzo:2000ez}.

Although the string inspired formalism allows for an elegant formulation of
amplitudes in terms of  a manifest Lorentz structure, one is in general not at
all dispensed from  doing cumbersome algebraic work.  The complexity of doing
tensor reduction in momentum space translates into the need to reduce Feynman
parameter integrals with nontrivial numerators to genuine $N$--point scalar
integrals. These can be expressed in terms of box, triangle and bubble scalar
integrals.   Substantial cancellations appear in all these steps and progress
in finding efficient calculational methods relies on a better understanding of
these mechanisms.

Here we will sketch the calculation of the six-point one-loop  amplitude in the
Yukawa model where all external legs are  massless scalars attached to a
massless fermion loop~\cite{Binoth:2001vm}.  The interest of this model is related to the fact that 
the appearance of tensor integrals can be completely avoided. This can
immediately be seen from the string inspired master formulas for one-loop
$N$--point functions derived in\,\cite{Schubert:2001he}. In this way one can
study the reduction mechanisms for scalar  6-point functions without the
additional complications arising from the tensor reduction.  Having understood
these mechanisms one can proceed towards the computation of gauge theory
amplitudes.



\subsubsection{Calculation of a hexagon amplitude in the Feynman diagrammatic
approach}
 
The amplitude $\Gamma_{\rm yuk}$ can be written as a sum over $6!$ permutations
of the external momentum vectors $p_1,\dots,p_6$,
\begin{equation}\label{Eq:start}
\Gamma^{\phi}_{\rm yuk}[p_1,p_2,p_3,p_4,p_5,p_6] = -\frac{g^6}{(4\pi)^{n/2}}
\frac{1}{6} 
 \sum\limits_{\pi\in S_6} {\cal A}(p_{\pi_1},p_{\pi_2},
                  p_{\pi_3},p_{\pi_4},p_{\pi_5},p_{\pi_6}).
\end{equation} 
Each permutation corresponds to a single Feynman diagram. The amplitude 
for the trivial permutation is given by,
\begin{eqnarray}\label{amp_graph}
{\cal A}(p_1,p_2,p_3,p_4,p_5,p_6) &=& \int \frac{d^nk}{i \pi^{n/2}}
\frac{{\rm tr}(q_1,q_2,q_3,q_4,q_5,q_6)}{q_1^2q_2^2q_3^2q_4^2q_5^2q_6^2},
\end{eqnarray}
where $q_j = k - r_j = k-p_1-\dots-p_j $.
%Note that all indices labelling momenta are
%understood to be mod 6 throughout this paper. 
Working out the trace gives a sum of products of terms 
$q_k\cdot q_j$ which can be written as 
inverse propagators which cancel directly. 
This means that each graph can simply be represented as
a linear combination of scalar integrals. 
For the trivial permutation we find
\begin{eqnarray}\label{amp_si}
2 \,{\cal A}(p_1,p_2,p_3,p_4,p_5,p_6) =
  4\, I_3^n(p_{12},p_{34},p_{56}) 
 +4\, I_3^n(p_{23},p_{45},p_{61})  \nonumber\\
 +{\rm tr}(p_1,p_2)\, I_4^n(p_1,p_2,p_{34},p_{56}) 
 +{\rm tr}(p_2,p_3)\, I_4^n(p_2,p_3,p_{45},p_{61}) \nonumber\\
 +{\rm tr}(p_3,p_4)\, I_4^n(p_3,p_4,p_{56},p_{12})  
 +{\rm tr}(p_4,p_5)\, I_4^n(p_4,p_5,p_{61},p_{23}) \nonumber\\
 +{\rm tr}(p_5,p_6)\, I_4^n(p_5,p_6,p_{12},p_{34}) 
 +{\rm tr}(p_6,p_1)\, I_4^n(p_6,p_1,p_{23},p_{45})  \nonumber\\
 +{\rm tr}(p_1,p_4)\, I_4^n(p_1,p_{23},p_4,p_{56}) 
 +{\rm tr}(p_2,p_5)\, I_4^n(p_2,p_{34},p_5,p_{61}) \nonumber\\
 +{\rm tr}(p_3,p_6)\, I_4^n(p_3,p_{45},p_6,p_{12})  \nonumber\\
 +{\rm tr}(p_1,p_2,p_3,p_4)\, I_5^n(p_{56},p_1,p_2,p_3,p_4)  
 +{\rm tr}(p_2,p_3,p_4,p_5)\, I_5^n(p_{61},p_2,p_3,p_4,p_5)  \nonumber\\
 +{\rm tr}(p_3,p_4,p_5,p_6)\, I_5^n(p_{12},p_3,p_4,p_5,p_6)  
 +{\rm tr}(p_4,p_5,p_6,p_1)\, I_5^n(p_{23},p_4,p_5,p_6,p_1)  \nonumber\\
 +{\rm tr}(p_5,p_6,p_1,p_2)\, I_5^n(p_{34},p_5,p_6,p_1,p_2)  
 +{\rm tr}(p_6,p_1,p_2,p_3)\, I_5^n(p_{45},p_6,p_1,p_2,p_3)  \nonumber\\
 +{\rm tr}(p_1,p_2,p_3,p_4,p_5,p_6)\, I_6^n(p_1,p_2,p_3,p_4,p_5,p_6) .
\end{eqnarray}
The arguments of the $N$--point scalar integrals are the momenta
of the external legs. We use the abbreviation 
$p_{ijk\dots}=p_i+p_j+p_k+\dots$.  
The spinor traces can be expressed by Mandelstam variables 
defined by the 9 cuts of the hexagon graph, but the form given above
is not only most compact but also most convenient to proceed.
%with $s_{l,l+1} = ( p_l + p_{l+1})^2$, 
%$s_{l-1,l,l+1} = ( p_{l-1} + p_l + p_{l+1})^2$. 
%The remaining traces are defined by cyclic relabelling.
%In the following the momenta inside the traces will be represented
%by their indices only.

We also note
that the amplitude is free of infrared poles. This can be seen
by power counting for the soft and collinear poles. 

To sketch the explicit calculation of the hexagon amplitude, 
we will draw special attention to the cancellation
mechanisms at work. First one has to
reduce hexagon and pentagon integrals to box 
integrals.
Then the explicit expressions for the box integrals are inserted. 
Finally the coefficients of  different terms are combined and simplified
by using linear relations for the reduction coefficients. We note
that in none of these steps the size of the expression 
will blow up, while this would surely be the case in a brute force approach.

The reduction formula for the hexagon integral reads,
\begin{eqnarray}
\label{hexa_red}
I_6^n(p_1,p_2,p_3,p_4,p_5,p_6) 
=\frac{1}{\det(\hat S)}
\bigl\{\hspace{5cm}\nonumber \\
\bigl[ {\rm tr}(123456){\rm tr}(3456)-2 s_{34}s_{45}s_{56} {\rm tr}(6123) \bigr]\; I_5^n(p_{12},p_3,p_4,p_5,p_6)\nonumber \\
+\bigl[{\rm tr}(123456){\rm tr}(4561)-2 s_{45}s_{56}s_{61} {\rm tr}(1234) \bigr]\; I_5^n(p_{23},p_4,p_5,p_6,p_1)\nonumber\\
+\bigl[{\rm tr}(123456){\rm tr}(5612)-2 s_{56}s_{61}s_{12} {\rm tr}(2345) \bigr] \;I_5^n(p_{34},p_5,p_6,p_1,p_2)\nonumber\\
+\bigl[{\rm tr}(123456){\rm tr}(6123)-2 s_{61}s_{12}s_{23} {\rm tr}(3456) \bigr]\; I_5^n(p_{45},p_6,p_1,p_2,p_3)\nonumber\\
+\bigl[{\rm tr}(123456){\rm tr}(1234)-2 s_{12}s_{23}s_{34} {\rm tr}(4561) \bigr]\;I_5^n(p_{56},p_1,p_2,p_3,p_4)\nonumber\\
+\bigl[{\rm tr}(123456){\rm tr}(2345)-2 s_{23}s_{34}s_{45} {\rm tr}(5612)
\bigr]\; I_5^n(p_{61},p_2,p_3,p_4,p_5) \bigr\},  
\end{eqnarray}
where the momenta inside the traces are represented
by their indices only. The coefficients in front of
the 5--point integrals are called $b_j$ ($j \in \{1,\dots ,6\}$) 
in the following.
They are defined by the linear equation,
\begin{eqnarray}\label{hexa_coeffs}
( \hat S \cdot b)_j = 1  &\Leftrightarrow& b_j = 
\sum\limits_{k=1}^6 \hat S^{-1}_{kj}\;\mbox{ where }\;\hat S_{kj}=(r_k-r_j)^2\\
\det(\hat S)&=& 4 s_{12}s_{23}s_{34}s_{45}s_{56}s_{61} 
- {\rm tr}(123456)^2.\nonumber
\end{eqnarray}
The traces allow for a compact notation for the coefficients $b_j$.
The  Gram matrix
$G_{kl} = 2 \,r_l \cdot r_k$ is related to 
$\hat S$ by $\hat S_{kl}=-G_{kl}+r_k^2+r_l^2$. 
For $N\ge 6$ and 4-dimensional external momenta one has det$(G)=0$, 
which leads to a non-linear constraint
between the Mandelstam variables. We note that this
constraint is represented {\em linearly} in terms of the
coefficients $b_j$. One has,
\begin{eqnarray}\label{hexa_cons}
\det(G) = 0   \hspace{1cm} \Leftrightarrow  
\hspace{1cm} \sum\limits_{j=1}^6 b_j = 0.
\end{eqnarray} 
By solving eq.~(\ref{hexa_coeffs}) with Cramer's rule
one sees that the constraint (\ref{hexa_cons}) relates
sums of determinants of 5$\times$5 matrices. In
terms of Mandelstam variables these are huge expressions 
just representing zero.
The guideline to keep the sizes of expressions under control
in calculations of multi--particle processes
is thus to use representations of amplitudes where the $b_j$
are kept manifestly and to use relations (\ref{hexa_coeffs}) and 
(\ref{hexa_cons}) to perform cancellations as far as possible. 


Applying the reduction formula (\ref{hexa_red}) above to reduce the 
hexagon, we observe that the coefficients of the hexagon and pentagon 
integrals in the amplitude combine in a nice way to form a 
resulting coefficient for a given pentagon which is again 
proportional to $b_j$. For example, the resulting coefficient of 
$I_5^n(p_{12},p_3,p_4,p_5,p_6)$ in (\ref{amp_si}) is 
\begin{eqnarray}
{\rm tr}(3456) + {\rm tr}(123456) \, b_1 = -2 s_{34}s_{45}s_{56} \, b_4\;,
\end{eqnarray}
analogous for all cyclic permutations. 
%Not only the size of the coefficients did not increase
%but also one can still use linear relations. 

Now we reduce the pentagons to boxes 
using the reduction formula given in\,\cite{Binoth:1999sp}. 
We obtain,
\begin{eqnarray}
{\cal A}(p_1,p_2,p_3,p_4,p_5,p_6) =
  \frac{2}{3} \, I_3^n(p_{12},p_{34},p_{56})  \hspace{6cm}\nonumber\\
 + s_{12}\, I_4^n(p_1,p_2,p_{34},p_{56}) 
+\frac{{\rm tr}(14)}{4}I_4^n(p_1,p_{23},p_{4},p_{56})\nonumber\\
+\frac{b_2}{2E_2}\Big\{
s_{23}s_{34}\,[\,{\rm tr}(1234)-2s_{12}\,(s_{234}-s_{23})]\,
I_4^n(p_2,p_{3},p_{4},p_{561})\nonumber\\
+s_{12}s_{23}\,[\,{\rm tr}(1234)-2s_{34}\,(s_{123}-s_{23})]\,
I_4^n(p_1,p_{2},p_{3},p_{456})\nonumber\\
+{\rm tr}(1234)\,E_2\,I_4^n(p_1,p_{23},p_{4},p_{56})\nonumber\\
+s_{34}\,[-s_{123}\,{\rm tr}(1234)-2s_{12}s_{23}\,(s_{123}-s_{56})]\,
I_4^n(p_{3},p_{4},p_{56},p_{12})\nonumber\\
+s_{12}\,[-s_{234}\,{\rm tr}(1234)-2s_{23}s_{34}\,(s_{234}-s_{56})]\,
I_4^n(p_1,p_2,p_{34},p_{56})\Big\}\nonumber\\
+\quad 5 \mbox{ cyclic permutations},
\label{amp_box}
\end{eqnarray}
where $E_1=s_{123}s_{345}-s_{12}s_{45}$. The $E_j$ 
for $j>1$ are defined by cyclic permutation.
Note that $E_{j}=E_{j+3}$.

The amplitude is now expressed in terms of four functions: The triangle 
with all three legs off-shell, box integrals with two off-shell legs at 
adjacent corners ($I_4^n(p_1,p_2,p_{34},p_{56})$ and 5 permutations), 
box integrals with two  off-shell legs at 
opposite corners ($I_4^n(p_1,p_{23},p_{4},p_{56})$ and 2 permutations), 
and box integrals with one off-shell leg 
($I_4^n(p_1,p_{2},p_{3},p_{456})$ and 5 permutations). 
%In order to see the cancellations,
 
We  now collect and combine the coefficients of particular terms in the 
cyclic sum. Again nontrivial cancellations
happen. In particular, using $\sum_{j=1}^6 b_j = 0$ and $\hat S\cdot b=1$ , 
the coefficients of the  box integrals with two  off-shell legs at 
opposite corners add up to zero and 
all dilogarithms related to box functions cancel! 
Hence the only terms which survive are the triangle graphs and
some logarithmic terms stemming from the finite parts of the 
box integrals, such that we finally obtain,
\begin{eqnarray}\label{amp_final}
{\cal A}(p_1,p_2,p_3,p_4,p_5,p_6) = G(p_1,p_2,p_3,p_4,p_5,p_6)
+ 5 \; \mbox{cyclic permutations},
\end{eqnarray}
with,
\begin{eqnarray} G(p_1,p_2,p_3,p_4,p_5,p_6) = 
 \frac{2}{3} \,I^n_3(p_{12},p_{34},p_{56})\hspace{6cm} \nonumber \\
+ \left\{  \frac{b_1}{E_1} 
\left[{\rm tr}(6123)-2s_{61}(s_{123}-s_{12}) \right] 
%\right. \hspace{2cm}\nonumber \\ \left. 
+\frac{b_2}{E_2} 
\left[ {\rm tr}(1234)-2s_{34}(s_{123}-s_{23}) \right] \right\} 
\hspace{1cm}\nonumber \\ \times
\log\left( \frac{s_{12}}{s_{123}} \right)
         \log\left( \frac{s_{23}}{s_{123}} \right) \nonumber \\
+\left\{ -b_1 + \frac{b_2}{2 E_2} 
\left[ {\rm tr}(1234) - 2 s_{34} (s_{123}-s_{23})  \right]
%\right. \hspace{3cm}\nonumber \\ \left.
+ \frac{b_6}{2 E_6}
\left[ {\rm tr}(5612)- 2 s_{56} (s_{345}-s_{61})\right]\right\}
\hspace{1cm}\nonumber \\ \times
\left[ \log\left( \frac{s_{12}}{s_{234}} \right) 
       \log\left( \frac{s_{56}}{s_{234}} \right)
     + \log\left( \frac{s_{34}}{s_{234}} \right) 
       \log\left( \frac{s_{12}}{s_{56}} \right) \right].\nonumber\\
       \label{ampsi}
\end{eqnarray}
Note that $G(p_1,p_2,p_3,p_4,p_5,p_6)$ has no spurious singularities. 
We checked that the numerator of expression (\ref{ampsi}) vanishes in 
the limits where its denominator  vanishes. 

Finally, the full amplitude is given by the sum over
permutations of the function $G$,
\begin{eqnarray}\label{Eq:final}
\Gamma_{\rm yuk}^{\phi}[p_1,p_2,p_3,p_4,p_5,p_6] 
= - \frac{g^6}{(4\pi)^2}\sum\limits_{\pi\in S_6}^{} G(p_{\pi_1},p_{\pi_2},
                  p_{\pi_3},p_{\pi_4},p_{\pi_5},p_{\pi_6}).
\end{eqnarray} 


\subsubsection{Summary}

The Yukawa model is a good testing ground to study 
nontrivial cancellations appearing in scalar integral reductions
without additional complications due to a nontrivial
tensor structure. Focusing on the massless case and $N=6$ 
we sketched the explicit calculation of the amplitude. 
It has been outlined how cancellations can be made manifest
at each step of the calculation by using linear relations
between reduction coefficients. 
With the present method there is no explosion of terms typical 
for multi--leg  calculations. The final answer is surprisingly compact
and contains -- apart from 3-point functions with 
3 off-shell legs -- only some products of logarithms.

As a next step more realistic examples have to be considered
including gauge bosons and a nontrivial infrared structure. 
It is justified to speculate that
the recombination of scalar integrals
will work similarly, such that
efficient algorithms to calculate 
six-point amplitudes at one loop are in reach.  
Work on this subject is in preparation.

\subsection{NNLO}
\label{sec:ho;nnlo}

\subsubsection{Motivation}
\label{sec:ho;nnlorationale}

There are many reasons why extending perturbative calculations to NNLO is
vital in reducing the theoretical uncertainty.  In the following we list five of
them.

\vspace{0.2cm}

\noindent 
{\it Renormalisation scale uncertainty}

\vspace{0.15cm}

\noindent 
In many cases, the uncertainty from the pdf's and from the choice of the
renormalisation scale $\mu_R$ give uncertainties that are as big or bigger than
the experimental errors.  Of course, the theoretical prediction should be
independent of $\mu_R$.   However, a scale dependence is introduced by
truncating the perturbative series. The change due to varying the scale is
formally higher order.
If an observable ${\cal O}bs$ is known to order $\alpha_s^N$ then,
$$
\frac{\partial}{\partial \ln(\mu_R^2)} 
\sum_{0}^N A_n(\mu_R) \alpha_s^n(\mu_R) = {\cal O} \left(\alpha_s^{N+1} \right).
$$
Often the uncertainty due to uncalculated higher
orders is estimated by varying the renormalisation scale upwards and downwards
by a factor of two around a typical hard scale in the process.  
However, the variation only produces copies of the lower order terms,
\\
e.g.  $$
{\cal O}bs = A_0 \alpha_s(\mu_R) + \left(A_1 + b_0 A_0
\ln\left({\mu_R^2\over \mu_0^2}\right)\right)
\alpha_s(\mu_R)^2.$$
$A_1$ will contain generally contain infrared logarithms and constants that are
not present in $A_0$ and  therefore {\em cannot be predicted} by varying
$\mu_R$. For example, $A_0$ may contain infrared logarithms $L$ up to $L^2$,
while $A_1$ would contain these logarithms up to $L^4$.  $\mu_R$ variation is
{\em only an estimate} of higher order terms A large variation probably means
that {\em predictable} higher order terms are large.

To illustrate the improvement in scale uncertainty that may occur at NNLO, let
us consider the production of a central jet in $p\bar p$ collisions.   The
renormalisation scale dependence is entirely predictable,
\begin{eqnarray*}
\frac{d\sigma}{dE_T} &=& \alpha^2_s(\mu_R) {A_0}\\
&+&  \alpha^3_s(\mu_R) \left ({A_1} + 2 b_0 L {A_0}\right)\\
&+&  \alpha^4_s(\mu_R) \left ({A_2} + 3 b_0 L {A_1} + (3b_0^2L^2+2 b_1 L) { A_0} \right)
\end{eqnarray*}
with $ L = \log(\mu_R/E_T)$.  $A_0$ and $A_1$ are the known LO and NLO coefficients
while $A_2$ is not presently known.  Inspection of Fig.~\ref{fig:nnlo} shows
that the scale dependence is systematically reduced by increasing the number of
terms in the perturbative expansion.   At NLO, there is always a turning point
where the prediction is insenstitive to small changes in $\mu_R$.  If this
occurs at a scale far from the typically chosen values of $\mu_R$, the
$K$-factor (defined as $ K = 1 + \alpha_s(\mu_R) {A_1}/{A_0}$) will be
large.  At NNLO the scale dependence is clearly significantly reduced, although
a more quantitative statement requires knowledge of $A_2$. 
\begin{figure}[t]
\label{fig:nnlo}
\begin{center}
\includegraphics[height=7cm]{nnlo1.ps}
\end{center}
\caption{Single jet inclusive distribution at $E_T = 100$~GeV and $0.1 < |\eta| < 0.7$ at
$\sqrt{s} = 1800$~GeV at LO (green), NLO (blue) and NNLO (red). The solid and dashed lines
how the NNLO prediction if $A_2=0$, $A_2= \pm A_1^2/A_0$ respectively. 
The same pdf's and $\alpha_s$ are used throughout.}
\end{figure}

\vspace{0.2cm}

\noindent
{\it Factorisation scale dependence}

\vspace{0.15cm}

\noindent
Similar qualitative arguments can be applied to the factorisation scale inherent
in perturbative predictions for quantities with initial state hadrons. 
Including the NNLO contribution reduces the uncertainty due to the truncation of
the perturbative series.

\vspace{0.2cm}

\noindent
{\it Jet algorithms}

\vspace{0.15cm}

\noindent
There is also a mismatch between the number of hadrons and the number of
partons in the event.   At LO each parton has to model a jet and there is no
sensitivity to the size of the jet cone.   At NLO two partons can combine
to  make a jet giving sensitivity to the shape and size of the jet cone.
Perturbation theory is starting to reconstruct the parton shower within the
jet. This is further improved at NNLO where up to three partons can form a
single jet, or alternatively two of the jets may be formed by two partons. This
may lead to a better matching of the jet algorithm between theory and 
experiment.

\vspace{0.2cm}

\noindent
{\it Transverse momentum of the incoming partons}

\vspace{0.15cm}

\noindent
At LO, the incoming particles have no transverse momentum with respect to the
beam so that the final state is produced at rest in the transverse plane.  At
NLO, single hard radiation off one of the incoming particles gives the final
state a transverse momentum kick even if no additional jet is observed. In some
cases, this is insufficient to describe the data and one appeals to the
intrinsic transverse motion of the partons confined in the proton to explain
the data.    However, at NNLO, double radiation from one particle or  single
radiation off each incoming particle gives more complicated transverse momentum
to the final state and may provide a better, and more theoretically motivated, 
description of the data.

\vspace{0.2cm}

\noindent
{\it Power corrections}
\begin{figure}[t]
\label{fig:powerc}
\begin{center}
\includegraphics[height=6cm]{powerc.ps}
\end{center}
\caption{The average value of $\langle 1-T\rangle$ given by Eq.~\ref{eq:omt}
showing the NLO prediction (dashed red), the NLO prediction with power correction of
$\lambda = 1$~GeV (solid blue) and an NNLO estimate with $a=3$ and a power correction
of $\lambda=0.5$~GeV (green dots).}
\end{figure}
Current comparisons of NLO predictions with experimental data generally reveal
the need for power corrections.   For example, in
electron-positron annihilation,  the experimentally measured average value of
1-Thrust lies well above the NLO predictions.   The difference can be accounted
for by a $1/Q$ power correction. While the form of the power correction can be
theoretically motivated, the magnitude is generally extracted from data and, to
some extent, can be attributed to uncalculated higher orders.   Including the
NNLO may therefore reduce the size of the phenomenological  power correction
needed to fit the data.

Before the calculation of the NNLO contribution it is not possible to make a more
quantitative statement. However to illustrate the qualitative point, 
let us take the simple example of an observable like
$\langle 1-T\rangle$ which can be modelled by the simplified series,
\begin{equation}
\label{eq:omt}
\langle 1-T\rangle = 0.33 \alpha_s(Q) + 1.0 \alpha_s(Q)^2 + a \alpha_s(Q)^3 +
\frac{\lambda}{Q},
\end{equation}
with $\alpha_s(Q) \sim 6\pi/23/\log(Q/\Lambda)$ and $\Lambda = 200$~MeV. Fig~\ref{fig:powerc} shows the NLO
perturbative prediction $a=0$,  $\lambda=0$ as well as the NLO prediction
combined with a power correction, $a=0$, $\lambda=1$~GeV which can be taken to
model the data.   If the NNLO coefficient turns out to be positive (which is by
no means guaranteed), then the size of the power correction needed to describe
the data will be reduced. For example, if we estimate the NNLO coefficient as
$a = 3$, which is large but perhaps not unreasonable, then the NLO prediction
plus power correction can almost exactly be reproduced with a power correction
of the same form, but $\lambda  = 0.6$~GeV.  We are effectively trading a
contribution of ${\cal O}(1/Q)$ for a contribution of $1/\log^3(Q/\Lambda)$.  
At present the data is insufficient to distinguish between these two functional
forms. 

\subsubsection{Parton densities at NNLO}
\label{sec:ho;nnlopdf}

Consistent NNLO predictions for processes involving hadrons in the  initial
state require not only the NNLO hard scattering  cross sections, but also
parton distribution functions which are accurate to this order. 

The evolution of parton distributions is governed at NNLO by the  three-loop
splitting functions, which are not fully known at present.  However, using the
available information on some of the lower Mellin 
moments~\cite{Larin:1994vu,Larin:1997wd,Retey:2000nq} and on the asymptotic
behaviour~\cite{Catani:1994sq},  
as well as some exactly known terms~\cite{Gracey:1994nn},
it is possible to  construct approximate expressions for these splitting 
functions~\cite{vanNeerven:2000wp}. These approximations (which are provided with an 
error band) can serve as a substitute until full results become 
available~\cite{Moch:2001fr}. 

The determination of NNLO parton distributions requires a global  fit to the
available data on a number of hard scattering  observables, with all
observables computed consistently at NNLO.  At present, the NNLO coefficient
functions are available only for  the inclusive  Drell-Yan
process~\cite{Hamberg:1991np,Harlander:2002wh} and for deep inelastic structure 
functions~\cite{Zijlstra:1992qd}. These two observables are by themselves  insufficient
to fully constrain all parton species. In the first  NNLO analysis of parton
distribution functions,  which was performed recently~\cite{Martin:2002dr}, these
processes were  therefore accompanied by several other observables only known
to NLO.  The resulting distributions~\cite{Martin:2002dr} illustrate some  important
changes in size and shape of the distribution functions  in going from NLO to
NNLO, visible in particular for the gluon distribution  function. 

It is clear that further progress on the determination of NNLO  parton
distributions requires a larger number of processes  (in particular jet
observables) to be treated consistently at NNLO. 



\subsubsection{Infrared limits of one-loop and tree-level processes}
\label{sec:ho;nnloIR}

For simplicity let us consider a process with no initial state partons and with
$m$ partons in the final state at LO.
The NNLO contribution to the cross section can be decomposed as 
\begin{equation}
\sigma^{NNLO} = \int_{m+2} d\sigma^{RR} + \int_{m+1} d\sigma^{VR} 
+ \int_{m}
d\sigma^{VV}.
\end{equation}
Here $\int_n$ denotes the $n$-particle final state, while $d\sigma^{XY}$ denotes
the fully differential cross section with double radiation ($XY=RR$), single
radiation from one-loop graphs ($XY=VR$) and the double virtual contribution
($XY=VV$) that includes both the square of the one-loop graphs and the
interference of tree and two-loop diagrams. 
After renormalisation of the virtual matrix elements, each of the contributions
is UV finite.   However,  each of these terms is separately infrared
divergent and this manifests itself as poles in $\epsilon$.
Nevertheless, the Kinoshita-Lee-Nauenberg theorem states that
the infrared singularities must cancel for sufficiently inclusive physical
quantities.  The problem is to isolate the infrared poles and analytically
cancel them before taking the $\e \to 0$ limit.
Establishing a strategy for doing this requires a good understanding of  
$ d\sigma^{RR}$ and $d\sigma^{VR}$ in the infrared region where the additional
radiated particle(s) are unresolved.
 
\vspace{0.2cm}

\noindent
{\it Double unresolved limits of tree amplitudes}

\vspace{0.15cm}

\noindent
The infrared singular regions of tree amplitudes can be divided into several
categories, 
\begin{enumerate} 
\item three collinear particles, 
\item two pairs of collinear particles, 
\item two soft particles, 
\item one soft and two collinear particles, 
\item a soft quark-antiquark pair. 
\end{enumerate} 
In each of these limits, the tree-level $m+2$ particle amplitudes factorise and
yield an infrared singular factor multiplying the tree-level $m$-particle
amplitude.   The various limits have been well studied.   

The limit where a quark, gluon and photon simultaneously become collinear was
first studied in \cite{Gehrmann-DeRidder:1998gf} and then extended for generic
QCD processes in Ref.~\cite{Campbell:1998hg}  by directly taking the limit of
tree-level matrix elements. These limits have been subsequently rederived using
general gauge invariant methods and extended to include azimuthal
correlations between the particles~\cite{Catani:1998nv,Catani:1999ss,DelDuca:1999ha}. 

When two independent pairs of particles  are collinear, the singular limits can
be treated independently, and just lead to a trivial extension of the NLO
result - the product of two single collinear splitting functions.

The limit where the momenta of two of the gluons simultaneously become soft has
also been studied~\cite{Berends:1989zn,Catani}. At the amplitude level, the singular
behaviour factorises in terms of a process independent soft two-gluon
current~\cite{Catani:1999ss} which is the generalisation of the one-gluon
eikonal current.

The soft-collinear limit occurs when the momentum of one gluon becomes soft
simultaneously with two other partons becoming collinear.   Factorisation
formulae in this limit have been provided in both the azimuthally averaged
case~\cite{Campbell:1998hg} and  including the angular
correlations~\cite{Catani:1998nv,Catani:1999ss,DelDuca:1999ha}.

Finally, when the momentum of both quark and antiquark of a $q \bar q$ pair
become soft, the tree-level amplitude again factorises~\cite{Catani:1999ss}.

Taken together, these factorisation formulae describe all of the cases where
$m$-partons are resolved in a $m+2$ parton final state. Techniques for
isolating the divergences have not yet been established and there is an
on-going effort to develop a set of {\em local} subtraction counter-terms that
can be analytically integrated over the infrared regions.


\vspace{0.2cm}

\noindent
{\it Single unresolved limits of one-loop amplitudes}

\vspace{0.15cm}

\noindent
The soft- and collinear-limits of one-loop QCD amplitudes has also been
extensively studied. In the collinear limit, the one-loop $m+1$ particle
process factorises as a  one-loop $m$ particle amplitude multiplied by a tree
splitting function together with a tree $m$ particle amplitude multiplied by a
one-loop splitting  function.  The explicit forms of the splitting amplitudes
were first determined to ${\cal O}(\e^0)$~\cite{Bern:1994zx}.  However, because
the integral over the infrared phase space generates additional poles, the
splitting functions have been determined to all orders in
$\e$~\cite{Bern:1998sc,Kosower:1999xi,Bern:1999ry,Kosower:1999rx}. Similarly,
when a gluon becomes soft, there is a factorisation of the one-loop amplitude
in terms of the one-loop soft
current~\cite{Bern:1998sc,Bern:1999ry,Catani:2000pi}. Because of the similarity
of the factorisation properties  with the single unresolved particle limits of 
tree-amplitudes, it is relatively straightforward to isolate the infrared poles
through the construction of a set of {\em local} subtraction counter-terms. In
Sec.~\ref{sec:ho;exercise}, we illustrate how the infrared singularities from
the single unresolved limits of one-loop amplitudes combine with the
predictable infrared pole structure of the virtual contribution for the
explicit example of $\gamma^* \to q\bar q$.


\subsubsection{Two-loop matrix elements for scattering processes}
\label{sec:ho;two-loop}

In recent years, considerable progress has been made on the  calculation of
two-loop  virtual corrections to the multi-leg matrix elements relevant for 
jet physics, which describe either $1\to 3$ decay or $2\to 2$ scattering 
reactions.  Much of this progress is due to a number of important technical 
breakthroughs related to the evaluation of the large number of  different
integrals appearing in the two-loop four-point amplitudes.

It should be recalled that  perturbative corrections to many inclusive quantities
have been computed  to the two- and three-loop level already several years ago. 
From the technical  point of view, these inclusive calculations correspond to the
computation of  multi-loop two-point functions, for which many elaborate
calculational tools have been developed.  Using dimensional
regularization~\cite{Bollini:1972ui,Ashmore:1972uj,Cicuta:1972jf,'tHooft:1972fi}  with $d\neq 4$
dimensions as regulator for ultraviolet and infrared divergences,  the large number
of different integrals appearing in  multi-loop two-point functions can be reduced
to a small number of  so-called {\em master integrals} by using 
integration-by-parts
identities~\cite{'tHooft:1972fi,Tkachov:1981wb,Chetyrkin:1981qh}. These identities 
exploit the fact that the integral over the total derivative of any  of the loop
momenta vanishes in dimensional regularization. 

Integration-by-parts identities can also be obtained for  integrals appearing
in amplitudes with more than two external legs;  for these amplitudes, another
class of identities exists due to  Lorentz invariance of the amplitudes. These 
Lorentz invariance identities~\cite{Gehrmann:1999as} rely on the fact that  an
infinitesimal Lorentz transformation commutes with the loop integral, thus
relating different integrals. Using integration-by-parts  and Lorentz
invariance identities, all  two-loop Feynman amplitudes for $2\to 2$ scattering
or $1\to 3$ decay  processes  can be expressed as linear combinations of a
small number of  master integrals,
which have  to be computed by some different method.  Explicit
reduction formulae for on-shell two-loop four-point  integrals were derived
in~\cite{Smirnov:1999wz,Anastasiou:1999bn,Anastasiou:2000mf}. Computer algorithms for the  automatic reduction of all
two-loop four-point integrals  were described in~\cite{Gehrmann:1999as,Laporta:2001dd}. 

A related development was the proof of the equivalence of  integration-by-parts
identities for integrals with the same  total number of external and loop
momenta~\cite{Baikov:2000jg}. Consequently,  much of the tools developed for
the computation of three-loop propagator
integrals~\cite{Gorishnii:1989gt,Chetyrkin:1996ia} can be readily applied to
two-loop  vertex functions. As a first application, the two-loop  QCD
corrections to Higgs boson production in gluon-gluon fusion were 
computed~\cite{Harlander:2000mg} in the limit of large top quark mass. This
result allowed the complete NNLO
description~\cite{Catani:2001ic,Harlander:2001is,Harlander:2002wh} of inclusive
Higgs production at hadron colliders. 

The master integrals relevant to $2\to 2$ scattering or $1\to 3$ decay 
processes are massless, scalar two-loop  four-point functions with all legs
on-shell or a single leg off-shell.  Several techniques for the computation of
those functions  have been proposed in the literature, such as the application
of a Mellin-Barnes transformation to all propagators~\cite{Smirnov:1999gc,Tausk:1999vh} 
or the negative dimension approach~\cite{Anastasiou:1999cx,Suzuki:2001yf}.  Both techniques
rely on an explicit integration  over the loop momenta, with differences mainly
in the representation used for the propagators. These techniques were used
successfully to compute a number of master integrals. Employing the
Mellin-Barnes method, the on-shell planar double box
integral~\cite{Smirnov:1999gc,Anastasiou:2000kp}, the on-shell non-planar double box
integral~\cite{Tausk:1999vh}  and two double box
integrals with one leg off-shell~\cite{Smirnov:2000vy,Smirnov:2000ie}
were computed. Most recently, the same method was used to  derive the on-shell
planar double box integral~\cite{Smirnov:2001cm} with one internal mass scale.  The
negative dimension approach has been applied~\cite{Anastasiou:1999cx}  to compute the
class of two-loop box integrals which correspond to a one-loop bubble insertion
in one of the propagators of the one-loop box. 

A method for the analytic computation of master integrals avoiding the explicit
integration over the loop momenta is to derive differential equations in 
internal propagator masses or in external momenta for the master integral,  and
to solve these with appropriate boundary conditions.  This method has first
been suggested by Kotikov~\cite{Kotikov:1991pm} to relate  loop integrals with
internal masses to massless loop integrals.  It has been elaborated in detail
and generalized to differential  equations in external momenta
in~\cite{Remiddi:1997ny}; first  applications were presented
in~\cite{Caffo:1998du,Caffo:1999nk}. The
computation of master integrals from differential equations proceeds as 
follows. Carrying out the derivative with respect to an external invariant  on
the master integral of a given topology, one obtains a linear combination of a
number of  more complicated integrals, which can however be reduced to the 
master integral itself plus simpler integrals by applying the reduction 
methods discussed above. As a result, one obtains  an inhomogeneous linear
first order differential equation in each invariant for the master integral.

The inhomogeneous term in these differential equations contains only 
topologies simpler than the topology under consideration, which are considered
to be known  if working in a bottom-up approach.  The master integral is then
obtained by matching the general  solution of its differential equation to an
appropriate boundary condition. Quite in general, finding a boundary condition
is  a  simpler problem than evaluating the whole integral, since it depends on
a smaller number of kinematical variables. In some cases, the boundary
condition can even be determined from the differential equation itself. 

Using the differential equation technique, one of the on-shell 
planar double box integrals~\cite{Gehrmann:2000xj} as well as the full set of 
planar and non-planar off-shell double box integrals~\cite{Gehrmann:2000zt,Gehrmann:2001ck}
 were derived. 

A strong check on all these computations of master integrals is  given by the
completely numerical calucations of~\cite{Binoth:2000ps}, which are  based on an iterated
sector decomposition to isolate the infrared pole  structure. The methods
of~\cite{Binoth:2000ps} were applied to confirm {\em all}  of the above-mentioned
calculations. 

The two-loop four-point functions with all legs on-shell can be expressed  in
terms of Nielsen's polylogarithms~\cite{nielsen,lewin,Kolbig:1986qt,bit}.  In
contrast, the closed analytic  expressions for  two-loop four-point functions
with one leg off-shell contain  two new classes of functions: harmonic 
polylogarithms~\cite{Remiddi:1999ew,Gehrmann:2001pz} and two-dimensional
harmonic  polylogarithms (2dHPL's)~\cite{Gehrmann:2001jv}.  Accurate numerical
implementations for these functions~\cite{Gehrmann:2001pz,Gehrmann:2001jv} are
available. 

\vspace{0.2cm}

\noindent
{\it $2 \to 2$ Processes with all legs on-shell}

\vspace{0.15cm}

\noindent
With the explicit solutions of the integration-by-parts and  Lorentz-invariance
identities for on-shell two-loop four-point 
functions~\cite{Smirnov:1999wz,Anastasiou:1999bn,Anastasiou:2000mf} and the
corresponding master  
integrals~\cite{Gonsalves:1983nq,Kramer:1987sr,Smirnov:1999gc,Tausk:1999vh,Anastasiou:2000kp,Gehrmann:2000xj}, all
necessary ingredients  for the computation of two-loop corrections to  $2 \to 2$
processes with all legs on-shell are now available. In fact, only half a  year
elapsed between the completion of the full set of master 
integrals~\cite{Anastasiou:2000kp,Gehrmann:2000xj} and the calculation of the 
two-loop QED corrections to Bhabha-scattering~\cite{Bern:2000ie}.  Subsequently,
results were obtained for the two-loop  QCD corrections to all parton-parton
scattering processes~\cite{Anastasiou:2000kg,Anastasiou:2000ue,Anastasiou:2001sv,
Glover:2001af}. For gluon-gluon scattering, the two-loop helicity amplitudes have
also  been derived~\cite{Bern:2000dn,Bern:2002tk}. Moreover, two-loop corrections were
derived  to processes involving two partons and two real
photons~\cite{Bern:2001df,Anastasiou:2002zn}. Finally, light-by-light scattering in
two-loop QED and QCD was  considered in~\cite{Bern:2001dg}, these results were
extended to supersymmetric QED in~\cite{Binoth:2002xg}. 

It turns out that supersymmetry can provide strong checks on the  consistency
of the matrix elements. Calculations in this framework  do however require
modifications to the dimensional reguarization scheme, which were discussed in
detail in~\cite{Bern:2002zk}. 

At the same order in perturbation theory
as the two-loop matrix elements (which are 
obtained by contracting the two-loop and the tree level amplitudes),
one also finds contributions from the square of the one-loop amplitude. 
The evaluation of these contributions uses well-known one-loop techniques and 
is straightforward. For parton-parton scattering, these one-loop
self-interference contributions were computed 
in~\cite{Anastasiou:2000mv,Anastasiou:2001sv,Glover:2001rd}. 

The results for the two-loop QED matrix element for  Bhabha
scattering~\cite{Bern:2000ie} were used in~\cite{Glover:2001ev} to extract the  single
logarithmic contributions to the Bhabha scattering cross section,  thus
improving considerably on the accuracy of earlier~\cite{Arbuzov:1997qd} double
logarithmic results.




\vspace{0.2cm}

\noindent
{\it $2 \to 2$ Processes with one off-shell leg}

\vspace{0.15cm}

\noindent
With the reduction algorithm~\cite{Gehrmann:1999as} and full set of two-loop master 
integrals~\cite{Gehrmann:2000zt,Gehrmann:2001ck} for two-loop four-point functions  now available, it
is possible to compute the two-loop corrections  a number of $1\to 3$ decay and
$2\to 2$ scattering reactions with one  off-shell leg. 


As a first result (and initiated during this workshop),
in~\cite{Garland:2001tf},  the two-loop QCD matrix element for  $e^+e^- \to
3$~jets and the  corresponding one-loop self-interference matrix element were
derived. The three jet production rate in electron-positron  collisions and
related event shape observables are in fact the  most precisely measured jet 
observables at present~\cite{Bethke:2000ai}. They will also play an important
role in future  QCD studies~\cite{Aguilar-Saavedra:2001rg} at a the proposed
high energy linear $e^+e^-$ collider. 

It is worthwhile to note that  besides its  phenomenological importance,  the
three-jet rate has also served as a theoretical testing ground for the
development of new techniques for higher order calculations in QCD: both the
subtraction~\cite{Ellis:1981wv} and the phase-space
slicing~\cite{Fabricius:1981sx} methods for the extraction of infrared
singularities from  NLO real radiation  processes were developed in the context
of the first three-jet calculations. The systematic formulation of 
phase-space  slicing~\cite{Giele:1992vf} as  well as the dipole
subtraction~\cite{Catani:1997vz} method were also first demonstrated for
three-jet observables, before being applied to other processes. It is very
likely that similar techniques at higher orders will first be developed in the
context of jet production in $e^+e^-$ annihilation, which in contrast to
hadron--hadron collisions or electron--proton scattering does not pose the
additional difficulty of the regularization of  initial state singularities.  



Processes related to  $e^+e^- \to 3$~jets  by crossing symmetry  are
$(2+1)$-jet production in  deep inelastic $ep$ scattering and
vector-boson-plus-jet production at hadron  colliders.  Crossing of the 
$e^+e^- \to 3$~jets two-loop matrix element  to the kinematic regions relevant
for these scattering processes  requires the analytic continuation of 2dHPL's
outside their range of allowed arguments. This topic is currently under
investigation. 



\subsubsection{Infrared structure for $\gamma^* \to 2$~jets at NNLO}
\label{sec:ho;exercise}

As an exercise in how to combine the various next-to-next-to-leading order
contributions, we consider the infrared singularity structure for  $\gamma^*
\to 2$~jets at NNLO.   This is the simplest process that we can imagine at
NNLO.   It involves the one- and two-loop two parton amplitudes, the one-loop
three parton amplitudes and the tree-level four parton amplitudes.    The aim
is to symbolically identify the origins of the infrared singularities of the
single and double radiation graphs.  Once this is achieved,  it should be
possible to construct {\em local} counter terms that can be subtracted
numerically from the three- and four-parton matrix elements to give infrared
finite contributions  and analytically integrated and combined with the
explicit singularity structure of the two-parton contribution.



\vspace{0.2cm}

\noindent
{\it The one- and two-loop two parton contribution}

\vspace{0.15cm}

\noindent
The renormalized $\gamma^* \to q\bar q$
amplitude can be written as
\begin{equation}
|\cmqq\rangle = \sqrt{4\pi\alpha}e_q  \left[
|\cmqq^{(0)}\rangle 
+ \asmuopi |\cmqq^{(1)}\rangle 
+ \asmuopi^2 |\cmqq^{(2)}\rangle 
+ {\cal O}(\alpha_s^3) \right] \;,
\end{equation}
where $\alpha$ denotes the electromagnetic coupling constant, 
$e_q$ the quark charge,
and the $|{\cal M}_{q\bar q}^{(i)}\rangle$ are the $i$-loop contributions to the 
renormalized amplitude. They are scalars in colour space. 

The squared amplitude, summed over spins, colours and quark flavours, 
is denoted by
\begin{equation}
\langle\cmqq|\cmqq\rangle = \sum |\cm (\gamma^* \to q\bar q)|^2 
= {\cal A}_{q\bar q}\; .
\end{equation}
The perturbative expansion of ${\cal A}_{q\bar q}$ at renormalization scale 
$\mu^2 = \sab$ reads,
\begin{equation}
{\cal A}_{q\bar q} = 4\pi\alpha\sum_q N e_q^2 \Bigg[\caqq{0}
+\asmuopi \caqq{2} +\asmuopi^2 \caqq{4} + {\cal O}(\alpha_s^3)\Bigg] \;,
\end{equation}
where $N$ is the number of colours
and where, 
\begin{eqnarray}
\label{eq:Aqq0}
\caqq{0} &=& \langle\cmqq^{(0)}|\cmqq^{(0)}\rangle 
= 4 (1-\e) \sab \;,\\
\label{eq:Aqq2}
\caqq{2} &=& 
\langle\cmqq^{(0)}|\cmqq^{(1)}\rangle +
\langle\cmqq^{(1)}|\cmqq^{(0)}\rangle \nonumber \\
\caqq{4} &=& 
\langle\cmqq^{(1)}|\cmqq^{(1)}\rangle +
\langle\cmqq^{(0)}|\cmqq^{(2)}\rangle +
\langle\cmqq^{(2)}|\cmqq^{(0)}\rangle \;,
\label{eq:Aqq4}
\end{eqnarray}
with~\cite{Kramer:1987sg,Matsuura:1988wt,Matsuura:1989sm},
\vspace{-0.5cm}
\begin{eqnarray}
\langle\cmqq^{(0)}|\cmqq^{(0)}\rangle 
&\equiv& \includegraphics[width=2.5cm]{m200.ps} , \\
\langle\cmqq^{(0)}|\cmqq^{(1)}\rangle 
&\equiv& \includegraphics[width=2.5cm]{m210.ps} \nonumber \\
&=& \includegraphics[width=2.5cm]{m200.ps} \Bigg\{ \frac{1}{2}\, C_F\Bigg[
-\frac{2}{\e^2}-\frac{3}{\e} +\left( \frac{7\pi^2}{6}-8\right) 
\nonumber \\ &&
\hspace{3.5cm}
+\left(-16+\frac{7\pi^2}{4} +\frac{14}{3}\zeta_3\right)\, \e 
\nonumber \\ &&
\hspace{3.5cm}
+\left(-32 +\frac{14\pi^2}{3}+ 7\zeta_3-\frac{73\pi^4}{720}\right) \, \e^2
+ {\cal O}(\e^3)\Bigg] \Bigg\},\\
\langle\cmqq^{(1)}|\cmqq^{(1)}\rangle 
&\equiv& \includegraphics[width=2.5cm]{m211.ps} \nonumber \\
&=& \includegraphics[width=2.5cm]{m200.ps} \Bigg\{  C_F^2\Bigg[
\frac{1}{\e^4} + \frac{3}{\e^3} 
+ \left(\frac{41}{4}-\frac{\pi^2}{6}\right)\frac{1}{\e^2}
 + \left(28 -\frac{\pi^2}{2}-\frac{14}{3}\zeta_3\right)\frac{1}{\e} 
\nonumber \\ && 
\hspace{3.5cm}
+ \left(         72
          - \frac{41\pi^2}{24}
          - \frac{7\pi^4}{120}
          - 14 \zeta_3 \right)
+ {\cal O}(\e)
\Bigg] \Bigg\} ,\\
\langle\cmqq^{(0)}|\cmqq^{(2)}\rangle 
&\equiv& \includegraphics[width=2.5cm]{m220.ps} \nonumber \\
&=& \includegraphics[width=2.5cm]{m200.ps} \Bigg\{ 
\frac{1}{2\e^4} \, C_F^2  + \frac{1}{\e^3}\left[
\frac{3}{2}\, C_F^2 +\frac{11}{8} C_FC_A - \frac{1}{2} C_F T_R N_F
\right]\nonumber \\
&&\hspace{3cm}
     + \frac{1}{\e^2}\left[ 
\left( \frac{41}{8} - \frac{13\pi^2}{12}\right) C_F^2
+ \left( \frac{4}{9} +\frac{\pi^2}{24} \right) C_FC_A
-\frac{2}{9} C_F T_R N_F \right]\nonumber \\
&&\hspace{3cm}
     + \frac{1}{\e}\Bigg[ 
\left( \frac{221}{16} - 3\pi^2 -\frac{16}{3}\zeta_3\right) C_F^2
+ \left( \frac{65}{108} +\frac{\pi^2}{12} \right) C_F T_R N_F
\nonumber \\
&&\hspace{3.5cm}
+ \left( -\frac{961}{432} -\frac{11\pi^2}{48} +\frac{13}{4}\zeta_3\right)
                       C_FC_A\Bigg]\nonumber \\
&&\hspace{3cm}
+ \Bigg[
\left( \frac{1151}{32} - \frac{475\pi^2}{48} 
-\frac{29}{2}\zeta_3+\frac{59\pi^4}{144}\right) C_F^2
\nonumber \\
&&\hspace{3.5cm}
+ \left( -\frac{51157}{2592} +\frac{1061\pi^2}{432} 
   +\frac{313}{36}\zeta_3 - \frac{\pi^4}{45}\right)      C_FC_A
\nonumber \\
&&\hspace{3.5cm}
+ \left( \frac{4085}{648} -\frac{91\pi^2}{108} +\frac{1}{9} \zeta_3\right)
 C_F T_R N_F \Bigg]
+ {\cal O}(\e)
\Bigg\}.
\end{eqnarray}
\vspace{0.5cm}

\noindent We see that the NLO part $\caqq{2}$ is proportional to $\CF$ while the NNLO
contribution $\caqq{4}$ contains $\CF^2$, $\CF\CA$ and $\CF T_R \NF$ parts.

Catani~\cite{Catani:1998bh} has shown how to organize the 
infrared pole structure of the one- and two-loop contributions 
renormalized in the 
\MSbar\ scheme in terms of the tree and renormalized one-loop amplitudes,
$\ket{\cmqq^{(0)}}$ and $\ket{\cmqq^{(1)}}$ respectively by introducing two infrared
operators, $\Ione$ and $\Htwo$.  

The form of $\Ione$ can be simply understood by considering the total cross
section to ${\cal O}(\alpha_s)$ which is finite,
\vspace{-0.5cm}
\begin{equation}
\includegraphics[width=2.5cm]{m210.ps} 
+
\includegraphics[width=2.5cm]{m201.ps} 
+
\includegraphics[width=2.5cm]{m300.ps} 
= {\rm finite}.
\label{eq:qqnlo}
\end{equation}
\vspace{0.5cm}

\noindent Taking soft and collinear limits of the real radiation graphs, 
and integrating over the unresolved phase space, we see that,
\vspace{-0.5cm}
\begin{eqnarray}
\includegraphics[width=2.5cm]{m300.ps} 
&\longrightarrow& 
\includegraphics[width=2.5cm]{m200I1a.ps}
+
\includegraphics[width=2.5cm]{m200I1b.ps}
+{\rm finite}
\nonumber \\
&=& 
-2 \Ione \includegraphics[width=2.5cm]{m200.ps} 
+{\rm finite}.
\label{eq:findi1}
\end{eqnarray}
\vspace{0.3cm}

\noindent
The symbol $\includegraphics[width=2cm]{dot.ps}$
 represents the taking of the
infrared limit and integrating over the infrared phase space.   This operation
can be carried out in many ways. The most commonly used procedures are 
the subtraction formalism~\cite{Ellis:1981wv},
phase space slicing method~\cite{Fabricius:1981sx}, as well as 
systematic procedures~\cite{Giele:1992vf} and improved 
formulations of these methods, such as the dipole 
subtraction formalism~\cite{Catani:1997vz}. 


In general, $\Ione$ contains colour correlations. However, 
for this particular process, there is only one colour structure present at
tree level which is simply $\delta_{ij}$. Adding higher loops does not
introduce additional colour structures into the amplitude, and $\Ione$ is therefore a
$1 \times 1$ matrix in colour space.  It is proportional to $\CF$ and
is given by
\begin{equation}
\Ione
=
- \frac{e^{\epsilon\gamma}}{2\Gamma(1-\epsilon)}  
\CF
\left(\frac{2}{\epsilon^2}+\frac{3}{\epsilon}\right)
\left(-\frac{\mu^2}{s_{12}+i0}\right)^{\epsilon}  ,\label{eq:I1}
\end{equation}
with $\mu^2 = \sab$.  The double pole is due to the soft gluon singularity while
the single pole is derived from the collinear quark-gluon splitting function.
Inserting Eq.~(\ref{eq:findi1}) into Eq.~(\ref{eq:qqnlo}) yields the infrared
singular behaviour of the one-loop amplitude,
\vspace{-0.5cm}
\begin{equation}
\includegraphics[width=2.5cm]{m210.ps} 
= 
\Ione\includegraphics[width=2.5cm]{m200.ps} 
+{\rm finite}.
\end{equation}
\vspace{0.5cm}
\noindent
According to Catani~\cite{Catani:1998bh}, the infrared behaviour of the square of the
one-loop
contribution is given by,
\vspace{-0.5cm}
\begin{equation}
\includegraphics[width=2.5cm]{m211.ps}=
2\Re \Biggl[
-\frac{1}{2}\left| \Ione \right | ^2 \includegraphics[width=2.5cm]{m200.ps}
+\Ione
\includegraphics[width=2.5cm]{m201.ps}  
\Biggr ] + {\rm
finite},
\label{eq:polesa}
\end{equation}
\vspace{0.5cm}
where $\Ione^\dag$ is obtained from Eq.~(\ref{eq:I1}) by reversing the sign of
$+i0$.  The two-loop contribution can be decomposed as,
\vspace{-0.5cm}
\begin{eqnarray}
\lefteqn{\includegraphics[width=2.5cm]{m220.ps} + 
\includegraphics[width=2.5cm]{m202.ps}}\nonumber \\
 &=& 2\Re \Biggl [
-\frac{1}{2} \left(\Ione\right)^2 
\includegraphics[width=2.5cm]{m200.ps}  
+ \Ione \includegraphics[width=2.5cm]{m210.ps}
\nonumber \\
&&
-\frac{\beta_0}{\e} 
\Ione \includegraphics[width=2.5cm]{m200.ps}
+ e^{-\epsilon \gamma } \frac{ \Gamma(1-2\epsilon)}{\Gamma(1-\epsilon)} 
\left(\frac{\beta_0}{\epsilon} + K\right)\Ionep
\includegraphics[width=2.5cm]{m200.ps}\nonumber \\
&&
+\Htwo \includegraphics[width=2.5cm]{m200.ps} \Biggr ] + {\rm finite},
\label{eq:polesb}
\end{eqnarray}
with 
\begin{equation}
\beta_0 = \frac{11 \CA - 4 T_R \NF}{6},\qquad
K = \left( \frac{67}{18} - \frac{\pi^2}{6} \right) \CA - 
\frac{10}{9} T_R \NF.
\end{equation}
We note that the first line of the RHS of Eq.~(\ref{eq:polesb}) is
 proportional to
$\CF^2$ and produces fourth order poles in $\e$, 
while the second line contains terms like $\CF\CA$ and $\CF\TF\NF$ and generates 
third order poles in $\e$.
The last term of Eq.~(\ref{eq:polesb}) involving $\Htwo$
produces only a single pole in $\epsilon$ and is given by 
\begin{equation}
\label{eq:htwo}
\bra{\cm^{(0)}}\Htwo\ket{\cm^{(0)}} 
=\frac{e^{\epsilon \gamma}}{4\,\epsilon\,\Gamma(1-\epsilon)} H^{(2)} 
\braket{\cm^{(0)}}{\cm^{(0)}} \;,  
\end{equation}
where the constant $H^{(2)}$ is renormalization-scheme-dependent.
As with the single pole parts of $\Ione$,
the process-dependent
$H^{(2)}$ can be constructed by counting the number of
radiating partons present in the event.
In our case, there is only a 
 quark--antiquark pair present in the final
state, so that 
\begin{equation}
H^{(2)} =  2H^{(2)}_{q},
\end{equation}
where in the \MSbar\ scheme
\begin{equation}
H^{(2)}_q =       
         C_F^2  \left(
          - \frac{3}{8}
          - 6 \zeta_3
          + \frac{\pi^2}{2}
          \right)
       + C_F C_A  \left(
            \frac{245}{216}
          + \frac{13}{2} \zeta_3
          - \frac{23\pi^2}{48}
          \right) 
       + C_F T_R N_F  \left(
          - \frac{25}{54}
          + \frac{\pi^2}{12}
          \right)\;.
\end{equation}
The singularities present in Eqs.~(\ref{eq:polesa}) and (\ref{eq:polesb}) must
be cancelled against the infrared poles present in the one-loop three parton
and tree-level four parton contributions.

\vspace{0.2cm}

\noindent
{\it The one-loop three parton contributions}

\vspace{0.15cm}

\noindent
The renormalized $\gamma^* \to q\bar q g$
amplitude can be written as
\begin{equation}
|\cmqqg\rangle = \sqrt{4\pi\alpha}e_q \sqrt{\asmuopi} \left[
|\cmqqg^{(0)}\rangle 
+ \asmuopi |\cmqqg^{(1)}\rangle 
+ {\cal O}(\alpha_s^2) \right] \;,
\end{equation}
where the $|{\cal M}_{q\bar qg}^{(i)}\rangle$ are the $i$-loop contributions to the 
renormalized amplitude.   

The squared three-parton amplitude, summed over spins, 
colours and quark flavours, 
is denoted by
\begin{equation}
\langle\cmqqg|\cmqqg\rangle = \sum |\cm (\gamma^* \to q\bar qg)|^2 
= {\cal A}_{q\bar qg}\; .
\end{equation}
The perturbative expansion of ${\cal A}_{q\bar qg}$ at renormalization scale 
$\mu^2 = \sab$ reads,
\begin{equation}
{\cal A}_{q\bar qg} = 4\pi\alpha\sum_q N e_q^2 \Bigg[ 
+\asmuopi \caqqg{2} +\asmuopi^2 \caqqg{4} + {\cal O}(\alpha_s^3)\Bigg] \;,
\end{equation}
where,\vspace{-0.5cm}
\begin{eqnarray}
\label{eq:Aqqg2}
\caqqg{2} &=& \langle\cmqqg^{(0)}|\cmqqg^{(0)}\rangle 
=\includegraphics[width=2.5cm]{m300.ps},\\
&& \nonumber \\
\label{eq:Aqqg4}
\caqqg{4} &=& 
\langle\cmqqg^{(0)}|\cmqqg^{(1)}\rangle +
\langle\cmqqg^{(1)}|\cmqqg^{(0)}\rangle \nonumber \\
&=& 
\includegraphics[width=2.5cm]{m310.ps}
+\includegraphics[width=2.5cm]{m301.ps}
-\frac{\beta_0}{\e} \includegraphics[width=2.5cm]{m300.ps}
\end{eqnarray}
\vspace{0.5cm}

\noindent 
$\caqqg{2}$ is proportional to $\CF$, while $\caqqg{4}$ contains $\CF^2$, 
$\CF\CA$ and $\CF T_R N_F$ terms~\cite{Ellis:1981wv,Fabricius:1981sx}.

The virtual unresolved contributions to two-jet final states at this order arise
from the one-loop corrections to $\gamma^*\to q \bar q g$ in the  limit where the
gluon becomes either collinear or soft. In both these  limits, one observes a
factorization of the {\it renormalized}  matrix element. In the collinear limit,
the one-loop $n+1$ parton  amplitude factorises as a one-loop $n$ parton amplitude
multiplied by a tree-level splitting function plus a tree-level $n$ parton
amplitude multiplied by a one-loop splitting 
function~\cite{Bern:1994zx,Bern:1998sc,Kosower:1999xi,Bern:1999ry,Kosower:1999rx}. 
In the soft limit, one observes the factorization  of the one-loop $n+1$ parton
amplitude into the one-loop $n$ parton amplitude  times the leading order soft
factor, plus the  $n$ parton tree level amplitude multiplied with the
soft gluon factor  at the one-loop
level~\cite{Bern:1998sc,Bern:1999ry,Catani:2000pi}. This one-loop correction to the
soft gluon factor contains only the colour factor $\CA$. In the infrared limit, we
therefore obtain the usual $\Ione$ factor  multiplying the one-loop two-parton
amplitude plus a new contribution $\dIone$ multiplying the tree-level contribution,
\vspace{-0.5cm} 
\begin{eqnarray} 
\includegraphics[width=2.5cm]{m310.ps} 
&\longrightarrow&  
\includegraphics[width=2.5cm]{m210I1a.ps} +
\includegraphics[width=2.5cm]{m210I1b.ps} \nonumber \\ &+&
\includegraphics[width=2.5cm]{m200dI1a.ps} +
\includegraphics[width=2.5cm]{m200dI1b.ps} +{\rm finite}. 
\end{eqnarray}
\vspace{0.5cm}

\noindent
The first term is proportional to $\CF^2$, while 
$\dIone$ contains terms proportional to $\CF$ multiplied 
by double poles in $\e$
and terms proportional to $\CA$.  Altogether, the infrared structure from the 
one-loop three parton contribution is given by,
\vspace{-0.5cm}
\begin{eqnarray}
\label{eq:polesc}
\lefteqn{\includegraphics[width=2.5cm]{m310.ps}
+\includegraphics[width=2.5cm]{m301.ps}
-\frac{\beta_0}{\e} \includegraphics[width=2.5cm]{m300.ps}}\nonumber \\
&=& 
2\Re \Biggl[
-2 \Ione \includegraphics[width=2.5cm]{m210.ps} 
-2 \dIone \includegraphics[width=2.5cm]{m200.ps} \nonumber \\
&& +\frac{\beta_0}{\e} \Ione \includegraphics[width=2.5cm]{m200.ps}\Biggr ]
+{\rm finite}.
\end{eqnarray}

\noindent
When Eq.~(\ref{eq:polesc}) is 
taken together with Eqs.~(\ref{eq:polesa}) and (\ref{eq:polesb}) we see
that the terms of the type 
\begin{displaymath}
\Ione \includegraphics[width=2.5cm]{m210.ps}
\end{displaymath} 
and
\begin{displaymath}
\frac{\beta_0}{\e} \Ione \includegraphics[width=2.5cm]{m200.ps}
\end{displaymath} 
precisely cancel.

\vspace{0.2cm}

\noindent
{\it The tree-level four parton contributions}

\vspace{0.15cm}

\noindent
The renormalized $\gamma^* \to q\bar q +2$~parton
amplitudes can be written as
\begin{eqnarray}
|\cmqqgg\rangle &=& \sqrt{4\pi\alpha}e_q  \asmuopi  \left[
|\cmqqgg^{(0)}\rangle 
+ {\cal O}(\alpha_s^2) \right] , \nonumber \\
|\cmqqqq\rangle &=& \sqrt{4\pi\alpha}e_q  \asmuopi  \left[
|\cmqqqq^{(0)}\rangle 
+ {\cal O}(\alpha_s^2) \right] , 
\end{eqnarray}
where the $|{\cal M}_{q\bar qpp}^{(i)}\rangle$ are the $i$-loop contributions to the 
renormalized amplitude for $p=g,~q$.   

The squared four-parton amplitudes, summed over spins, 
colours and quark flavours, 
are denoted by $\cmqqgg$ and $\cmqqqq$, which have the perturbative expansions,
\begin{equation}
{\cal A}_{q\bar qpp} = 4\pi\alpha\sum_q N e_q^2 \Bigg[ 
\asmuopi^2 {\cal A}_{q\bar qpp}^{(4)} + {\cal O}(\alpha_s^3)\Bigg] \;,
\end{equation}
where,\vspace{-0.5cm}
\begin{eqnarray}
\label{eq:Aqqgg4}
\caqqgg{4} &=& \langle\cmqqgg^{(0)}|\cmqqgg^{(0)}\rangle 
=\includegraphics[width=2.5cm]{m400gg.ps},\\
&& \nonumber \\
\label{eq:Aqqqq4}
\caqqqq{4} &=& \langle\cmqqqq^{(0)}|\cmqqqq^{(0)}\rangle 
=\includegraphics[width=2.5cm]{m400qq.ps}.
\end{eqnarray}
\vspace{0.5cm}

\noindent The infrared singularity structure of these matrix elements is well
documented in the double soft, triple collinear and independent double 
collinear
limits~\cite{Gehrmann-DeRidder:1998gf,Campbell:1998hg,Catani:1998nv,Catani:1999ss,DelDuca:1999ha}.  
At present, it is not clear how to efficiently isolate these singularities. Up
to now, the only calculation involving these elements  was carried out for the
case of  photon-plus-one-jet final  states in electron--positron
annihilation~\cite{Gehrmann-DeRidder:1998gf,Gehrmann-DeRidder:1997wx}, which
was performed in the  hybrid subtraction method, which is an extension of the 
phase space slicing procedure of~\cite{Giele:1992vf}.  These double unresolved
factors have also been used to compute the general form of 
logarithmically-enhanced contributions (up to next-to-next-to-leading
logarithmic accuracy) of the transverse momentum distributions of high-mass
systems in hadronic collisions~\cite{deFlorian:2000pr,deFlorian:2001zd}.


\vspace{0.2cm}

\noindent
{\it Summary}

\vspace{0.15cm}

\noindent
The infrared singular structure of the two-loop amplitudes is well
documented.   Concerning the singular structure of the 
one-loop three parton amplitudes, it
seems straightforward to evaluate and to construct counter terms.  There is a
precise cancellation with some of the singular terms in the one- and two-loop
two parton contributions.   However, more work is required to disentangle the
singular structure of the tree-level  four-parton contributions.

\subsection{The high energy limit of qcd at two loops}
\label{sec:ho;highenergy}

In this section we give a summary of the present status of
the analytic structure of QCD amplitudes in the limit of forward
and backward scattering. The gluon and quark Regge trajectories are
evaluated at two loop accuracy.

\subsubsection{Forward Scattering}

The scattering of two particles, with momenta $p_1 p_2\to p_3 p_4$
in their centre-of-mass frame, is parametrised by the invariants
$s=(p_1+p_2)^2$, $t=(p_2-p_3)^2$ and $u=(p_2-p_4)^2$. Since 
$t = - s (1-\cos\theta)/2$, where $\theta$ is
the scattering angle, the kinematic region where
the squared centre-of-mass energy is much greater than
the momentum transfer, $s\gg|t|$, defines the forward scattering.
It has been known for a long time that in the $s\gg|t|$ limit
the scattering amplitudes of a gauge
theory take a simple analytic form. In fact, in that limit a generic
scattering process is dominated by the exchange of the particle of
highest spin in the $t$ channel. That is a gluon in QCD
(which we take as representative of $SU(N)$
Yang-Mills theories), a photon in QED, and a graviton in quantum gravity.

In the scattering of two partons in QCD in the limit $s\gg|t|$,
the processes which yield the largest contribution are
quark-quark, quark-gluon and gluon-gluon scattering. The other
processes, like $q\ \bar{q}\to g\ g$, which do not feature gluon
exchange in the $t$ channel, are subleading in this limit.
The leading processes all have the same analytic form.
In fact, the tree amplitude for parton-parton scattering 
$i_a(p_2)\, j_b(p_1)\to i_{a'}(p_3)\,j_{b'}(p_4)$, with $i, j$ either a
quark or a gluon, may be written as~\cite{Kuraev:1976ge}
%
\begin{equation}
\cM^{(0) aa'bb'}_{ij\to ij} = 2  s
\left[\gs\, (T^c_r)_{aa'}\, C^{i(0)}(p_2,p_3) \right]
{1\over t} \left[\gs\, (T^c_r)_{bb'}\, C^{j(0)}(p_1,p_4) 
\right]\, ,\label{elas}
\end{equation}
%
where $a, a', b, b'$ label the parton colours, and $r$ represents either
the fundamental $(F)$ or the adjoint $(G)$ representations of $SU(N)$, with
$(T^c_G)_{ab} = i f^{acb}$ and ${\rm tr}(T^c_F T^d_F) = \delta^{cd}/2$.
The coefficient function $C^{i(0)}$ is process dependent; in helicity
space it just contributes a phase 
factor~\cite{DelDuca:1995zy,DelDuca:1999ha}: its square is 1.
By removing the colour factor and replacing the
strong with the electromagnetic coupling, \eqn{elas} holds also for
the forward scattering of two charged leptons in QED with exchange 
of a photon in the $t$ channel.

The amplitude (\ref{elas}) constitutes the leading term of the BFKL 
theory~\cite{Kuraev:1976ge,Kuraev:1977fs,Balitsky:1978ic}, which models
strong-interaction processes with two large and disparate scales,
by resumming the radiative corrections to parton-parton
scattering. This is achieved to leading logarithmic (LL) accuracy, in
$\ln(s/|t|)$, through the BFKL equation, {\em i.e.}
an integral equation obtained by computing the one-loop LL
corrections to the gluon exchange in the $t$ channel. These corrections are:
the emission of a gluon along the ladder~\cite{Lipatov:1976zz}, and 
the one-loop {\em gluon Regge trajectory} (see \eqn{exp1loop}).
To see how the latter comes about, we generalise
\eqn{elas} to include the virtual radiative corrections
and write the high-energy amplitude for parton-parton 
scattering, with exchange of a colour octet in the $t$
channel, as~\cite{Fadin:1993wh,Lipatov:1989bs}
%
\begin{eqnarray}
\lefteqn{\cM^{aa'bb'}_{ij\to ij} } \nonumber\\ 
&=& s
\left[\gs\, (T^c_{r_A})_{aa'}\, C^{i}(p_2,p_3) \right]
{1\over t} \left[\left({-s\over -t}\right)^{\alpha(t)} +
\left({s\over -t}\right)^{\alpha(t)}  \right]
\left[\gs\, (T^c_{r_A})_{bb'}\, C^{j}(p_1,p_4) \right]
\nonumber\\ &+& s
\left[\gs\, (T^c_{r_S})_{aa'}\, C^{i}(p_2,p_3) \right]
{1\over t} \left[\left({-s\over -t}\right)^{\alpha(t)} -
\left({s\over -t}\right)^{\alpha(t)}  \right] 
\left[\gs\, (T^c_{r_S})_{bb'}\, C^{j}(p_1,p_4) \right]\, ,\label{elasb}
\end{eqnarray}
%
with 
%
\begin{equation}
(T^c_{G_A})_{aa'} = i f^{aca'} \quad (T^c_{G_S})_{aa'} = d^{aca'}
\quad T^c_{F_A} = T^c_F \quad T^c_{F_S} = \sqrt{N_c^2-4\over N_c^2} T^c_F\, .
\end{equation}
%
The first (second) line of \eqn{elasb} corresponds to the exchange of a
colour octet in the $t$ channel, of negative (positive) signature under 
$s\leftrightarrow u$ exchange.
The function $\alpha(t)$ in \eqn{elasb} is the gluon Regge trajectory. 
It has the perturbative expansion
%
\begin{equation}
\alpha(t) = \tilde\gs^2(t) \alpha^{(1)} + 
\tgs^4(t) \alpha^{(2)} + \ord(\tgs^6)\,
,\label{alphb}
\end{equation}
% \begin{equation}
% \alpha(t) = \tilde{g_{\scriptscriptstyle S}}^2(t) \alpha^{(1)} + 
% \tilde{g_{\scriptscriptstyle S}}^4(t) \alpha^{(2)} + {\cal O}(\tgs^6)\,
% ,\label{alphb}
% \end{equation}
%
while the coefficient function can be written as
%
\begin{equation}
C^{i} = C^{i (0)}(1 + \tgs^2(t) C^{i (1)} + \tgs^4(t) C^{i (2)}) + 
\ord(\tgs^6)\, .\label{fullv}
\end{equation}
%
In \eqns{alphb}{fullv}, we rescaled the coupling
%
\begin{equation}
\tgs^2(t) = \gs^2 \left({\mu^2\over -t}\right)^{\eps}
{1\over (4\pi)^{2-\epsilon}}\, {\Gamma(1+\epsilon)\,
\Gamma^2(1-\epsilon)\over \Gamma(1-2\epsilon)}\, ,\label{rescal}
\end{equation}
%
and used dimensional regularisation in $d=4-2\epsilon$ dimensions.
Note that in multiplying \eqn{elasb} by the tree amplitude, the second
line of \eqn{elasb} contributes only to quark-quark scattering,
since the colour factor $f^{ada'}$, occurring in any tree amplitude 
with an external gluon, acts as an $s$-channel projector.
Then we write the projection of the amplitude (\ref{elasb}) on the tree
amplitude as an expansion in $\tgs^2(t)$
%
\begin{equation}
\cM^{aa'bb'}_{ij\to ij} \cM^{(0) aa'bb'}_{ij\to ij}
= |\cM^{(0) aa'bb'}_{ij\to ij}|^2 \left( 1 +
\tgs^2(t)\ M^{(1) aa'bb'}_{ij\to ij} + \tgs^4(t) M^{(2) aa'bb'}_{ij\to ij}
+ \ord(\tgs^6) \right)\, ,\label{elasexpand}
\end{equation}
%
with $i, j = g, q$. 
%
\begin{figure}[t]
(a)
\begin{center}
~\includegraphics[width=2.5cm]{A1.ps}
\end{center}
(b)
\begin{center}
~\includegraphics[width=2.5cm]{C1up.ps}\hspace{2cm} 
~\includegraphics[width=2.5cm]{C1down.ps} 
\end{center}
\caption{Schematic one-loop expansion of the factorised form for the 
high energy limit of the parton-parton scattering amplitude. 
The pairs of concentric circles represent the one-loop corrections to the 
coefficient function and Regge trajectory and the individual diagrams
represent  terms that contribute at (a) leading  and (b) next-to-leading
logarithmic order. }
\label{fig:oneloop}
\end{figure}
%
The one-loop coefficient of \eqn{elasexpand} is, 
%
\begin{equation}
M^{(1) aa'bb'}_{ij\to ij}
= \alpha^{(1)} \ln\left({s\over -t}\right) +\ C^{i(1)} + C^{j(1)}
- i{\pi\over 2} \left( 1 + K {N_c^2-4\over N_c^2} \right)
\alpha^{(1)}\, ,\label{exp1loop}
\end{equation}
%
where $K = 1$ for quark-quark scattering, and $K = 0$ in the other cases.
Schematically, this is illustrated in Fig.~\ref{fig:oneloop}.
In \eqn{exp1loop}
we used the usual prescription $\ln(-s) = \ln(s) - i\pi$, for $s > 0$.
The one-loop gluon Regge trajectory, $\alpha^{(1)} = 2\CA/\eps$, 
Fig.~\ref{fig:oneloop}(a), is independent of the type of parton 
undergoing the scattering process (it is {\em universal}).
Conversely, the one-loop coefficient 
functions, $C^{i (1)}$, Fig.~\ref{fig:oneloop}(b),
are process and IR-scheme dependent (see 
Ref.~\cite{DelDuca:2001gu} and references therein).
They can be used to construct the next-to-leading order (NLO)
impact factors, to be used in conjunction with the BFKL resummation
at next-to-leading-log (NLL) accuracy~\cite{Fadin:1998py}.
Note that \eqn{exp1loop} forms a system of three equations (given
by gluon-gluon, quark-quark and quark-gluon scattering)
and only two unknowns: the one-loop coefficients $C^{g(1)}$ and 
$C^{q(1)}$. Thus, we can use two equations to determine 
$C^{g(1)}$ and $C^{q(1)}$, and the third to show that high-energy 
factorisation holds to one-loop accuracy. Finally,
note that in \eqn{exp1loop} the contribution of the positive signature
gluon appears in the imaginary part, thus it can contribute only to
a next-to-next-to-leading order (NNLO) calculation.
%
\begin{figure}[t]
(a)
\begin{center}
~\includegraphics[width=2.5cm]{A1A1.ps}
\end{center}
(b)
\begin{center}
~\includegraphics[width=2.5cm]{A2.ps}\hspace{1cm} 
~\includegraphics[width=2.5cm]{C1upA1.ps}\hspace{1cm} 
~\includegraphics[width=2.5cm]{C1downA1.ps} 
\end{center}
(c)
\begin{center}
~\includegraphics[width=2.5cm]{C2up.ps}\hspace{1cm}  
~\includegraphics[width=2.5cm]{C2down.ps}\hspace{1cm}  
~\includegraphics[width=2.5cm]{C1C1.ps} 
\end{center}
\caption{ Schematic two-loop expansion of the factorised form for the 
high energy limit of the parton-parton scattering amplitude. 
The combinations of ovals and circles represent the one-loop and 
two-loop corrections to the coefficient function and Regge trajectory and
the individual diagrams represent  terms that contribute at 
(a) leading, (b) next-to-leading and (c) next-to-next-to-leading
logarithmic order. }
\label{fig:twoloop}
\end{figure}
%
The two-loop coefficient of \eqn{elasexpand} is,
%
\begin{eqnarray}
M^{(2) aa'bb'}_{ij\to ij}
&=& {1\over 2} \left(\alpha^{(1)}\right)^2 \ln^2\left({s\over -t}\right) 
\nn\\  &+& \left[ \alpha^{(2)} + \left(C^{i(1)} + C^{j(1)}
\right) \alpha^{(1)} - 
i{\pi\over 2} \left( 1 + K {N_c^2-4\over N_c^2} \right)
\left(\alpha^{(1)}\right)^2
\right] \ln\left({s\over -t}\right) \nn\\ &+& 
\left[ C^{i(2)} + C^{j(2)} + 
C^{i(1)}\ C^{j(1)} - 
{\pi^2\over 4} \left( 1 + K {N_c^2-4\over N_c^2} \right)
\left(\alpha^{(1)}\right)^2 \right] \nn\\ 
&-& i{\pi\over 2} \left( 1 + K {N_c^2-4\over N_c^2} \right)
\left[ \alpha^{(2)} + \left(C^{i(1)} + 
C^{j(1)} \right) \alpha^{(1)} \right]
\, .\label{exp2loop}
\end{eqnarray}
%
Schematically, this is illustrated in Fig.~\ref{fig:twoloop}.
The first line of \eqn{exp2loop} is just the exponentiation of the one-loop
trajectory (Fig.~\ref{fig:twoloop}(a)). The second line of \eqn{exp2loop}
allows to determine $\alpha^{(2)}$, the two-loop gluon 
trajectory (the first diagram in Fig.~\ref{fig:twoloop}(b)), which
provides the virtual corrections to the NLL kernel of the BFKL equation. 
The two-loop gluon trajectory is universal, and has been computed in 
Ref.~\cite{Fadin:1995xg,Fadin:1996tb,Fadin:1996km,Blumlein:1998ib}.
Recently it was re-evaluated in a completely independent 
way~\cite{DelDuca:2001gu}, by taking the high energy limit of the two-loop
amplitudes for parton-parton 
scattering~\cite{Anastasiou:2001sv,Anastasiou:2000kg,
Anastasiou:2000ue,Glover:2001af}. 
The (unrenormalised) two-loop gluon Regge trajectory is
%
\begin{equation}
\alpha^{(2)} = \CA \left[ \beta_0 {2\over \eps^2} + K{2\over \eps} 
+ \CA \left({404\over 27} - 2\zeta_3\right) 
+ \NF \left(-{56\over 27}\right)\right]\, ,\label{eq:2loop}
\end{equation}
%
where
%
\begin{equation} 
\beta_0= {(11\CA-2\NF)\over 6}, \qquad\qquad K = \left({67\over 18} -
{\pi^2\over 6} \right) \CA - {5\over 9} \NF\, .\label{eq:beta}
\end{equation}
%
Note that at two loop accuracy, the gluon trajectory vanishes
in the Abelian limit, $\CA\to 0$. That is in agreement with our 
expectations: in the forward scattering of two charged leptons in QED, a log
of the type $\ln(s/t)$ occurs for the first time at the three loop
level. Thus we expect that only in $\alpha^{(3)}$ it occurs a term
which is independent of $\CA$.
However, note that in QED there is no analog of \eqn{elasb}, because
the photon does not reggeise~\cite{mandel}.

The third and fourth lines of \eqn{exp2loop} are respectively
the real (Fig.~\ref{fig:twoloop}(c)) and the imaginary parts of the constant 
term. Note that the gluon of
positive signature does not contribute to the BFKL resummation at LL and NLL
accuracy. In principle, the two-loop
coefficient functions $C^{i(2)}$ could be used to construct
the NNLO impact factors, if the BFKL resummation held to 
next-to-next-to-leading-log (NNLL) accuracy.
The real constant term of \eqn{exp2loop}
can be used to evaluate the two-loop coefficient functions $C^{g(2)}$
and $C^{q(2)}$. Like in \eqn{exp1loop}, we have three equations 
(given by the two-loop amplitudes for gluon-gluon, 
quark-quark and quark-gluon scattering~\cite{Anastasiou:2001sv,
Anastasiou:2000kg,Anastasiou:2000ue,Glover:2001af})
and only two unknowns. Thus we can
use two equations to evaluate $C^{g(2)}$ and $C^{q(2)}$, and the third
to check if high-energy factorisation holds at NNLO accuracy.
In fact, it was found that factorisation is 
violated~\cite{DelDuca:2001gu}. That, in turn, voids the evaluation 
of $C^{g(2)}$ and $C^{q(2)}$. The source of the violation is not clear 
at present. It might be due to the (yet unknown) exchange
of three or more reggeised gluons in \eqn{elasb}.

\subsubsection{Backward scattering}

In the scattering of two particles in their centre-of-mass frame,
$u = - s (1+\cos\theta)/2$, thus the
kinematic region where $s\gg |u|$ defines the backward scattering.
In this limit, we consider a scattering process with matter exchange, 
namely quark-gluon scattering, which proceeds via
the exchange of a quark in the crossed channel\footnote{By crossing symmetry,
this is equivalent to quark pair annihilation (or creation from) two
gluons at small angles.}. In this limit, the tree amplitude for quark-gluon 
scattering $q_a(p_2)\ g_b(p_1)\to q_{a'}(p_3)\ g_{b'}(p_4)$ can be written as
%
\begin{equation}
{\cM}^{(0)}_{qg\to gq} = -2i \left[\gs\, (T^b)_{a'i} C^{(0)}_{gq}(p_1,p_3) 
\right] \sqrt{s\over - u} \left[\gs\, (T^{b'})_{ia}\, C^{(0)}_{qg}(p_2,p_4) 
\right]\, ,\label{elasq}
\end{equation}
%
where $a, a'$ ($b, b'$) label the colour of the quarks (gluons). In
helicity space, the coefficient functions $C^{(0)}$ just contribute a
phase factor~\cite{Bogdan:2002sr}. 

Generalising \eqn{elasq} to include virtual radiative corrections,
it was shown that the quark reggeises as well.
In fact, in the limit $s\gg |u|$ the amplitude for quark-gluon
scattering with exchange in the $u$ channel of a colour triplet, with positive 
signature under $s\leftrightarrow t$ exchange, is~\cite{Fadin:1977jr}
%
\begin{equation}
{\cM}_{qg\to gq} = -i \left[\gs\, (T^b)_{a'i} C_{gq}(p_1,p_3) 
\right] \sqrt{s\over - u} \left[ \left({s\over - u}\right)^{\delta(u)}
+ \left({- s\over - u}\right)^{\delta(u)} \right]
\left[\gs\, (T^{b'})_{ia}\, C_{qg}(p_2,p_4) \right]\, .\label{allooppm}
\end{equation}
%
The function $\delta(u)$ is the {\em quark Regge trajectory}.
It has the perturbative expansion
%
\begin{equation}
\delta(u) = \tilde\gs^2(u) \delta^{(1)} + 
\tgs^4(u) \delta^{(2)} + \ord(\tgs^6)\, ,\label{delta}
\end{equation}
%
with $\tilde\gs^2$ given in \eqn{rescal}. 
The coefficient function $C$ has the same expansion as $C^i$ in
\eqn{fullv}, but for replacing $\tilde\gs^2(t)$ with $\tilde\gs^2(u)$,
and like in \eqn{elasexpand}, we can write the interference of \eqn{allooppm}
with the tree amplitude (\ref{elasq}) as an expansion in $\tilde\gs^2(u)$
%
\begin{equation}
\cM \cM^{(0) *} = |\cM^{(0)}|^2
\left( 1 + \tgs^2(u)\ M^{(1)} + \tgs^4(u) M^{(2)}
+ \ord(\tgs^6) \right)\, .\label{elasexpandpm}
\end{equation}
%
The one-loop coefficient of \eqn{elasexpandpm} is
%
\begin{equation}
M^{(1)}
= \delta^{(1)} \ln\left({s\over -u}\right) +\ 2 C^{(1)}
- i{\pi\over 2} \delta^{(1)}\, .\label{exp1loopplus}
\end{equation}
%
The one-loop quark trajectory is 
$\delta^{(1)} = 2\CF/\eps$~\cite{Fadin:1977jr}. Note that
$\delta^{(1)}$ does not vanish in the Abelian limit $\CF\to 1$.
In fact, in electron-photon scattering in QED,
the electron reggeises~\cite{gellmann,McCoy:1976fj}.
The coefficient function $C^{(1)}$
can be found in Ref.~\cite{Bogdan:2002sr,Fadin:2001dc}.
The two-loop coefficient of \eqn{elasexpandpm} is
%
\begin{eqnarray}
M^{(2)}
&=& {1\over 2} \left(\delta^{(1)}\right)^2 \ln^2\left({s\over -u}\right) 
+ \left[ \delta^{(2)} + 2 C^{(1)} \delta^{(1)} - 
i{\pi\over 2} \left(\delta^{(1)}\right)^2
\right] \ln\left({s\over -u}\right) \nn\\ &+& 
\left[ 2 C^{(2)} + \left( C^{(1)}\right)^2 -
{\pi^2\over 4} \left(\delta^{(1)}\right)^2  \right]
- i{\pi\over 2} \left[ \delta^{(2)} + 2 C^{(1)} \delta^{(1)} \right]
\, .\label{exp2loopplus}
\end{eqnarray}
%
The first term is just the exponentiation of the one-loop trajectory. 
Using the two-loop amplitude for quark-gluon 
scattering~\cite{Anastasiou:2001sv}, the second term of \eqn{exp2loopplus}
allows to determine $\delta^{(2)}$, the (unrenormalised)
two-loop quark Regge trajectory~\cite{Bogdan:2002sr}
%
\begin{equation}
\delta^{(2)} = \CF \left[ \beta_0 {2\over \eps^2} + K {2
\over \eps} + \CA \left({404\over 27} -2 \zeta_3\right)
+ \NF \left(-{56\over 27}\right) 
+ (\CF-\CA) \left(16\zeta_3\right) \right]\, ,\label{twolooptraj}
\end{equation}
%
with $\beta_0$ and $K$ as in \eqn{eq:beta}.
Note that \eqn{twolooptraj} has the remarkable feature that by
mapping $\CF\to\CA$, the two-loop gluon trajectory (\ref{eq:2loop})
is obtained. Since the forward and backward scattering are
seemingly unrelated, there is presently no understanding of why that
occurs.

The second line of \eqn{exp2loopplus} displays respectively the real
and the imaginary parts of the constant term.
If \eqn{allooppm} held at NNLO, by comparing the
two-loop amplitude for quark-gluon scattering~\cite{Anastasiou:2001sv} 
with the real part of the constant term of \eqn{exp2loopplus}
one could derive the two-loop coefficient function $C^{(2)}$.
However, \eqn{allooppm} is not expected to hold at NNLO, because
of possible unknown Regge cut contributions.

\subsubsection{Summary}

In conclusion, we have given a summary of the present status of the analytic
structure of QCD amplitudes in the limit of forward and backward scattering. We
have displayed the gluon and quark Regge trajectories at two loop accuracy.
They are strikingly similar: the gluon Regge trajectory can be obtained from
the quark  trajectory by mapping $\CF\to\CA$. 


\subsection{CONCLUSION AND OUTLOOK}
\label{sec:ho;conc}

In this working group we addressed issues related to multiparticle states at
NLO and precision observables at NNLO. The motivation for this is very
straightforward: we need more accurate theoretical predictions in order to make
a sensible comparison with experiment, either to establish the existence of new
physics or to make more precise determinations of the parameters of the theory.
During the course of the workshop there
has been significant progress in this field.   In particular, work
initiated at the workshop have yielded the following identifiable results,
\begin{itemize}
\item 
The first Feynman diagrammatic calculation of a hexagon amplitude~\cite{Binoth:2001vm}. This was
performed   in the Yukawa model where all external legs are  massless scalars
attached to a massless fermion loop. This is a good model to study how the
enormous cancellations in the scalar integral reductions occur without the
complications of tensor structures in the numerator. The next step is to  apply
these methods to  realistic examples  including gauge bosons and a nontrivial
infrared structure.  However, it is reasonable to expect that the recombination
of scalar integrals will work similarly, such that efficient algorithms to
calculate  six-point amplitudes at one loop are in reach.  
\item 
The first evaluation of two-loop four point amplitudes with one off-shell
leg~\cite{Garland:2001tf} for the case $\gamma^* \to q\bar q g$.  This process 
is important in determining the strong coupling constant from hadronic data in
electron-positron annihilation.  By analytical continuation it should be
possible to extend these results to describe processes in hadron-hadron
collisions, $p\bar p \to V + {\rm jet}$,  and deep inelastic scattering,
$\gamma^* p \to (2+1)$~jets. 
\item 
The first determination of the two-loop quark Regge
trajectory~\cite{Bogdan:2002sr} together with confirmation of the two-loop
gluon Regge trajectory~\cite{DelDuca:2001gu}.   These trajectories control the
behaviour of the scattering amplitudes in
the high energy limit - or equivalently forward or backward scattering.
\end{itemize}

While the benefits of NLO calculations are generally well appreciated,  there
is some debate as to the motivation for NNLO calculations.   In 
Sec.~\ref{sec:ho;nnlorationale} we give reasons why NNLO
predictions should qualitatively improve the agreement between theory and
experiment, reduced renormalisation and factorisation scale dependence, 
reduced phenomenological power correction, better matching between the
parton-level and hadron-level jet algorithm and a better description of the
transverse momentum of the final state.

The basic ingredients for computing NNLO corrections are  by now either known
or conceptually understood.  However a method of how to combine the
individually infrared singular parts in a way that allows the construction of
general purpose NNLO parton level Monte Carlo programs has not been
established. This is clearly the next major task and requires considerable
effort.   During the workshop some preparatory steps were made in this 
direction for the relatively simple case of $\gamma^* \to q\bar q$ and the
general structure of the infrared singularities for this process is discussed
in Sec.~\ref{sec:ho;exercise}.  We expect that by the time of the next workshop,
there will be significant progress to report.

\vskip1cm
\noindent

\subsection*{Acknowledgements}

We would like to thank the organizers of ``Les Houches 2001" for creating a
stimulating and productive environment for the workshop. We would like to thank
our many and various collaborators for many stimulating and  thought provoking
discussions; T.B., J.-Ph.G. and G.H. thank C.~Schubert, V.D.D and E.W.N.G.
thank V.~Fadin and A.~Bogdan and    T.G. and E.W.N.G. thank L.~Garland,
A.~Koukoutsakis and E.~Remiddi. 

%\bibliography{bibho}
%\end{document}

