%Title: Fermionic Corrections to the Three-Loop Matching Coefficient of the Vector Current
%Author: P. Marquard, J.H. Piclum, D. Seidel, M. Steinhauser
%Published: <EM> Nuclear Physics </EM> <B> B758 </B> (2006) 144-160. 
%hep-ph/0607168
%- {{{ Header:

\documentclass[12pt]{article}
\usepackage{a4wide,epsfig,amsmath,amssymb,cite}
\usepackage{scalefnt}
\voffset0cm
\hoffset0cm
\oddsidemargin0cm
\evensidemargin0cm
\topmargin0cm
\textwidth16.cm
\textheight22cm

\renewcommand{\textfraction}{0}
\renewcommand{\topfraction}{1}
\renewcommand{\bottomfraction}{1}

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

\newcommand{\note}[1]{{\tiny (note)}\marginpar {\scriptsize #1}}

%%%\newcommand{\qsla}{q\!\!\!/\,\,}
\newcommand{\Lmu}{L_\mu}
\newcommand{\lmu}[1]{L_{#1}}
\newcommand{\order}[1]{\mathcal{O}\left( #1 \right)}

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

\begin{document}

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

\title{\vskip-3cm{\baselineskip14pt
\begin{flushleft}
\normalsize DESY 06-107 \\
\normalsize SFB/CPP-06-31 \\
\normalsize TTP06-22 \\
\end{flushleft}}
\vskip1.5cm
Fermionic Corrections to the Three-Loop Matching Coefficient of the
Vector Current
}

\author{\small 
  P. Marquard$^{(a)}$, J.H. Piclum$^{(a,b)}$, D. Seidel$^{(a)}$,
  M. Steinhauser$^{(a)}$
\\
{\small\it (a) Institut f{\"u}r Theoretische Teilchenphysik,
  Universit{\"a}t Karlsruhe (TH)}\\
{\small\it 76128 Karlsruhe, Germany}
\\
{\small\it (b) II. Institut f\"ur Theoretische Physik, 
  Universit\"at Hamburg}\\
{\small\it 22761 Hamburg, Germany}
}

\date{}

\maketitle

\thispagestyle{empty}

\begin{abstract}
In this paper we consider the matching coefficient of the vector
current between Quantum Chromodynamics (QCD) and Non-Relativistic QCD
(NRQCD) to three-loop order in perturbation theory.
We evaluate the fermionic corrections containing a closed massless
fermion loop.
The results are building blocks both for the bottom and top quark system
at threshold. We explain in detail the methods used for the evaluation
of the Feynman diagrams, classify the occurring master integrals and
provide results for the latter. The numerical effects are significant. 
They have the tendency to improve the behaviour of the perturbative
series --- both for the bottom and top quark system.

\medskip

\noindent
PACS numbers: 12.38.Bx 14.65.Ha
\end{abstract}

\newpage

%- }}}
%- {{{ Introduction:

\section{Introduction}

One of the main goals of a future international linear collider (ILC) is
the precise measurement of the top quark threshold. 
Next to a precise extraction of the strong coupling an unrivaled
determination of the top quark mass and its width is possible.
This would open up a new chapter in the electroweak precision physics
which leads to very strong checks of the Standard Model or
possible extensions.


The theoretical calculations are based on an effective
theory~\cite{Caswell:1985ui,Bodwin:1994jh} (for a review
see~\cite{Brambilla:2004wf})
which enables the simultaneous expansion in the two small
parameters present in threshold phenomena, the strong coupling
$\alpha_s$ and the velocity of the heavy quarks, $v$. 
Furthermore, a resummation of correction terms like $(\alpha_s/v)^n$
to all orders is possible in perturbation theory.
In a first step
the effective theory is constructed in a matching procedure where it
is required that the Green functions agree with the corresponding ones
in full QCD in the limit of infinitely heavy top quark mass.
In a second step the Green function in the effective theory has to be
evaluated leading to the desired cross section.

Experimental studies~\cite{Martinez:2002st} have shown that an
uncertainty of less 
than 3\% is necessary for the theoretical predictions for the cross
section of the reaction $e^+ e^-\to t\bar{t}$ in the threshold region,
i.e. for $\sqrt{s}\approx 350$~GeV.
This has not yet been achieved.

As far as the strong coupling is concerned
currently the analysis is complete to next-to-next-to-leading
order (NNLO)~\cite{Hoang:2000yr} demonstrating the urgent need for the 
third-order corrections.
This huge enterprise has been started already some years ago and many 
building blocks have been evaluated starting from the determination of
the third-order Hamiltonian~\cite{Kniehl:2002br}, third-order
corrections to the energy levels and wave
functions~\cite{Penin:2002zv,Penin:2005eu,Beneke:2005hg} up to the
evaluation of the contributions from the Coulomb potential to the
non-relativistic Green function~\cite{Beneke:2005hg}.
Next to fixed-order corrections also the resummation of
next-to-next-to leading logarithmic terms has been considered and quite
some progress has been
achieved~\cite{Pineda:2001ra,Pineda:2001et,Hoang:2003ns,Penin:2004ay}.

In order to reach a theory-uncertainty below 3\% it is also necessary to
worry about electroweak effects which can be parametrically 
as big as the third-order QCD terms. First results in this context
are available~\cite{Hoang:2004tg,Eiras:2006xm} which show that indeed numerical
effects in the region of a few percent are possible.

In this paper a further building block needed for the N$^3$LO
analysis is provided: the fermionic corrections to the matching
coefficient of the vector current. More precisely we compute the
three-loop contributions to the matching coefficient containing at
least one closed loop of massless fermions which we enumerate by $n_l$.

The paper is organized as follows: in the next section we describe in
detail the theoretical framework and the methods used for the
calculation. In particular, we classify the occurring integrals and
discuss their evaluation. Sections~\ref{sec::z2} and~\ref{sec::match} 
contain a detailed discussion about the 
wave function and matching coefficient up to order 
$\alpha_s^3 n_l$, respectively. 
In Section~\ref{sec::phen} we briefly discuss the
phenomenological impact of our result and Section~\ref{sec::concl} contains
our conclusions and outlook.
In the Appendix a complete list of master integrals is provided ---
both for the ``natural'' and the $\epsilon$-finite basis.

%- }}}
%- {{{ Method:

\section{Method}

The matching coefficient establishing the relation between the quark
current in the full and effective theory
constitutes a building block for all threshold phenomena involving the
coupling of a photon or $Z$ boson to heavy quarks.
We define the currents through
\begin{eqnarray}
  j_v^\mu &=& \bar{Q} \gamma^\mu Q\,,
  \nonumber\\
  \tilde{j}^i &=& \phi^\dagger \sigma^i \chi\,,
\end{eqnarray}
and generically denote the heavy quark by $Q$.
$\phi$ and $\chi$ are two-component Pauli spinors for quark and
anti-quark, respectively.
The definition of the matching coefficient is established through
\begin{eqnarray}
  j^k_v &=& c_v(\mu) \tilde{j}^k 
  +\frac{d_v(\mu)}{6m_Q^2}\phi^\dagger\sigma^k \vec{D}{}^2\chi
  + \ldots
  \label{eq::def_of_cv}
  \,,
\end{eqnarray}
where $k=1,2,3$ denotes the spacial components. 
$\vec{D}$ contains the space-like
components of the gauge-covariant derivative involving the gluon fields and
the ellipsis stand for operators of higher mass dimension. 
The second term on the r.h.s. of Eq.~(\ref{eq::def_of_cv}) is already of
NNLO with $d_v=1+{\cal O}(\alpha_s)$.
The Wilson coefficients $c_v$ and $d_v$ may be expressed as series
in $\alpha_s$ and represent the contributions from the hard modes
which have been integrated out in order to arrive at the effective theory.
Note that both the higher order terms in the inverse heavy quark
mass and the time-like component is not of interest for the aim of this paper.
The purpose of this paper is the evaluation of the fermionic
correction to the non-singlet contribution of $c_v$ containing at least
one closed massless quark loop.
The one-loop corrections to $d_v$ contain no contribution proportional
to $n_l$.

In order to compute $c_v$ it is obvious to consider the 
$Q\bar{Q}\gamma$ vertex which in the following is denoted by $\Gamma_v$.
This effectively transforms Eq.~(\ref{eq::def_of_cv}) into
\begin{eqnarray}
  Z_2 \Gamma_v &=& c_v \tilde{Z}_2 \tilde{Z}_v^{-1} \tilde{\Gamma}_v
  + \ldots
  \,,
  \label{eq::def_cv}
\end{eqnarray}
where we have on the left- and right-hand side quantities of the full and
effective theory, respectively,
and the ellipses denote terms suppressed by inverse powers
of the heavy quark mass. 
In Eq.~(\ref{eq::def_cv}) it is understood that 
the couplings and masses in $\Gamma_v$ are renormalized.

Since $c_v$ takes into account the degrees of freedom which have been
integrated out from the full theory it depends in our case only 
on $\alpha_s$ and --- for dimensional reasons --- on $\ln\mu^2/m_Q^2$.
In particular, $c_v$ does not explicitly depend on the momenta of the external
particles. Thus, it is useful to apply the so-called threshold
expansion~\cite{Beneke:1997zp,Smirnov:2002pj} to Eq.~(\ref{eq::def_cv})
which has the consequence that $\Gamma_v$ has to be evaluated for 
$s=4m_Q^2$ since all except the hard region cancel in Eq.~(\ref{eq::def_cv}).
Furthermore, on the right-hand side only tree contributions have to be
considered. In particular we have $\tilde{Z}_2=1$.

Starting from two-loop order the matching procedure
exhibits infra-red divergences which are compensated by ultra-violet
divergences of the effective theory rendering physical quantities
finite. In Eq.~(\ref{eq::def_cv}) the renormalization constant
$\tilde{Z}_v$ which generates the anomalous dimension of $\tilde{j}_v$
takes over this part. 
Note that the vector current in the full theory does not get
renormalized. 

\begin{figure}
  \leavevmode
  \epsfxsize=\textwidth
 \epsffile{dias.eps}
\caption{\label{fig::sample} Feynman diagrams contributing to the 
  matching coefficient. Bold lines denote heavy quarks with mass $m_Q$,
  thin lines denote massless quarks and curly lines denote gluons. In
  (d) and (h) mass counterterm diagrams are shown.}
\end{figure}

Sample Feynman diagrams contributing to the one-, two- and three-loop 
fermionic part of $c_v$
are shown in Fig.~\ref{fig::sample}. Due to the special kinematic
situation with $s=q^2=4 m_Q^2$ and on-shell heavy quark lines 
it is possible to perform a partial fractioning in the integrands
corresponding to the various types of diagrams and map all occurring
integrals to one of the following functions:
\begin{eqnarray}
  J^{(2)}_\pm(n_1,\dots,n_5) &=& \left(
  \frac{\mu^{2\epsilon}}{i\pi^{d/2}} \right)^2 \int \frac{{\rm d}^d k\;
  {\rm d}^d l}{ (k^2)^{n_1} (l^2)^{n_2} ((k-l)^2)^{n_3} (k^2 + q\cdot
  k)^{n_4} (l^2 \pm q\cdot l)^{n_5}}\nonumber
  \,,
  \\
  L^{(2)}_\pm(n_1,\dots,n_5) &=& \nonumber \\
  \lefteqn{\!\!\!\!\!\! \left(
  \frac{\mu^{2\epsilon}}{i\pi^{d/2}} \right)^2 \int \frac{{\rm d}^d k\;
  {\rm d}^d l}{ (k^2)^{n_1} (l^2)^{n_2} ((k+l)^2 + q\cdot (k+l))^{n_3}
  (k^2 + q\cdot k)^{n_4} (l^2 \pm q\cdot l)^{n_5}} \,,} &&\nonumber
  \\
  J^{(3)}_\pm(n_1,\dots,n_9) &=& \left(
  \frac{\mu^{2\epsilon}}{i\pi^{d/2}} \right)^3 \int \frac{{\rm d}^d k\;
  {\rm d}^d l\; {\rm d}^d p}{ (k^2)^{n_1} (l^2)^{n_2} (p^2)^{n_3}
  ((k-l)^2)^{n_4} ((l-p)^2)^{n_5} ((p-k)^2)^{n_6}}\nonumber \\
  && \times \frac{(l^2 + q\cdot l)^{-n_8}}{ (k^2 + q\cdot k)^{n_7} (p^2
  \pm q\cdot p)^{n_9}}\nonumber
  \,,
  \\
  L^{(3,n_l)}_\pm(n_1,\dots,n_9) &=& \left(
  \frac{\mu^{2\epsilon}}{i\pi^{d/2}} \right)^3 \int \frac{{\rm d}^d k\;
  {\rm d}^d l\; {\rm d}^d p}{ (k^2)^{n_1} (l^2)^{n_2} ((k+l)^2 + q\cdot
  (k+l))^{n_3} (k^2 + q\cdot k)^{n_4}} \nonumber \\
  && \times \frac{(p^2 + q\cdot p)^{-n_9}}{ (l^2 \pm q\cdot l)^{n_5}
  (p^2)^{n_6} ((p+k)^2)^{n_7} ((p-l)^2)^{n_8}}
  \,.
  \label{eq::integrals}
\end{eqnarray}
For convenience we have also listed the two-loop functions
$J^{(2)}_\pm$ and $L^{(2)}_\pm$ originally defined in
Ref.~\cite{Beneke:1997zp}. The three-loop functions 
$J^{(3)}_\pm$ and $L^{(3,n_l)}_\pm$
contain irreducible
scalar products which are shown as numerators in the three-loop
integrals in Eq.~(\ref{eq::integrals}). The corresponding indices can
only adopt values less or equal to zero. Furthermore, only two out of
the three indices $n_6$, $n_7$ and $n_8$ in $L^{(3,n_l)}_\pm$ can have
positive values.
Note that the integrals $J_+^{(2)}$, $L_+^{(2)}$ and $J_+^{(3)}$ are
actually two-point functions whereas the integrals $J_-^{(2)}$,
$L_-^{(2)}$, $J_-^{(3)}$ and $L^{(3,n_l)}_\pm$ correspond to vertices.

All Feynman diagrams are generated with {\tt
QGRAF}~\cite{Nogueira:1991ex}. The various topologies are identified 
with the help of {\tt q2e} and {\tt
exp}~\cite{Harlander:1997zb,Seidensticker:1999bb} 
which also adapts the notation in order to match the one of
Eq.~(\ref{eq::integrals}).

In a next step the reduction of the various functions to so-called master
integrals has to be achieved. For this step we have chosen two
approaches: the Laporta method~\cite{Laporta:1996mq,Laporta:2001dd} 
and Baikov's method~\cite{Baikov:1996rk,Baikov:1996iu} in
the formulation of~\cite{Smirnov:2003kc}. The application of these
methods reduces the three-loop integrals to twelve master
integrals. Since some of the master integrals are only known numerically
it is very useful to construct a so-called $\epsilon$-finite
basis~\cite{Chetyrkin:2006dh}. 
Details of this procedure and the results for the master integrals are
given in Appendix~\ref{masters}.
Let us in the following provide more details on our implementation of 
each of the two reduction methods.

%- {{{ Laporta:

\subsection{Implementation of the Laporta algorithm}

The Laporta algorithm is based on the
integration-by-parts (IBP) relations~\cite{Chetyrkin:1981qh} where in a first step
a huge system of equations is created by inserting numerical integer values 
into the relations. After assigning a weight specifying the complexity
to each occurring integral the system of equations is solved
step-by-step. In this way an arbitrary integral is expressed in terms
of master integrals which cannot be further reduced.

To date there is only one publicly available computer code, 
{\tt AIR}~\cite{Anastasiou:2004vj}, where the Laporta algorithm is
implemented. We have successfully applied {\tt AIR} to the two-loop
diagrams. However, the three-loop integrals of
Eq.~(\ref{eq::integrals}) cannot be treated with the help of {\tt AIR}. 

For this reason we have chosen to write a new implementation of Laporta's
algorithm, {\tt Crusher}~\cite{PMDS}. It is written in {\tt C++} and uses
{\tt GiNaC}~\cite{Bauer:2000cp} for simple
manipulations like taking derivatives of polynomial quantities. In the
practical implementation of the Laporta algorithm one of the most
time-consuming operations is the simplification of the coefficients
appearing in front of the individual integrals. This task is performed
with the help of {\tt Fermat}~\cite{fermat} where a special interface
has been used (see Ref.~\cite{Tentyukov:2006ys}).
The main features of the
implementation are the automated generation of the IBP identities
and a complete symmetrization of the diagrams. Using {\tt Crusher} we
solved approximately forty million equations to reduce all integrals to master
integrals. 

%- }}}

%- {{{ Baikov:

\subsection{Baikov's method}
Baikov's method is also based on the IBP
relations. Here, however, one explicitly constructs the coefficient
functions of the master integrals as parametric integrals over some
polynomial, which encodes the topology of the initial integral. The
integrations can then be performed as Cauchy integrations around the
origin or as contour integrations between the roots of the
polynomial. The results are rational functions depending on the
dimension and the kinematical invariants.

In some cases, however, it is not possible to perform the integrations
in this way. In these cases one has to solve the
integration-by-parts relations for the parametric integrals. They are
similar to the relations for the initial integral, however, they depend on
less indices and are thus in general significantly simpler. 
Since the relations can still be quite complicated, it is
convenient to again use Laporta's algorithm. In the current calculation
the program {\tt AIR} was used for this task.

In order to perform the calculations, the algorithms for each integral
were implemented in {\tt FORM}~\cite{Vermaseren:2000nd}. In a
first step, however, everything was implemented in {\tt
  Mathematica}. This is very convenient since it is possible
to produce almost all needed {\tt Mathematica} code automatically. Once
the program is finished, it can be easily translated to {\tt FORM} and
the {\tt Mathematica} version can be used for debugging.

The two-loop calculation was done entirely with this method. In the
three-loop case, however, it turned out that a complete solution of the recurrence
relations for the coefficient functions is not possible with {\tt
  AIR}. The calculation was therefore done with the Laporta method
described above. We used Baikov's method only for the $n_l^2$ part and
to cross-check some of the coefficient functions.

%- }}}

%- }}}
%- {{{ Z2:

\section{\label{sec::z2}Wave function renormalization constant}

According to Eq.~(\ref{eq::def_cv}) the wave function renormalization
constant in the on-shell scheme $Z_2$ constitutes a crucial input for
the computation of $c_v$. It has been computed to two-
and three-loop approximation in Refs.~\cite{Broadhurst:1991fy} 
and~\cite{Melnikov:2000zc}, respectively.
We have repeated the calculation of the two-loop and fermionic
three-loop contributions which are needed for the present calculation
and find complete agreement with the literature. For completeness we
repeat the result for $Z_2$ which can be cast into the form
\begin{eqnarray}
  Z_2 &=& 1 + \frac{\alpha_s(\mu)}{\pi} \left(\frac{e^{\gamma_E}}{4 \pi}
  \right)^{-\epsilon} \delta Z_2^{(1)} +
  \left(\frac{\alpha_s(\mu)}{\pi}\right)^2 \left(\frac{e^{\gamma_E}}{4 \pi}
  \right)^{-2\epsilon} \delta Z_2^{(2)} \nonumber \\
  && + \left(\frac{\alpha_s(\mu)}{\pi}\right)^3 \left(\frac{e^{\gamma_E}}{4
  \pi} \right)^{-3\epsilon} \left( \delta Z_2^{(3,nl)} + \mbox{non-$n_l$
  terms} \right) + \mathcal{O}\left(\alpha_s^4\right) \,,
  \label{eq::Z2}
\end{eqnarray}
with
\begin{eqnarray}
  \delta Z_2^{(1)} &=& - C_F \left[ \frac{3}{4\epsilon} + 1 +
  \frac{3}{4} \Lmu + \left( 2 + \frac{1}{16} \pi^2  + \Lmu + \frac{3}{8}
  \Lmu^2 \right) \epsilon \right. \nonumber \\
  && \left. + \left( 4 + \frac{1}{12} \pi^2 - \frac{1}{4} \zeta(3) +
  \left( 2 + \frac{1}{16}\pi^2 \right) \Lmu + \frac{1}{2} \Lmu^2 +
  \frac{1}{8} \Lmu^3 \right) \epsilon^2 \right] \,,
  \label{eq::Z21l}
\end{eqnarray}
%%%%% 2-loop
\begin{eqnarray}
%%% C_A C_F
  \delta Z_2^{(2)} &=& \left[ \frac{11}{32
  \epsilon^2} - \frac{127}{192 \epsilon} - \frac{1705}{384} +
  \frac{5}{16} \pi^2 - \frac{1}{2} \pi^2 \ln 2 + \frac{3}{4} \zeta(3)
  - \frac{215}{96} \Lmu - \frac{11}{32} \Lmu^2 \right. \nonumber \\
  && + \left( -\frac{9907}{768} + \frac{769}{1152} \pi^2 - \frac{23}{8}
  \pi^2 \ln 2 + \pi^2 \ln^2 2 + \frac{129}{16} \zeta(3) - \frac{7}{40}
  \pi^4 + \frac{1}{2} \ln^4 2 + 12 a_4 \right. \nonumber \\
  && \left.\left. + \left( -\frac{2057}{192} +
  \frac{109}{192} \pi^2 - \pi^2 \ln 2 + \frac{3}{2} \zeta(3) \right) \Lmu -
  \frac{259}{96} \Lmu^2 - \frac{11}{32} \Lmu^3\right) \epsilon \right]
  C_AC_F \nonumber \\
%%% C_F^2
  && + \left[ \frac{9}{32 \epsilon^2} + \left( \frac{51}{64} +
  \frac{9}{16} \Lmu \right) \frac{1}{\epsilon} +
  \frac{433}{128} - \frac{49}{64} \pi^2 + \pi^2 \ln 2 -
  \frac{3}{2} \zeta(3) + \frac{51}{32} \Lmu + \frac{9}{16} \Lmu^2
  \right. \nonumber \\
  && + \left( \frac{211}{256} - \frac{339}{128} \pi^2 +
  \frac{23}{4} \pi^2 \ln 2 - 2 \pi^2 \ln^2 2 - \frac{297}{16} \zeta(3) +
  \frac{7}{20} \pi^4 - \ln^4 2 - 24 a_4 \right. \nonumber \\
  && \left.\left. + \left( \frac{433}{64} - \frac{49}{32} \pi^2 + 2
  \pi^2 \ln 2 - 3 \zeta(3) \right) \Lmu + \frac{51}{32} \Lmu^2 +
  \frac{3}{8} \Lmu^3 \right) \epsilon \right] C_F^2 \nonumber \\
%%% C_F T
  && + \left[ \left( \frac{1}{16} + \frac{1}{4} \Lmu \right)
  \frac{1}{\epsilon} + \frac{947}{288} - \frac{5}{16} \pi^2
  +\frac{11}{24} \Lmu + \frac{3}{8} \Lmu^2 + \left( \frac{17971}{1728} -
  \frac{445}{288} \pi^2 \right.\right. \nonumber \\ 
  && \left.\left. + 2
  \pi^2 \ln 2 - \frac{85}{12} \zeta(3) + \left( \frac{1043}{144} -
  \frac{29}{48} \pi^2 \right) \Lmu + \frac{5}{8} \Lmu^2 + \frac{7}{24}
  \Lmu^3 \right) \epsilon \right] C_FT \nonumber \\
%%% C_F T n_l
  && + \left[ -\frac{1}{8 \epsilon^2} + \frac{11}{48 \epsilon} +
  \frac{113}{96} + \frac{1}{12} \pi^2 + \frac{19}{24} \Lmu + \frac{1}{8}
  \Lmu^2 + \right. \nonumber \\
  && \left. + \left( \frac{851}{192} + \frac{127}{288} \pi^2 +
  \zeta(3) + \left( \frac{145}{48} + \frac{3}{16} \pi^2 \right) \Lmu +
  \frac{23}{24} \Lmu^2 + \frac{1}{8} \Lmu^3 \right) \epsilon
  \right] C_FTn_l \,
  \label{eq::Z22l}
\end{eqnarray}
and
%%%%% 3-loop C_F T n_l
\begin{eqnarray}
%%% C_A
  \delta Z_2^{(3,nl)} &=&  C_FTn_l \left\{ \left[ \frac{11}{72
  \epsilon^3} - \frac{169}{432 \epsilon^2} + \left( \frac{313}{1296} +
  \frac{1}{4} \zeta(3) \right) \frac{1}{\epsilon}
  \right.\right. \nonumber \\
  && + \frac{111791}{15552} + \frac{13}{48} \pi^2 + \frac{47}{36}
  \pi^2 \ln 2 - \frac{2}{9} \pi^2 \ln^2 2 - \frac{35}{72} \zeta(3) +
  \frac{19}{1080} \pi^4  - \frac{1}{9} \ln^4 2 - \frac{8}{3} a_4
  \nonumber \\
  && \left. +
  \left( \frac{169}{27} - \frac{1}{18} \pi^2 + \frac{1}{3} \pi^2 \ln 2 +
  \frac{1}{4} \zeta(3) \right) \Lmu + \frac{469}{288} \Lmu^2 +
  \frac{11}{72} \Lmu^3 \right] C_A \nonumber \\
%%% C_F
  && + \left[ \frac{3}{32 \epsilon^3} + \left( -\frac{19}{192} +
  \frac{3}{32} \Lmu \right) \frac{1}{\epsilon^2} +
  \left( -\frac{235}{384} - \frac{7}{128} \pi^2 - \frac{1}{4} \zeta(3) -
  \frac{41}{64} \Lmu - \frac{3}{64} \Lmu^2 \right) \frac{1}{\epsilon}
  \right. \nonumber \\
  && -\frac{3083}{2304} + \frac{2845}{2304} \pi^2 -
  \frac{47}{18} \pi^2 \ln 2 + \frac{4}{9} \pi^2 \ln^2 2 +
  \frac{473}{96} \zeta(3) - \frac{229}{2160} \pi^4 + \frac{2}{9}
  \ln^4 2 + \frac{16}{3} a_4 \nonumber \\
  &&\left. + \left( -\frac{1475}{384} + \frac{133}{384} \pi^2 -
  \frac{2}{3} \pi^2 \ln 2 + \frac{1}{4} \zeta(3) \right) \Lmu -
  \frac{179}{128} \Lmu^2 - \frac{11}{64} \Lmu^3 \right] C_F \nonumber \\
%%% T
  && + \left[ \left( \frac{1}{36} + \frac{1}{12} \Lmu \right)
  \frac{1}{\epsilon^2} + \left( -\frac{5}{216} + \frac{1}{144} \pi^2 -
  \frac{1}{9} \Lmu + \frac{1}{24} \Lmu^2 \right) \frac{1}{\epsilon}
  \right. \nonumber \\
  && \left. - \frac{4721}{1296} + \frac{19}{54} \pi^2 - \frac{1}{36}
  \zeta(3) + \left( -\frac{329}{108} + \frac{25}{144} \pi^2 \right)
  \Lmu - \frac{7}{12} \Lmu^2 - \frac{5}{72} \Lmu^3 \right] T \nonumber \\
%%% T n_l
  && + \left[ -\frac{1}{36 \epsilon^3} + \frac{11}{216 \epsilon^2} +
  \frac{5}{1296 \epsilon} - \frac{5767}{7776} - \frac{19}{108} \pi^2 -
  \frac{7}{18} \zeta(3) \right. \nonumber \\
  && \left.\left. - \left( \frac{167}{216} + \frac{1}{18} \pi^2 \right)
  \Lmu - \frac{19}{72} \Lmu^2 - \frac{1}{36} \Lmu^3 \right] Tn_l \right\}
  \label{eq::Z23lnl}
  \,,
\end{eqnarray}
where $\Lmu = \ln(\mu^2/m_Q^2)$,
$C_F=(N_c^2-1)/(2N_c)$, $C_A=N_c$, $T=1/2$ and $n_l=n_f-1$ is the number of
light-quark flavours. $N_c=3$ and $\alpha_s$ is the strong coupling renormalized
in the $\overline{\rm MS}$ scheme with $n_f$ active flavours.
$\zeta(3)$ is Riemann's $\zeta$ function with the
value $\zeta(3)=1.202057\ldots$, $\gamma_E$ is Euler's constant with the
value $\gamma_E = 0.577216\ldots$ and $a_4 = {\rm Li}_4(1/2) =
0.517479\ldots$.
Note that $\delta Z_2^{(1)}$ and $\delta Z_2^{(2)}$ are given to order
$\epsilon^2$ and $\epsilon$, respectively, which is needed for the
three-loop calculation of $c_v$.

The integrals needed for $Z_2$ constitute a subset of the ones needed
for $c_v$. In fact, from Eq.~(\ref{eq::integrals}) only the integral types
with a ``$+$'' sign contribute. 
We also want to mention that our calculation poses the first
independent check of the $n_l$ part of Ref.~\cite{Melnikov:2000zc}.

%- }}}
%- {{{ Matching coefficients:

\section{\label{sec::match}Matching coefficient}

It is convenient to cast the perturbative expansion of the
matching coefficient in the form
\begin{eqnarray}
  c_v &=& 1 + \frac{\alpha_s(\mu)}{\pi} c_v^{(1)}
  + \left(\frac{\alpha_s(\mu)}{\pi}\right)^2 c_v^{(2)}
  + \left(\frac{\alpha_s(\mu)}{\pi}\right)^3 c_v^{(3)}
  + {\cal O}(\alpha_s^4)
  \,,
  \label{eq::cvdef}
\end{eqnarray}
where the one-~\cite{KalSar} and
two-loop~\cite{Czarnecki:1997vz,Beneke:1997jm} (see also
Ref.~\cite{Kniehl:2006qw}) terms are given by
\begin{eqnarray}
  c_v^{(1)}&=&-2C_F
  \nonumber\\
  c_v^{(2)}&=&\left(-\frac{151}{72}
  +\frac{89}{144}\pi^2
  -\frac{5}{6}\pi^2\ln2-\frac{13}{4}\zeta(3)\right)C_AC_F
  \nonumber\\&&\mbox{}
  +\left(\frac{23}{8}-\frac{79}{36}\pi^2
  +\pi^2\ln2-\frac{1}{2}\zeta(3)\right)C_F^2
  +\left(\frac{22}{9}-\frac{2}{9}\pi^2\right)C_FT
  \nonumber\\
  &&+\frac{11}{18}C_FTn_l
    -\frac{1}{2}\left[4\beta_0+\pi^2\left(\frac{1}{2}C_A +
  \frac{1}{3}C_F\right)\right]C_F\ln\frac{\mu^2}{m_Q^2}
    \,,
  \label{eq::cv12}
\end{eqnarray}
with $\beta_0=(11 C_A/3-4Tn_f/3)/4$. Note that the terms
proportional to $\beta_0 \ln\frac{\mu^2}{m_Q^2}$ are connected to the 
choice $\alpha_s(\mu)$ in Eq.~(\ref{eq::cvdef}) whereas the ones
proportional to $\pi^2 \ln\frac{\mu^2}{m_Q^2}$ originate from the
separation of hard and soft scales in the construction of NRQCD.

The decomposition of $c_v^{(3)}$ according to the colour structures 
is given by
\begin{eqnarray}
  c_v^{(3)} &=& C_F T n_l\left(
  C_F c_{FFL} + C_A c_{FAL} + T c_{FHL} + T n_l c_{FLL} 
  \right)
  + \mbox{non-$n_l$ terms}
  \label{eq::cv3lnl}
  \,.
\end{eqnarray}

In the notation of Eq.~(\ref{eq::def_cv}) it is understood that
in $\Gamma_v$ the renormalization of the coupling constant and heavy
quark mass is already performed. The renormalization of the
coupling constant can be done straightforwardly by replacing the bare
coupling with the renormalization constant times the renormalized
coupling. The mass renormalization, however, is more intricate since we
are dealing with on-shell integrals. It is therefore convenient to
calculate the counterterms directly by considering one- and two-loop
diagrams with zero-momentum insertions in the massive fermion lines. Sample
diagrams are shown in Fig.~\ref{fig::sample}(d) and (h). The vertex
denoted by a cross is replaced by the mass renormalization
constant in the on-shell scheme. One-loop diagrams with two insertions
are not needed at the order considered in this paper.

Starting from three-loop level the wave function renormalization
constant begins to depend on the QCD gauge parameter,
$\xi$~\cite{Melnikov:2000zc}.
However, it is interesting to note that the colour structures entering
our result (cf. Eq.~(\ref{eq::Z23lnl})) are independent of $\xi$. As a
consequence the quantity $\Gamma_v$ also has to be 
independent of $\xi$ which is indeed the case in our calculation.
Actually, the $\xi$ dependence remaining in the sum of all genuine
three-loop diagrams contributing to $\Gamma_v$ is canceled by the
contribution where a two-loop 
mass counterterm is inserted in the one-loop vertex diagram
(cf. Fig.~\ref{fig::sample}(d)).
Note that the mass counterterm contribution from the two-loop
fermionic diagram (cf. Fig.~\ref{fig::sample}(h)) is $\xi$ independent.
The cancellation of the gauge parameter serves as a welcome check for
the correctness of our result.

As already mentioned above, after taking into account all counterterm
contributions there are still infra-red divergences left in the
quantity $Z_2\Gamma_v$ and thus in $c_v\tilde{Z}_v^{-1}$. 
This is a consequence of the threshold expansion
accompanied by dimensional regularization which is used to extract the 
matching coefficient.
Alternatively it would have been
possible to choose a cut-off for the momentum integrations which
finally results in a factorization scale separating the hard and soft momenta.
In our approach the infra-red poles are absorbed into
$\tilde{Z}_v$ which we define in the $\overline{\rm MS}$ scheme.
As a consequence the matching coefficient itself is finite but scale
dependent. 
In the physical quantities this scale dependence gets canceled 
against the corresponding contributions from the effective theory.
Of course, one has to make sure that for the subtraction of the poles
in the effective theory the same scheme is used as in the full theory.

For the renormalization constant $\tilde{Z}_v$ we obtain
\begin{eqnarray}
  \tilde{Z}_v &=& 1 + \left(\frac{\alpha_s(\mu)}{\pi}\right)^2 \left(
  \frac{1}{12} C_F^2 + \frac{1}{8} C_FC_A \right) \frac{\pi^2}{\epsilon}
  \nonumber \\
  && + \left(\frac{\alpha_s(\mu)}{\pi}\right)^3 C_FTn_l \left[ \left(
  \frac{1}{54} C_F + \frac{1}{36} C_A \right) \frac{\pi^2}{\epsilon^2} -
  \left( \frac{25}{324} C_F + \frac{37}{432} C_A \right)
  \frac{\pi^2}{\epsilon} \right] \nonumber \\
  && + \dots \,, \label{eq::Zv}
\end{eqnarray}
where the ellipses stand for non-$n_l$ and $\order{\alpha_s^4}$ terms.

Our final result for $c_v^{(3)}$ reads
\begin{eqnarray}
  c_{FFL} &=& 46.7(1) + \left( -\frac{17}{12} + \frac{61}{36} \pi^2 -
  \frac{2}{3} \pi^2 \ln 2 + \frac{1}{3} \zeta(3) \right) \Lmu +
  \frac{1}{18} \pi^2 \Lmu^2\,, \nonumber\\
  c_{FAL} &=& 39.6(1) + \left( \frac{181}{54} - \frac{67}{432} \pi^2 +
  \frac{5}{9} \pi^2 \ln 2 + \frac{13}{6} \zeta(3) \right) \Lmu +
  \left( \frac{11}{9} + \frac{1}{12} \pi^2 \right) \Lmu^2\,,
  \nonumber\\
  c_{FHL} &=& -\frac{557}{162} + \frac{26}{81} \pi^2
  + \left( -\frac{55}{27} + \frac{4}{27} \pi^2 \right) \Lmu -
  \frac{4}{9} \Lmu^2 \,, \nonumber\\
  c_{FLL} &=& -\frac{163}{162} - \frac{4}{27} \pi^2
  - \frac{11}{27} \Lmu - \frac{2}{9} \Lmu^2\,.
  \label{eq::cv3}
\end{eqnarray}
The uncertainties assigned to the numerical constants in $c_{FFL}$ and
$c_{FAL}$ are based on a conservative estimate. Note that the precision
of these quantities is more than enough for all phenomenological
applications.
Inserting the numerical values for the
colour factors we obtain for $\mu=m_Q$
\begin{eqnarray}
  c_v^{(3)} &\approx& -0.823\, n_l^2 +121.\, n_l + \mbox{non-$n_l$ terms}
  \,.
  \label{eq::cv3num}
\end{eqnarray} 
The coefficients in Eq.~(\ref{eq::cv3}) correspond to an expansion
parameter $\alpha_s(\mu)$, as given in Eq.~(\ref{eq::cvdef}). Choosing
instead $\alpha_s(m_Q)$ leads to 
\begin{eqnarray}
  \bar{c}_{FFL} &=& 46.7(1) + \frac{25}{108}\pi^2 \Lmu -
  \frac{1}{18} \pi^2 \Lmu^2\,, \nonumber\\
  \bar{c}_{FAL} &=& 39.6(1) + \frac{37}{144}\pi^2 \Lmu - \frac{1}{12}\pi^2
  \Lmu^2\,, \nonumber\\
  \bar{c}_{FHL} &=& -\frac{557}{162} + \frac{26}{81} \pi^2
  \,, \nonumber\\
  \bar{c}_{FLL} &=& -\frac{163}{162} - \frac{4}{27} \pi^2\,.
  \label{eq::cv3mq}
\end{eqnarray} 
The dependence on $L_\mu$ in Eq.~(\ref{eq::cv3mq}) is canceled against
contributions from the effective theory. They agree 
with the ones of Ref.~\cite{Kniehl:2002yv}\footnote{Note that in
  Ref.~\cite{Kniehl:2002yv} there is a typo in the coefficient of the
  $\Lmu^2$ term. The ``$-3/2$'' should read ``$+1$''.}.

%- }}}
%- {{{ Phenomenological applications:

\section{\label{sec::phen}Phenomenological applications}

Let us in this section estimate the numerical effect of our new terms
on the bottom and top system. In particular we consider the 
decay of the $\Upsilon(1S)$ bound state to leptons and the production
of top quark pairs close to threshold.

Next to the matching coefficient considered in the previous Sections 
a crucial ingredient for these quantities 
is the wave function at the origin. Currently the second order is
known
completely~\cite{Kuhn:1998uy,Penin:1998zh,Penin:1998kx,Penin:1998ik,Melnikov:1998ug}
and at order $\alpha_s^3$ the
quadratically~\cite{Kniehl:1999mx,Manohar:2000kr} and 
linearly~\cite{Kniehl:2002yv,Hoang:2003ns} enhanced logarithms
and the corrections proportional to
$\beta_0^3$~\cite{Penin:2005eu,Beneke:2005hg} are available. 
It is convenient to introduce the quantity
$\rho_1=|\psi_1(0)|^2/|\psi_1^C(0)|^2$ where the Coulomb wave function is
given by 
$\left|\psi^C_n(0)\right|^2=C_F^3\alpha_s^3m_q^3/(8\pi n^3)$.
Let us for completeness list the results for principle quantum number
$n=1$
\begin{eqnarray}
  \rho_1 &=& 1+\frac{\alpha_s(\mu_s)}{\pi}
  \left[\left(4-\frac{2}{3}\pi^2\right)\beta_0+\frac{3}{4}a_1\right]
  \nonumber \\
  &&{}+\left(\frac{\alpha_s(\mu_s)}{\pi}\right)^2\left\{
    \left[-C_AC_F+\left(-2+\frac{2}{3}S(S+1)\right)C_F^2\right]
  \pi^2\ln(C_F\alpha_s(\mu_s)) \right. \nonumber\\&&{}
    +\left(-\frac{5}{3}\pi^2+20\zeta(3)+\frac{1}{9}\pi^4\right)\beta_0^2
    +\left(4-\frac{2}{3}\pi^2\right)\beta_1
    +\left(\frac{5}{2}-\frac{2}{3}\pi^2\right)\beta_0a_1
    \nonumber \\&&{}
    +\frac{3}{16}a_1^2+\frac{3}{16}a_2
    +\left.\frac{9}{4}\pi^2C_AC_F
    +\left(\frac{33}{8}-\frac{13}{9}S(S+1)\right)\pi^2C_F^2\right\}
    \nonumber\\&&{}
    +\left(\frac{\alpha_s(\mu_s)}{\pi}\right)^3
    \bigg\{ \pi^2 {\cal C}_2 \ln^2\left(C_F\alpha_s(\mu_s)\right)
    +\pi^2 {\cal C}_1\ln\left(C_F\alpha_s(\mu_s)\right)
    +{\cal C}_0^{\beta_0^3} 
    +\ldots
    \bigg\}
    \,,
    \label{eq::rho1}
\end{eqnarray}
with 
\begin{eqnarray}
  {\cal C}_2 &=&
  \left(-2C_AC_F+
  \left(-4+\frac{4}{3}S(S+1)\right)C_F^2\right)\beta_0
  -\frac{2}{3}C_A^2C_F
  \nonumber\\&&{}\mbox{}
  +\left(-\frac{41}{12}+\frac{7}{12}S(S+1)\right)C_AC_F^2
  -\frac{3}{2}C_F^3
  \,,
\end{eqnarray}
\begin{eqnarray}
  {\cal C}_1&=&\left[\left(-3+\frac{2}{3}\pi^2\right)C_AC_F+
    \left(\frac{4}{3}\pi^2-\left(\frac{10}{9}+\frac{4}{9}\pi^2
    \right) S(S+1)\right) C_F^2\right]\beta_0
  \nonumber\\&&{}
  +\left[-\frac{3}{4}C_AC_F+\left(-\frac{9}{4}
    +\frac{2}{3}S(S+1)\right) C_F^2\right]a_1
  +\frac{1}{4}C_A^3 + \left(\frac{59}{36}-4\ln2\right)C_A^2C_F
  \nonumber \\&&{}
  +\left(\frac{143}{36} - {4}\ln 2-\frac{19}{108} S(S+1)\right)C_AC_F^2
  +\left(-\frac{35}{18}+8\ln2-\frac{1}{3} S(S+1)\right)C_F^3
\nonumber\\&&{}
  +\left(-\frac{32}{15}+2\ln2+\left(1-\ln2\right)S(S+1)\right) C_F^2T
  +\frac{49}{36}C_AC_FTn_l
\nonumber\\&&{}
  +\left(\frac{8}{9}-\frac{10}{27}S(S+1)\right)C_F^2Tn_l\,
  \label{cln}
\end{eqnarray}
and
\begin{eqnarray}
  {\cal C}_0^{\beta_0^3} &=&\beta_0^3
  \Bigg[
    -{20}+\frac{22}{3}\pi^2 + 112\zeta(3) - \frac{7}{5}\pi^4 -
  {12\pi^2\zeta(3)} - 40\zeta(5) - 16\zeta(3)^2
  \nonumber\\&&\mbox{}
  + \frac{4}{105}\pi^6
  \Bigg]
  \,,
\end{eqnarray}
where $\mu_s=C_F m_q\alpha_s(\mu_s)$ is the soft scale. 
It is straightforward to obtain the result for general $\alpha_s(\mu)$
using standard renormalization group analyses.
In Eqs.~(\ref{eq::rho1})--(\ref{cln}) $S$ is the spin quantum number which is equal
to one in our applications. The ellipses in Eq.~(\ref{eq::rho1})
represent yet unknown corrections like, e.g., the pure ultrasoft
contributions. 

For completeness we provide the
one- and two-loop coefficients of the $\beta$ function and the static
potential which are given by
\begin{eqnarray}
  \beta_0 &=& \frac{1}{4} \left( \frac{11}{3} C_A - \frac{4}{3} Tn_l
  \right)
  \,,\nonumber\\
  \beta_1 &=& \frac{1}{16} \left( \frac{34}{3} C_A^2 - \frac{20}{3}
  C_ATn_l - 4 C_FTn_l \right)
  \,,\nonumber\\
  a_1 &=& \frac{31}{9} C_A - \frac{20}{9} Tn_l
  \,,\nonumber\\
  a_2 &=& \left( \frac{4343}{162} + 4 \pi^2 + \frac{22}{3} \zeta(3) -
  \frac{1}{4} \pi^4 \right) C_A^2 - \left( \frac{1798}{81} +
  \frac{56}{3} \zeta(3) \right) C_ATn_l \nonumber\\
  && - \left( \frac{55}{3} - 16 \zeta(3) \right) C_FTn_l + \left(
  \frac{20}{9} Tn_l \right)^2
  \,.
\end{eqnarray}

\subsection{Bottom system}

The leptonic decay of the $\Upsilon(1S)$ state can be cast in the
form~\cite{Kuhn:1998uy,Penin:1998zh,Penin:1998kx,Penin:1998ik,Melnikov:1998ug} 
\begin{eqnarray}
  \Gamma(\Upsilon(1S)\to l^+ l^-) &=& \Gamma^{\rm LO} \rho_1
  \left[c_v^2(m_b)+\frac{C_F^2\alpha_s^2(\mu_s)}{12}c_v(m_b)\left(d_v(m_b)+3\right)
    \right]
  + \ldots
  \,,
  \label{eq::gamll}
\end{eqnarray}
with $\Gamma^{\rm LO}=4\pi
N_cQ_b^2\alpha^2|\psi_1^C(0)|^2/\left(3m_b^2\right)$, $Q_b=-1/3$, $N_c =
3$, $\alpha$ is Sommerfeld's fine-structure constant
and nonperturbative contributions to Eq.~(\ref{eq::gamll}) are ignored.

Inserting the perturbative expansion for $\rho_1$ and $c_v$ we obtain
\begin{eqnarray}
  \Gamma_1 
  &\approx& \Gamma_1^{\rm LO}
  \left(1
  -1.70\,\alpha_s(m_b)
  -7.98\,\alpha_s^2(m_b)
  + 30.0\,\alpha_s^3(m_b)|_{n_l}
  +\ldots\right)
  \nonumber\\
  &&{}\times\Big[ 1
  -0.30\,\alpha_s(\mu_s)
  + \alpha_s^2(\mu_s) \left( 17.2 -5.19 \ln\alpha_s(\mu_s) \right)
  \nonumber\\&&{}\left.
  + \alpha_s^3(\mu_s) \left(
  -14.4 \ln^2\alpha_s(\mu_s)
  +0.17 \ln\alpha_s(\mu_s)
  -34.9 |_{\beta_0^3} \right)
  +\ldots\right]
  \,,
  \label{eq::gamser}
\end{eqnarray}
where in the matching coefficient $\mu=m_b$ has been chosen and 
the corresponding strong coupling is defined with five and
$\alpha_s(\mu_s)$ with four active flavours.

Starting from $\alpha_s(M_Z)=0.118$ we have used the program {\tt
  RunDec}~\cite{Chetyrkin:2000yt} to obtain
$\alpha_s(m_b)=0.2096$ and $\alpha_s(\mu_s)=0.2967$
with $m_b=5.3$~GeV and $\mu_s=2.0967$~GeV which leads to 
\begin{eqnarray}
  \Gamma_1 &\approx& \Gamma_1^{\rm LO}(1 -0.446_{\rm NLO}
  +1.75_{\rm NNLO}-1.20_{\rm N^3LO^\prime}+\ldots)
  \,,
  \label{eq::gamnum}
\end{eqnarray}
where Eq.~(\ref{eq::gamser}) is expanded and terms of order
$\alpha_s^4$ are dropped consistently. The prime reminds that the
third-order corrections are not complete. Apart form the new
contribution to $c_v$ and the known third-order corrections to $\rho_1$
we have also included all interference terms which are proportional
to powers of $n_l$.
Note that the new corrections are responsible for the
reduction of N$^3$LO$^\prime$ terms from $-1.67$ to $-1.20$ which
amounts to about 47\% of the Born cross section.
However, in total the fermionic corrections tend to further reduce the
strong increase of the perturbative coefficients leading to an overall
correction factor of approximately 10\%.

We have performed the numerical analysis also for $m_b=4.8$~GeV. There
are changes in the individual contributions to $\Gamma_1$
(c.f. Eq.~(\ref{eq::gamnum})) of the order of a few 
per cent, however, the final correction factor remains the same.

\subsection{Top system}

In the top-quark case, the nonperturbative effects are negligible.
However, the effect of the top-quark total decay width $\Gamma_t$ has to be
properly taken into account \cite{Fadin:1987wz}, as it is relatively large and
smears out the Coulomb-like resonances below threshold.
The NNLO analysis of the cross
section \cite{Hoang:2000yr} shows that only the ground-state pole gives rise to a
prominent resonance.

A crucial quantity in connection to the threshold production of top
quark pairs is the peak of the normalized cross section  
$R=\sigma(e^+e^-\to t\bar t)/\sigma(e^+e^-\to\mu^+\mu^-)$.
It is dominated by the contribution from the would-be toponium
ground-state, which can be cast in the from
\begin{eqnarray}
  R_1(e^+e^-\to t\bar{t}) &=& R_1^{\rm LO}\rho_1
  \left[c_v^2(m_t)+\frac{C_F^2\alpha_s^2(\mu_s)}{12}c_v(m_t)\left(d_v(m_t)+3\right)
    \right]
  + \ldots
  \,,
  \label{eq::ree}
\end{eqnarray}
with the leading order term
$R_1^{\rm LO}=6\pi N_cQ_t^2|\psi_1^C(0)|^2/\left(m_t^2
\Gamma_t\right)$.
The contributions from the higher Coulomb-like poles and the continuum are not
included in Eq.~(\ref{eq::ree}).

The analog equations to Eqs.~(\ref{eq::gamser}) and~(\ref{eq::gamnum})
read
\begin{eqnarray}
  R_1   &\approx& R_1^{\rm LO}
  \left(1
  -1.70\,\alpha_s(m_t)
  -7.89\,\alpha_s^2(m_t)
  + 37.2\,\alpha_s^3(m_t)|_{\rm n_l}
  +\ldots\right)
  \nonumber\\
  &&{}\times\Big[1
  -0.43\,\alpha_s(\mu_s)
  + \alpha_s^2(\mu_s) \left( 16.1 - 5.19 \ln\alpha_s(\mu_s) \right)
  \nonumber\\&&{}\left.
  + \alpha_s^3(\mu_s) \left(
  -13.8 \ln^2\alpha_s(\mu_s)
  +2.06 \ln\alpha_s(\mu_s)
  -27.2 |_{\beta_0^3} \right)
  +\ldots\right]
  \,,
  \label{eq::rser}
\end{eqnarray}
where $\alpha_s(m_t)$ and $\alpha_s(\mu_s)$ are defined with six and five
active flavours, respectively. Using again $\alpha_s(M_Z)=0.118$ 
one gets $\alpha_s(m_t)=0.1075$
and $\alpha_s(\mu_s)=0.1398$ with
$m_t=175$~GeV and $\mu_s=32.625$~GeV and finally
\begin{eqnarray}
  R_1 &\approx& R_1^{\rm LO}(1 -0.243_{\rm NLO}
  +0.435_{\rm NNLO}-0.195_{\rm N^3LO^\prime}+\ldots)
  \,.
  \label{eq::rnum}
\end{eqnarray}
The fermionic corrections to $c_v$ are responsible for a
reduction of the third-order coefficient from $-0.268$ to
$-0.195$ and thus amount to moderate 7\% of the leading order term.
Similarly as for the bottom quark case also for the top quark the
perturbative series is alternating and the third-order coefficient 
tends to stabilize the expansion.
It is interesting to note that after the inclusion of our new
terms the total corrections amount to less than 1\%.

%- }}}
%- {{{ Conclusions:

\section{\label{sec::concl}Conclusions and outlook}

This paper deals with the question of establishing a relation for the
vector current between full QCD and NRQCD. The corresponding matching
coefficient, $c_v$, constitutes a building block in all threshold phenomena
involving a vector coupling.
The main result of this paper is the fermionic 
three-loop contribution to $c_v$ which contains a closed light fermion
loop. The numerical effect of the new terms is relatively big and 
amounts to about 47\%
for the bottom and to about 7\% for the top quark system.
However, their inclusion leads to a reduction of the
overall corrections. E.g., in the case of the top quark the corrections
to the normalization at the peak of the production cross section
amount to less than 1\% after including all currently known perturbative terms. 

Several steps still have to be taken in order to arrive at the
complete analysis of the top-anti-top quark system at threshold. Next
to the non-fermionic pieces to $c_v$ the major building blocks which
are still missing are the ultra-soft corrections and the three-loop
static potential.

%- }}}
%- {{{ Acknowledgements:

\vspace*{1em}

\noindent
{\large\bf Acknowledgements}\\
We would like to thank A.A.~Penin and V.A.~Smirnov for useful comments
and discussions.
J.H.P. would like to thank S.~Bekavac for discussions about
Mellin-Barnes integrals and  M.~Kalmykov for useful advice on two-loop
sunset integrals. We thank M.~Tentyukov for providing the interface to
{\tt Fermat}.
This work was supported by the ``Impuls- und Vernetzungsfonds'' of the
Helmholtz Association, contract number VH-NG-008 and 
by the DFG through SFB/TR~9. The Feynman diagrams were
drawn with {\tt JaxoDraw} \cite{Binosi:2003yf}

%- }}}
%- {{{ Appendix:

\begin{appendix}

\section{Results for the master integrals\label{masters}}
\begin{figure}
  \leavevmode
  \epsfxsize=\textwidth
 \epsffile{masters.eps}
\caption{\label{fig::masters} Three-loop master integrals. Bold lines
  denote massive lines with mass $2m_Q$, thin lines denote massive lines
  with mass $m_Q$ and dashed lines denote massless lines. All external
  lines are on-shell. A dot on a line denotes a squared propagator.}
\end{figure}

The master integrals for the three-loop functions of
Eq.~(\ref{eq::integrals}) are shown in Fig.~\ref{fig::masters}. For the
self-energy integrals $I_1$--$I_6$ the external momentum squared equals
the square of the heavy quark mass, $q^2 = m_Q^2$. For the integrals
$I_7$ and $I_8$ we have $q^2 = 4 m_Q^2$.
The momenta $q_1$ and $q_2$ flowing out of the vertex diagrams on the
right fulfill $q_1^2 = q_2^2 = m_Q^2$ with $(q_1 + q_2)^2 = 4 m_Q^2$.
In Minkowski space the results for the master integrals read
\begin{eqnarray}
  I_1 &=& J_+^{(3)}(0,1,0,1,1,0,1,0,1) \nonumber \\
  &=& m_Q^2 \left( \frac{\mu^2}{m_Q^2}  e^{-\gamma_E}\right)^{3 \epsilon}
  \left\{ \frac{1}{3\epsilon^3} +
  \frac{5}{3\epsilon^2} + \left(4 + \frac{3}{4}\pi^2\right)
  \frac{1}{\epsilon} - \frac{10}{3} + \frac{11}{4} \pi^2 + \frac{25}{3}
  \zeta(3) \right. \nonumber \\
  && + \left( -\frac{302}{3} + 2\pi^2 + \frac{89}{3}\zeta(3) +
  \frac{3059}{1440} \pi^4 \right)\epsilon + \left( -734 -
  \frac{69}{2}\pi^2 + 16\zeta(3) + \frac{9751}{1440} \pi^4
  \right. \nonumber \\
  && \left.\left. + \frac{107}{4}\pi^2\zeta(3) +
  \frac{2309}{5}\zeta(5) \right) \epsilon^2 + \order{\epsilon^3}
  \right\} \,,
\\
  I_2 &=& J_+^{(3)}(0,1,1,1,0,1,1,0,0) \nonumber \\
  &=& -m_Q^2 \left( \frac{\mu^2}{m_Q^2} \right)^{3 \epsilon}
  \frac{\Gamma^4(1-\epsilon) \Gamma^2(\epsilon)}{\Gamma^2(2-2\epsilon)}
  \frac{\Gamma(3\epsilon-1) \Gamma(3-6\epsilon)}{\Gamma(3-4\epsilon)} \,,
\end{eqnarray}
\begin{eqnarray}
  I_3 &=& L_+^{(3,n_l)}(0,0,1,1,1,1,0,1,0) \nonumber \\
  &=& m_Q^2 \left( \frac{\mu^2}{m_Q^2} e^{-\gamma_E}\right)^{3 \epsilon}
  \left\{ \frac{2}{3\epsilon^3} +
  \frac{10}{3\epsilon^2} + \left( \frac{26}{3} + \frac{1}{2}\pi^2 \right)
  \frac{1}{\epsilon} + 2 + \frac{9}{2}\pi^2 + \frac{14}{3}\zeta(3)
  \right. \nonumber \\
  && + \left( -\frac{398}{3} + \frac{53}{2}\pi^2 - 16\pi^2 \ln 2 +
  \frac{238}{3}\zeta(3) + \frac{287}{720}\pi^4 \right) \epsilon
  + \left( -1038 + \frac{259}{2}\pi^2 \right. \nonumber \\
  && - 160\pi^2 \ln 2 + \frac{128}{3} \pi^2 \ln^2 2 +
  \frac{1862}{3}\zeta(3) + \frac{323}{720}\pi^4 +
  \frac{7}{2}\pi^2\zeta(3) + \frac{478}{5}\zeta(5) \nonumber \\
  && \left. + \frac{64}{3}\ln^4 2 + 512 a_4 \right) \epsilon^2 +
  \order{\epsilon^3} \bigg\} \,,
\\
  I_4 &=& J_+^{(3)}(0,0,0,1,1,0,1,0,1) \nonumber \\
  &=& m_Q^4 \left( \frac{\mu^2}{m_Q^2} \right)^{3 \epsilon}
  \frac{\Gamma^2(1-\epsilon) \Gamma(\epsilon) \Gamma^2(2\epsilon-1)
  \Gamma(3\epsilon-2)}{\Gamma(4\epsilon-2) \Gamma(2-\epsilon)} \,,
\\
  I_5 &=& J_+^{(3)}(0,1,0,1,0,0,1,0,1) \nonumber \\
  &=& m_Q^4 \left( \frac{\mu^2}{m_Q^2} \right)^{3 \epsilon}
  \Gamma(\epsilon-1) \frac{\Gamma^2(1-\epsilon)
  \Gamma(\epsilon)}{\Gamma(2-2\epsilon)} \frac{\Gamma(2\epsilon-1)
  \Gamma(3-4\epsilon)}{\Gamma(3-3\epsilon)} \,,
\\
  I_6 &=& J_+^{(3)}(0,0,1,1,1,0,1,0,0) \nonumber \\
  &=& m_Q^4 \left( \frac{\mu^2}{m_Q^2} \right)^{3 \epsilon}
  \frac{\Gamma^3(1-\epsilon) \Gamma(2\epsilon-1) \Gamma(3\epsilon-2)
  \Gamma(5-6\epsilon)}{\Gamma(3-3\epsilon) \Gamma(4-4\epsilon)} \,,
\\
  I_7 &=& J_-^{(3)}(0,0,0,1,1,0,1,0,1) \nonumber \\
  &=& m_Q^4 \left( \frac{\mu^2}{m_Q^2} e^{-\gamma_E} \right)^{3 \epsilon}
  \left\{ \frac{1}{3\epsilon^3} +
  \frac{1}{2\epsilon^2} + \left(-\frac{17}{36} + \frac{1}{12} \pi^2
  \right)\frac{1}{\epsilon} - 6.6827387(1) \right. \nonumber \\
  && - 56.300353(1)\epsilon - 209.48231(1)\epsilon^2 + \order{\epsilon^3}
  \bigg\} \,,
\\
  I_8 &=& J_-^{(3)}(0,0,0,2,1,0,1,0,1) \nonumber\\
  &=& -m_Q^2 \left( \frac{\mu^2}{m_Q^2} e^{-\gamma_E} \right)^{3 \epsilon}
  \left\{ \frac{1}{3\epsilon^3} +
  \frac{2}{3\epsilon^2} + \left( -\frac{8}{3} + \frac{5}{6}\pi^2 \right)
  \frac{1}{\epsilon} + 10.797602(1)
  \right. \nonumber\\
  && + 62.250613(1)\epsilon + \order{\epsilon^2} \bigg\} \,,
\\
  I_9 &=& J_-^{(3)}(0,1,0,1,1,0,1,0,1) \nonumber \\
  &=& m_Q^2 \left( \frac{\mu^2}{m_Q^2} e^{-\gamma_E} \right)^{3 \epsilon}
  \left\{ \frac{1}{3\epsilon^3} +
  \frac{5}{3\epsilon^2} + \left( \frac{10}{3} + \frac{3}{4}\pi^2 \right)
  \frac{1}{\epsilon} + 33.8328(4) + 152.870(4)\epsilon \right. \nonumber \\
  &&+ \order{\epsilon^2} \bigg\} \,,
\end{eqnarray}
\begin{eqnarray}
  I_{10} &=& L_+^{(3,n_l)}(0,0,1,1,1,0,1,1,0) \nonumber\\
  &=& m_Q^2 \left( \frac{\mu^2}{m_Q^2} e^{-\gamma_E} \right)^{3 \epsilon}
  \left\{ \frac{2}{3\epsilon^3} +
  \frac{10}{3\epsilon^2} + \left(8 + \frac{1}{2} \pi^2 \right)
  \frac{1}{\epsilon} + 52.5698(4) + 145.087(4)\epsilon \right. \nonumber \\
  && + 562.250(14)\epsilon^2 + \order{\epsilon^3}
  \bigg\}  \,,
\\
  I_{11} &=& J_-^{(3)}(0,1,0,1,0,1,1,0,1) \nonumber \\
  &=& m_Q^2 \left( \frac{\mu^2}{m_Q^2} e^{-\gamma_E} \right)^{3 \epsilon}
  \left\{ \frac{1}{2\epsilon^3} +
  \frac{11}{6\epsilon^2} + \left( -\frac{1}{6} + \frac{29}{24}\pi^2
  \right) \frac{1}{\epsilon} + 34.791(4) + 154.08(2)\epsilon
  \right. \nonumber \\
  &&+ \order{\epsilon^2} \bigg\} \,,
\\
  I_{12} &=& J_-^{(3)}(0,1,0,1,0,1,1,0,2) \nonumber \\
  &=& \left( \frac{\mu^2}{m_Q^2} e^{-\gamma_E} \right)^{3\epsilon}
  \left\{ \frac{1}{6\epsilon^3} + \frac{1}{2\epsilon^2} +
  \left( \frac{1}{6} + \frac{7}{24}\pi^2 \right) \frac{1}{\epsilon} +
  7.024(4) + 32.04(2)\epsilon \right. \nonumber \\
  && + \order{\epsilon^2} \bigg\} \,,
\end{eqnarray}
where $\zeta(5) = 1.036927\ldots$ and $a_4 = {\rm Li}_4(1/2) =
0.517479\ldots$.

$I_1$ was calculated in Ref.~\cite{Melnikov:2000zc}. We have repeated
the calculation described in the reference with the program
\texttt{XSummer}~\cite{Moch:2005uc} and find complete agreement.  The
results for $I_3$--$I_6$ can be found in Ref.~\cite{Laporta:1996mq}. We
have checked the result for $I_3$ numerically with the help of a one-fold
Mellin-Barnes~\cite{Smirnov:1999gc,Tausk:1999vh} representation and find
complete agreement. Note, that $I_2, I_4, I_5$ and $I_6$ are quite
simple to evaluate and are available for general $\epsilon$.
Nevertheless the result for $I_6$ as given in
Ref.~\cite{Melnikov:2000zc} is not correct whereas the correct result
can be found in Ref.~\cite{Laporta:1996mq}. The pole part of $I_7$
agrees with the result of Ref.~\cite{Groote:2004qq}.

The remaining master integrals, $I_7$--$I_{12}$, have been evaluated
with the help of the Mellin-Barnes method where the 
evaluation of the integrals has been performed with the program
\texttt{MB}~\cite{Czakon:2005rk}.
$I_7$ and $I_8$ can be expressed in terms of a one-fold Mellin-Barnes
representation and are thus known with a quite high precision. On the
other hand, $I_{10}$ is represented by a two-fold and 
$I_9$, $I_{11}$ and $I_{12}$ even by a three-fold integration which
results in less accurate results. The quoted uncertainties in the above
equations correspond to twice the Vegas error given by {\tt MB} for the
multi-dimensional integrals and to a conservative estimate in case of
the one-dimensional integrals $I_7$ and $I_8$. Note that the latter
errors are negligible.

\begin{figure}
 \centerline{\epsfig{file=epfin,scale=0.65}}
\caption{\label{fig::epfin} $\epsilon$-finite master integrals. The same
coding as in Fig.~\ref{fig::masters} is adopted.}
\end{figure}

As can be seen in the above results some of the master integrals are
needed to higher order in $\epsilon$ which on one hand makes the
calculation very tedious and on the other hand leads to less accurate
results. For this reason we decided to change basis and switch --- at
least for some of the integrals --- to the so-called
$\epsilon$-finite master integrals~\cite{Chetyrkin:2006dh} which have
the advantage that the 
coefficient function is finite and thus the integral itself is only
needed to order $\epsilon^0$.

Since $I_1$--$I_6$ are known analytically, a replacement is, of
course, not necessary. Furthermore, for $I_7$ and $I_8$ the numerical
precision is 
sufficient for our calculation. 
As far as the remaining four integrals are concerned, we
found it convenient to replace $I_9$ and $I_{12}$ by the integrals
shown in Fig.~\ref{fig::epfin}. 
Their numerical evaluation with the help of \texttt{MB} is
straightforward leading to the results
\begin{eqnarray}
  I_9^f &=& J_-^{(3)}(0,1,1,1,1,0,1,0,1) \nonumber \\
  &=& \left( \frac{\mu^2}{m_Q^2}  e^{-\gamma_E}\right)^{3\epsilon}
  \left\{ \frac{1}{6\epsilon^3} +
  \frac{3}{2\epsilon^2} + \left( \frac{55}{6} + \frac{3}{8} \pi^2
  \right) \frac{1}{\epsilon} + 64.678(8) + \order{\epsilon} \right\} \,,
\\
  I_{12}^f &=& J_-^{(3)}(0,1,0,1,1,1,1,0,1) =  \left(
  \frac{\mu^2}{m_Q^2}  e^{-\gamma_E}\right)^{3\epsilon}
  \left\{ \frac{2.4041(4)}{\epsilon} + 8.1(2) + \order{\epsilon} \right\} \,.
\end{eqnarray}
Switching from $I_9$ and $I_{12}$ to $I_9^f$ and $I_{12}^f$ reduces
the number of coefficients which are only known numerically from
17 to 14. In particular the $1/\epsilon$-poles of $I_9$ and
$I_{12}$ have been determined in analytical form from the
$\epsilon$-finite integrals via the corresponding IBP relations.

A further reduction is achieved after exploiting the locality and
incorporating the knowledge of the linear $\ln(\mu^2/m_Q^2)$ term in
$c_v$~\cite{Kniehl:2002yv} which we have checked numerically in a first
step. In this way the analytical result for the $1/\epsilon$-pole
of $I_{11}$ has been determined. In total ten (eleven) numerical
coefficients contribute to our final result given in
Eqs.~(\ref{eq::cv3}) and~(\ref{eq::cv3num}) using the $\epsilon$-finite
(``normal'') basis.

In the end it turned out, that the numerical precision is better using
the ``normal'' basis. The reason is, that the uncertainty of the
$\epsilon^0$-coefficient of $I_{12}^f$ is by far the largest. Still, it
was useful to consider the $\epsilon$-finite basis since it enabled us
to determine analytical results for two integral coefficients of the
``normal'' basis.

\end{appendix}

%- }}}
%- {{{ bibliography:

\begin{thebibliography}{99}

%\cite{Caswell:1985ui}
\bibitem{Caswell:1985ui}
  W.~E.~Caswell and G.~P.~Lepage,
  %``Effective Lagrangians For Bound State Problems In QED, QCD, And Other Field
  %Theories,''
  Phys.\ Lett.\ B {\bf 167} (1986) 437.
  %%CITATION = PHLTA,B167,437;%%

%\cite{Bodwin:1994jh}
\bibitem{Bodwin:1994jh}
  G.~T.~Bodwin, E.~Braaten and G.~P.~Lepage,
  %``Rigorous QCD analysis of inclusive annihilation and production of heavy
  %quarkonium,''
  Phys.\ Rev.\ D {\bf 51} (1995) 1125
  [Erratum-ibid.\ D {\bf 55} (1997) 5853]
  [arXiv:hep-ph/9407339].
  %%CITATION = HEP-PH 9407339;%%

%\cite{Brambilla:2004wf}
\bibitem{Brambilla:2004wf}
  N.~Brambilla {\it et al.},
  %``Heavy quarkonium physics,''
  arXiv:hep-ph/0412158.
  %%CITATION = HEP-PH 0412158;%%

%\cite{Martinez:2002st}
\bibitem{Martinez:2002st}
  M.~Martinez and R.~Miquel,
  %``Multi-parameter fits to the t anti-t threshold observables at a future  e+ e-
  %linear collider,''
  Eur.\ Phys.\ J.\ C {\bf 27} (2003) 49
  [arXiv:hep-ph/0207315].
  %%CITATION = HEP-PH 0207315;%%

%\cite{Hoang:2000yr}
\bibitem{Hoang:2000yr}
  A.~H.~Hoang {\it et al.},
  %``Top-antitop pair production close to threshold: Synopsis of
  %recent NNLO results,''
  Eur.\ Phys.\ J.\ directC {\bf 2} (2000) 1
  [arXiv:hep-ph/0001286].
  %%CITATION = HEP-PH 0001286;%%

%\cite{Kniehl:2002br}
\bibitem{Kniehl:2002br}
  B.~A.~Kniehl, A.~A.~Penin, V.~A.~Smirnov and M.~Steinhauser,
  %``Potential NRQCD and heavy-quarkonium spectrum at
  %next-to-next-to-next-to-leading order,''
  Nucl.\ Phys.\ B {\bf 635} (2002) 357
  [arXiv:hep-ph/0203166].
  %%CITATION = HEP-PH 0203166;%%

%\cite{Penin:2002zv}
\bibitem{Penin:2002zv}
  A.~A.~Penin and M.~Steinhauser,
  %``Heavy quarkonium spectrum at O(alpha(s)**5 m(q)) and bottom / top  quark
  %mass determination,''
  Phys.\ Lett.\ B {\bf 538} (2002) 335
  [arXiv:hep-ph/0204290].
  %%CITATION = HEP-PH 0204290;%%

%\cite{Penin:2005eu}
\bibitem{Penin:2005eu}
  A.~A.~Penin, V.~A.~Smirnov and M.~Steinhauser,
  %``Heavy quarkonium spectrum and production / annihilation rates to  order
  %beta(0)**3 alpha(s)**3,''
  Nucl.\ Phys.\ B {\bf 716} (2005) 303
  [arXiv:hep-ph/0501042].
  %%CITATION = HEP-PH 0501042;%%

%\cite{Beneke:2005hg}
\bibitem{Beneke:2005hg}
  M.~Beneke, Y.~Kiyo and K.~Schuller,
  %``Third-order Coulomb corrections to the S-wave Green function,
  %energy levels and wave functions at the origin,''
  Nucl.\ Phys.\ B {\bf 714} (2005) 67
  [arXiv:hep-ph/0501289].
  %%CITATION = HEP-PH 0501289;%%

%\cite{Pineda:2001ra}
\bibitem{Pineda:2001ra}
  A.~Pineda,
  %``Renormalization group improvement of the NRQCD Lagrangian and heavy
  %quarkonium spectrum,''
  Phys.\ Rev.\ D {\bf 65} (2002) 074007
  [arXiv:hep-ph/0109117].
  %%CITATION = HEP-PH 0109117;%%

%\cite{Pineda:2001et}
\bibitem{Pineda:2001et}
  A.~Pineda,
  %``Next-to-leading-log renormalization-group running in heavy-quarkonium
  %creation and annihilation,''
  Phys.\ Rev.\ D {\bf 66} (2002) 054022
  [arXiv:hep-ph/0110216].
  %%CITATION = HEP-PH 0110216;%%

%\cite{Hoang:2003ns}
\bibitem{Hoang:2003ns}
  A.~H.~Hoang,
  %``Three-loop anomalous dimension of the heavy quark pair production current
  %in non-relativistic QCD,''
  Phys.\ Rev.\ D {\bf 69} (2004) 034009
  [arXiv:hep-ph/0307376].
  %%CITATION = HEP-PH 0307376;%%

%\cite{Penin:2004ay}
\bibitem{Penin:2004ay}
  A.~A.~Penin, A.~Pineda, V.~A.~Smirnov and M.~Steinhauser,
  %``Spin dependence of heavy quarkonium production and annihilation
  %rates: Complete next-to-next-to-leading logarithmic result,''
  Nucl.\ Phys.\ B {\bf 699} (2004) 183
  [arXiv:hep-ph/0406175].
  %%CITATION = HEP-PH 0406175;%%

%\cite{Hoang:2004tg}
\bibitem{Hoang:2004tg}
  A.~H.~Hoang and C.~J.~Reisser,
  %``Electroweak absorptive parts in NRQCD matching conditions,''
  Phys.\ Rev.\ D {\bf 71} (2005) 074022
  [arXiv:hep-ph/0412258].
  %%CITATION = HEP-PH 0412258;%%

%\cite{Eiras:2006xm}
\bibitem{Eiras:2006xm}
  D.~Eiras and M.~Steinhauser,
  %``Complete Higgs mass dependence of top quark pair threshold production to
  %order alpha alpha(s),''
  arXiv:hep-ph/0605227.
  %%CITATION = HEP-PH 0605227;%%

%\cite{Beneke:1997zp}
\bibitem{Beneke:1997zp}
  M.~Beneke and V.~A.~Smirnov,
  %``Asymptotic expansion of Feynman integrals near threshold,''
  Nucl.\ Phys.\ B {\bf 522} (1998) 321
  [arXiv:hep-ph/9711391].
  %%CITATION = HEP-PH 9711391;%%

%\cite{Smirnov:2002pj}
\bibitem{Smirnov:2002pj}
  V.~A.~Smirnov,
  ``Applied asymptotic expansions in momenta and masses,''
  Springer (2002).
  %\href{http://www.slac.stanford.edu/spires/find/hep/www?irn=4841620}{SPIRES entry}

%\cite{Nogueira:1991ex}
\bibitem{Nogueira:1991ex}
  P.~Nogueira,
  %``Automatic Feynman graph generation,''
  J.\ Comput.\ Phys.\  {\bf 105} (1993) 279.
  %%CITATION = JCTPA,105,279;%%

%\cite{Harlander:1997zb}
\bibitem{Harlander:1997zb}
  R.~Harlander, T.~Seidensticker and M.~Steinhauser,
  %``Complete corrections of O(alpha alpha(s)) to the decay of the Z
  %boson  into
  %bottom quarks,''
  Phys.\ Lett.\ B {\bf 426} (1998) 125
  [hep-ph/9712228].
  %%CITATION = HEP-PH 9712228;%%

%\cite{Seidensticker:1999bb}
\bibitem{Seidensticker:1999bb}
  T.~Seidensticker,
  %``Automatic application of successive asymptotic expansions of
  %Feynman
  %diagrams,''
  hep-ph/9905298.
  %%CITATION = HEP-PH 9905298;%%

%\cite{Laporta:1996mq}
\bibitem{Laporta:1996mq}
  S.~Laporta and E.~Remiddi,
  %``The analytical value of the electron (g-2) at order alpha~3 in QED,''
  Phys.\ Lett.\ B {\bf 379} (1996) 283
  [arXiv:hep-ph/9602417].
  %%CITATION = HEP-PH 9602417;%%

%\cite{Laporta:2001dd}
\bibitem{Laporta:2001dd}
  S.~Laporta,
  %``High-precision calculation of multi-loop Feynman integrals by  difference
  %equations,''
  Int.\ J.\ Mod.\ Phys.\ A {\bf 15} (2000) 5087
  [arXiv:hep-ph/0102033].
  %%CITATION = HEP-PH 0102033;%%

%\cite{Baikov:1996rk}
\bibitem{Baikov:1996rk}
  P.~A.~Baikov,
  %``Explicit solutions of the 3--loop vacuum integral recurrence relations,''
  Phys.\ Lett.\ B {\bf 385} (1996) 404
  [arXiv:hep-ph/9603267].
  %%CITATION = HEP-PH 9603267;%%

%\cite{Baikov:1996iu}
\bibitem{Baikov:1996iu}
  P.~A.~Baikov,
  %``Explicit solutions of the multi-loop integral recurrence relations and  its
  %application,''
  Nucl.\ Instrum.\ Meth.\ A {\bf 389} (1997) 347
  [arXiv:hep-ph/9611449].
  %%CITATION = HEP-PH 9611449;%%

%\cite{Smirnov:2003kc}
\bibitem{Smirnov:2003kc}
  V.~A.~Smirnov and M.~Steinhauser,
  %``Solving recurrence relations for multi-loop Feynman integrals,''
  Nucl.\ Phys.\ B {\bf 672} (2003) 199
  [arXiv:hep-ph/0307088].
  %%CITATION = HEP-PH 0307088;%%

%\cite{Chetyrkin:2006dh}
\bibitem{Chetyrkin:2006dh}
  K.~G.~Chetyrkin, M.~Faisst, C.~Sturm and M.~Tentyukov,
  %``e-finite basis of master integrals for the integration-by-parts method,''
  Nucl.\ Phys.\ B {\bf 742} (2006) 208
  [arXiv:hep-ph/0601165].
  %%CITATION = HEP-PH 0601165;%%

%\cite{Chetyrkin:1981qh}
\bibitem{Chetyrkin:1981qh}
  K.~G.~Chetyrkin and F.~V.~Tkachov,
  %``Integration By Parts: The Algorithm To Calculate Beta Functions In 4 Loops,''
  Nucl.\ Phys.\ B {\bf 192} (1981) 159.
  %%CITATION = NUPHA,B192,159;%%

%\cite{Anastasiou:2004vj}
\bibitem{Anastasiou:2004vj}
  C.~Anastasiou and A.~Lazopoulos,
  %``Automatic integral reduction for higher order perturbative  calculations,''
  JHEP {\bf 0407} (2004) 046
  [arXiv:hep-ph/0404258].
  %%CITATION = HEP-PH 0404258;%%

\bibitem{PMDS}
  P.~Marquard and D.~Seidel,
  unpublished.

%\cite{Bauer:2000cp}
\bibitem{Bauer:2000cp}
  C.~Bauer, A.~Frink and R.~Kreckel,
  %``Introduction to the GiNaC Framework for Symbolic Computation within the C++
  %Programming Language,''
  arXiv:cs.sc/0004015.
  %%CITATION = CS-SC 0004015;%%

\bibitem{fermat} R.~H.~Lewis, Fermat's User Guide,
  http://www.bway.net/\~{}lewis.

%\cite{Tentyukov:2006ys}
\bibitem{Tentyukov:2006ys}
  M.~Tentyukov and J.~A.~M.~Vermaseren,
  %``Extension of the functionality of the symbolic program FORM by external
  %software,''
  arXiv:cs.sc/0604052.
  %%CITATION = CS-SC 0604052;%%

%\cite{Vermaseren:2000nd}
\bibitem{Vermaseren:2000nd}
  J.~A.~M.~Vermaseren,
  %``New features of FORM,''
  arXiv:math-ph/0010025.
  %%CITATION = MATH-PH 0010025;%%

%\cite{Broadhurst:1991fy}
\bibitem{Broadhurst:1991fy}
  D.~J.~Broadhurst, N.~Gray and K.~Schilcher,
  %``Gauge invariant on-shell Z(2) in QED, QCD and the effective field theory of
  %a static quark,''
  Z.\ Phys.\ C {\bf 52} (1991) 111.
  %%CITATION = ZEPYA,C52,111;%%

%\cite{Melnikov:2000zc}
\bibitem{Melnikov:2000zc}
  K.~Melnikov and T.~van Ritbergen,
  %``The three-loop on-shell renormalization of QCD and QED,''
  Nucl.\ Phys.\ B {\bf 591} (2000) 515
  [arXiv:hep-ph/0005131].
  %%CITATION = HEP-PH 0005131;%%

\bibitem{KalSar}
  G. K{\"a}llen and A. Sarby,
  K.\ Dan.\ Vidensk.\ Selsk.\ Mat.-Fis.\ Medd.\ {29},  N17 (1955) 1.

%\cite{Czarnecki:1997vz}
\bibitem{Czarnecki:1997vz}
  A.~Czarnecki and K.~Melnikov,
  %``Two-loop QCD corrections to the heavy quark pair production cross  section
  %in e+ e- annihilation near the threshold,''
  Phys.\ Rev.\ Lett.\  {\bf 80} (1998) 2531
  [arXiv:hep-ph/9712222].
  %%CITATION = HEP-PH 9712222;%%

%\cite{Beneke:1997jm}
\bibitem{Beneke:1997jm}
  M.~Beneke, A.~Signer and V.~A.~Smirnov,
  %``Two-loop correction to the leptonic decay of quarkonium,''
  Phys.\ Rev.\ Lett.\  {\bf 80} (1998) 2535
  [arXiv:hep-ph/9712302].
  %%CITATION = HEP-PH 9712302;%%

%\cite{Kniehl:2006qw}
\bibitem{Kniehl:2006qw}
  B.~A.~Kniehl, A.~Onishchenko, J.~H.~Piclum and M.~Steinhauser,
  %``Two-loop matching coefficients for heavy quark currents,''
  Phys.\ Lett.\ B {\bf 638} (2006) 209
  [arXiv:hep-ph/0604072].
  %%CITATION = HEP-PH 0604072;%%

%\cite{Kniehl:2002yv}
\bibitem{Kniehl:2002yv}
  B.~A.~Kniehl, A.~A.~Penin, M.~Steinhauser and V.~A.~Smirnov,
  %``Heavy-quarkonium creation and annihilation with O(alpha(s)**3 ln(alpha(s)))
  %accuracy. ((V)),''
  Phys.\ Rev.\ Lett.\  {\bf 90} (2003) 212001
  [Erratum-ibid.\ {\bf 91} (2003) 139903(E)]
  [arXiv:hep-ph/0210161].
  %%CITATION = HEP-PH 0210161;%%

%\cite{Kuhn:1998uy}
\bibitem{Kuhn:1998uy}
  J.~H.~K{\"u}hn, A.~A.~Penin and A.~A.~Pivovarov,
  %``Coulomb resummation for b anti-b system near threshold and precision
  %determination of alpha(s) and m(b),''
  Nucl.\ Phys.\ B {\bf 534} (1998) 356
  [arXiv:hep-ph/9801356].
  %%CITATION = HEP-PH 9801356;%%

%\cite{Penin:1998zh}
\bibitem{Penin:1998zh}
  A.~A.~Penin and A.~A.~Pivovarov,
  %``Next-to-next-to-leading order vacuum polarization function of heavy  quark
  %near threshold and sum rules for b anti-b system,''
  Phys.\ Lett.\ B {\bf 435} (1998) 413
  [arXiv:hep-ph/9803363].
  %%CITATION = HEP-PH 9803363;%%

%\cite{Penin:1998kx}
\bibitem{Penin:1998kx}
  A.~A.~Penin and A.~A.~Pivovarov,
  %``Bottom quark pole mass and $|$V(cb)$|$ matrix element from R(e+ e- $\to$ b
  %anti-b) and Gamma(sl)(b $\to$ c l nu/l) in the next-to-next-to-leading
  %order,''
  Nucl.\ Phys.\ B {\bf 549} (1999) 217
  [arXiv:hep-ph/9807421].
  %%CITATION = HEP-PH 9807421;%%

%\cite{Penin:1998ik}
\bibitem{Penin:1998ik}
  A.~A.~Penin and A.~A.~Pivovarov,
  %``Top quark threshold production in gamma gamma collision in the
  %next-to-leading order,''
  Nucl.\ Phys.\ B {\bf 550} (1999) 375
  [arXiv:hep-ph/9810496].
  %%CITATION = HEP-PH 9810496;%%

%\cite{Melnikov:1998ug}
\bibitem{Melnikov:1998ug}
  K.~Melnikov and A.~Yelkhovsky,
  %``The b quark low-scale running mass from Upsilon sum rules,''
  Phys.\ Rev.\ D {\bf 59} (1999) 114009
  [arXiv:hep-ph/9805270].
  %%CITATION = HEP-PH 9805270;%%

%\cite{Kniehl:1999mx}
\bibitem{Kniehl:1999mx}
  B.~A.~Kniehl and A.~A.~Penin,
  %``Order alpha(s)**3 ln**2(1/alpha(s)) corrections to  heavy-quarkonium
  %creation and annihilation,''
  Nucl.\ Phys.\ B {\bf 577} (2000) 197
  [arXiv:hep-ph/9911414].
  %%CITATION = HEP-PH 9911414;%%

%\cite{Manohar:2000kr}
\bibitem{Manohar:2000kr}
  A.~V.~Manohar and I.~W.~Stewart,
  %``Running of the heavy quark production current and 1/v potential in QCD,''
  Phys.\ Rev.\ D {\bf 63} (2001) 054004
  [arXiv:hep-ph/0003107].
  %%CITATION = HEP-PH 0003107;%%

%\cite{Chetyrkin:2000yt}
\bibitem{Chetyrkin:2000yt}
  K.~G.~Chetyrkin, J.~H.~K{\"u}hn and M.~Steinhauser,
  %``RunDec: A Mathematica package for running and decoupling of the strong
  %coupling and quark masses,''
  Comput.\ Phys.\ Commun.\  {\bf 133} (2000) 43
  [arXiv:hep-ph/0004189].
  %%CITATION = HEP-PH 0004189;%%

%\cite{Fadin:1987wz}
\bibitem{Fadin:1987wz}
  V.~S.~Fadin and V.~A.~Khoze,
  %``Threshold Behavior Of Heavy Top Production In E+ E- Collisions,''
  JETP Lett.\  {\bf 46} (1987) 525
  [Pisma Zh.\ Eksp.\ Teor.\ Fiz.\  {\bf 46} (1987) 417].
  %%CITATION = JTPLA,46,525;%%

%\cite{Binosi:2003yf}
\bibitem{Binosi:2003yf}
  D.~Binosi and L.~Theussl,
  %``JaxoDraw: A graphical user interface for drawing Feynman diagrams,''
  Comput.\ Phys.\ Commun.\  {\bf 161} (2004) 76
  [arXiv:hep-ph/0309015].
  %%CITATION = HEP-PH 0309015;%%

%\cite{Moch:2005uc}
\bibitem{Moch:2005uc}
  S.~Moch and P.~Uwer,
  %``XSummer: Transcendental functions and symbolic summation in Form,''
  Comput.\ Phys.\ Commun.\  {\bf 174} (2006) 759
  [arXiv:math-ph/0508008].
  %%CITATION = MATH-PH 0508008;%%

%\cite{Smirnov:1999gc}
\bibitem{Smirnov:1999gc}
  V.~A.~Smirnov,
  %``Analytical result for dimensionally regularized massless on-shell  double
  %box,''
  Phys.\ Lett.\ B {\bf 460} (1999) 397
  [arXiv:hep-ph/9905323].
  %%CITATION = HEP-PH 9905323;%%

%\cite{Tausk:1999vh}
\bibitem{Tausk:1999vh}
  J.~B.~Tausk,
  %``Non-planar massless two-loop Feynman diagrams with four on-shell legs,''
  Phys.\ Lett.\ B {\bf 469} (1999) 225
  [arXiv:hep-ph/9909506].
  %%CITATION = HEP-PH 9909506;%%

%\cite{Groote:2004qq}
\bibitem{Groote:2004qq}
  S.~Groote, J.~G.~K{\"o}rner and A.~A.~Pivovarov,
  %``Laurent series expansion of sunrise-type diagrams using configuration  space
  %techniques,''
  Eur.\ Phys.\ J.\ C {\bf 36} (2004) 471
  [arXiv:hep-ph/0403122].
  %%CITATION = HEP-PH 0403122;%%

%\cite{Czakon:2005rk}
\bibitem{Czakon:2005rk}
  M.~Czakon,
  %``Automatized analytic continuation of Mellin-Barnes integrals,''
  arXiv:hep-ph/0511200.
  %%CITATION = HEP-PH 0511200;%%

\end{thebibliography}

%- }}}

\end{document}

