%Title: Next-to-next-to-leading order vacuum polarization function of heavy quark near threshold and sum rules for $b \bar b$ system.
%Author: A.A.Penin and A.A.Pivovarov
%Published: * Phys. Lett. * ** B435 ** (1998) 413-419.
%hep-ph/9803363
\documentstyle[12pt]{article}
\renewcommand{\baselinestretch}{1.6}
\textwidth 6in
\textheight 9.0in
\topmargin = 0mm
\evensidemargin = 0mm
\oddsidemargin = 0mm
\begin{document}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\al}{\alpha}
\title
{
\vspace*{-20mm}
\bf
\begin{flushright}
{\large \bf TTP/98-13\\[-5mm]
MZ-TH/98-11\\[-5mm]
March 1998\\[15mm]}
\end{flushright}
\bf
Next-to-next-to-leading order vacuum polarization function of heavy quark
near threshold and sum rules for
$b \bar b$ system.
}
\author{
A.A.Penin\thanks{On leave from Institute for Nuclear Research,
Moscow, Russia}\\[-5mm]
{\small {\em Institut f{\"u}r Theoretische Teilchenphysik
Universit{\"a}t Karlsruhe}}\\[-5mm]
{\small {\em D-76128 Karlsruhe, Germany}}\\
and\\
A.A.Pivovarov{$^*$}\\[-5mm]
{\small {\em
Institut f\"ur Physik, Johannes-Gutenberg-Universit\"at}}\\[-5mm]
{\small {\em Staudinger Weg 7, D-55099 Mainz, Germany}}
}
\date{}
\maketitle
\begin{abstract}
A correlator of the vector current of a heavy quark
is computed analytically near threshold
in the next-to-next-to-leading order in perturbative and relativistic
expansion that
includes $\al_s^2$, $\al_sv$ and $v^2$ corrections
in the coupling constant and velocity of the heavy quark
to the nonrelativistic Coulomb approximation.
Based on this result, the numerical values of the $b$-quark pole mass
and the strong coupling constant are determined
from the analysis of sum rules for
the $\Upsilon$ system.\\[2mm]
PACS numbers: 11.65.Fy, 12.38.Bx, 12.38.Cy, 11.55.Hx
\end{abstract}
\thispagestyle{empty}
\newpage
\noindent
Insufficiency of the ordinary PT for description the near threshold
behavior of vacuum polarization function was noted long ago in the context
of Coulombic resummation in nonrelativistic QED \cite{Braun,Barb}.
Recently a considerable progress has been made in studying the
near threshold production of heavy quark-antiquark pair within
perturbation theory of QCD with resummation of threshold singularities.
Both perturbative and relativistic corrections have been taken into
account in the
next-to-next-to-leading order
in the coupling constant and velocity of the heavy quark
to the leading nonrelativistic approximation based on Coulomb
potential \cite{KPP,Hoang,Mel}.
This theoretical development provides more accurate
description of the heavy quark vacuum polarization
function in the threshold region necessary for
such applications as the top
quark production \cite{top} and the precise quantitative
investigation of the $\Upsilon$
system \cite{NSVZ,Vol2}. In the latter case
higher order corrections to leading Coulomb behavior in the threshold
region
are essential both
numerically for extracting the
$b$-quark mass and the strong coupling constant \cite{KPP}
and qualitatively for justifying the perturbative expansion
around Coulomb solution.
The analytical calculation of
the next-to-next-to-leading order corrections
has not
been completed yet though some results are available\footnote{
Semi-analytical analysis of the complete next-to-next-to-leading
order corrections to the heavy quark polarization function
near the two-particle threshold has been done
in the context of the photon mediated t quark pair production
\cite{Hoang,Mel}.}.
In this paper we present the complete analytical expression for
a correlator of the
vector current of heavy
quarks near threshold
in the next-to-next-to-leading order
resumming all ${\cal O}[(\al /v)^n\times(\al_s^2,\al_sv,v^2)]$
terms, with $v$ being the heavy quark velocity.
The correlator is further used
for determination of the bottom quark pole mass $m_b$
and the strong coupling constant $\al_s$ from sum rules for
the $\Upsilon$ system.
We study the near
threshold behavior of the polarization function $\Pi(s)$
of the $b$-quark vector current $j_\mu=\bar b\gamma_\mu b$
\[
\left(q_\mu q_\nu-g_{\mu \nu}q^2\right)\Pi(q^2)=
i\int dxe^{iqx}\langle 0|Tj_{\mu}(x)j_{\nu}(0)|0\rangle
\]
within the nonrelativistic expansion \cite{CasLep} which
in the next-to-next-to-leading order reads
\be
\Pi(s)
={N_c\over 2 m_b^2}\left(
C_h(\al_s)G(0,0,k)+{4\over 3}{k^2\over m_b^2}G_C(0,0,k)\right)
\label{thmom}
\ee
with $k=\sqrt{m_b^2-s/4}$ being a natural energy variable near
threshold.
First term in brackets gives the representation for the correlator
within NRQCD with
$C_h(\al_s)$ being a
perturbative coefficient matching
correlators of relativistic and nonrelativistic vector currents.
The coefficient
$C_h(\al_s)$ is computable in full QCD and by now is
known to the second order in $\al_s$ expansion
\[
C_h(\al_s)
=1-C_h^1C_F{\al_s\over \pi}+C_h^2C_F\left({\al_s\over \pi}\right)^2
\]
with
$
C_h^1=4$ \cite{KS}
and
\[
C_h^2=\left({39\over 4}-\zeta(3)+{4\pi^2\over 3}
\ln{2}-{35\pi^2\over 18}\right)C_F-\left({151\over 36}-{13\over 2}
\zeta(3)+{8\pi^2\over 3} \ln{2}-{179\pi^2\over 72}\right)C_A
\]
\be
+
\left({44\over 9}-{4\pi^2\over 9}+{11\over 9}n_f\right)T_F
+2\left(\beta_0+{\pi^2\over 3}C_F+{\pi^2\over 2}C_A
\right)\ln{\left(m_b\over\mu\right)}
\label{ch}
\ee
in $\overline{\rm MS}$ renormalization scheme \cite{ch}.
Here the group invariants for QCD are $C_A=3$, $C_F=4/3$, $T_F=1/2$,
and $\gamma_E=0.577216\ldots$ is the Euler constant, $\zeta(z)$ is
the Riemann $\zeta$-function, $n_f$ is the number of light flavors, and
$
\beta_0=11-2 n_f/3$.
The quantity $G({\bf x},{\bf y},k)$ is the nonrelativistic
Green function (GF) of the following
Schr{\"o}dinger equation
\[
\left(-{\Delta_{\bf x}\over m_b}-{\Delta_{\bf x}^2\over 4m_b^3}
+V_C(x)+{\al_s\over 4\pi}V_1(x)+\left({\al_s\over 4\pi}\right)^2
V_2(x)\right.
\]
\be
\left.
+V_{NA}(x)+V_{BF}({\bf x},{\bf s})
+ {k^2\over m_b}\right)
G({\bf x},{\bf y},k)=\delta({\bf x}-{\bf y})
\label{Schr}
\ee
where $V_C(x)=-C_F\al_s/x$ is the Coulomb
potential which is supposed to dominate the whole QCD interaction in the
energy region of interest, $x=|{\bf x}|$,
$V_{NA}(x)=-C_AC_F\al_s^2/(m_bx^2)$ is the non-Abelian potential
of quark-antiquark interaction \cite{Gup}, $V_{BF}({\bf x},{\bf s})$
is the standard Breit-Fermi potential (up to the color factor $C_F$)
containing the quark
spin operator ${\bf s}$, {\it e.g.}~\cite{Landau}.
The terms $V_i$ ($i=1,2$) represent
first and second order perturbative QCD corrections to the Coulomb
potential \cite{Fish,Peter}
$$
V_1(x)=V_C(x)(C_0^1+C_1^1\ln(x\mu)),
$$
\be
V_2(x)=V_C(x)(C_0^2+C_1^2\ln(x\mu)+C_2^2 \ln^2(x\mu)),
\label{potcorr1}
\ee
where
\[
C_0^1=a_1+2\beta_0\gamma_E,\qquad C_1^1=2\beta_0,
\]
\[
C_0^2=\left({\pi^2\over 3}+\gamma_E^2\right)\beta_0^2
+2(\beta_1+2\beta_0a_1)\gamma_E+a_2,
\]
\[
C_1^2=2(\beta_1+2\beta_0a_1)+8\beta_0^2\gamma_E,
\qquad
C_2^2=4\beta_0^2,
\]
$$
a_1={31\over 9}C_A-{20\over 9}T_Fn_f,
$$
\[
a_2= \left({4343\over 162}+6\pi^2-{\pi^4\over 4}
+{22\over3}\zeta(3)\right)C_A^2-
\left({1798\over 81} + {56\over 3}\zeta(3)\right)C_AT_Fn_f
\]
$$
-\left({55\over 3} - 16\zeta(3)\right)C_FT_Fn_f
+\left({20\over 9}T_Fn_f\right)^2,
$$
\[
\beta_1={34\over 3}C_A^2-{20\over 3}C_AT_Fn_f-4C_FT_Fn_f.
\]
The second term in eq.~(\ref{thmom}) is generated by the
operator of dimension four in
the nonrelativistic expansion of the
vector current (see, for example, \cite{Bod}). It contains
the GF of the pure Coulomb Schr{\"o}dinger equation \cite{Schw}
at the origin
\be
G_C(0,0,k)=-{C_F\al_sm_b^2\over 4\pi}\left({k\over C_F\al_sm_b}+
\ln\left({k\over \mu}\right)+\gamma_E+\Psi_1\left(1- {C_F\al_sm_b\over 2k}
\right)\right)
\label{G0}
\ee
where $\Psi_1(x)={\Gamma'(x)/\Gamma(x)}$
and $\Gamma(x)$ is the Euler $\Gamma$-function.
The solution to eq.~(\ref{Schr}) can be
found within the standard nonrelativistic
perturbation theory
around the Coulomb GF $G_C({\bf x},{\bf y},k)$.
The leading order
corrections to the Coulomb GF at the origin
due to
$\Delta^2$, $V_{NA}$ and $V_{BF}$ terms are known
analytically
\cite{Hoang,Mel}\footnote{The term $V_{NA}$ can be fully
accounted for the
Coulomb GF
because the corresponding differential
equation is exactly solvable in
standard special functions. Numerically this is not important for
applications though.}.
After including these corrections
the approximate GF of eq.~(\ref{Schr}) at the origin
takes the form \cite{Hoang}
\[
G(0,0,k)=-
{C_F\al_sm_b^2\over 4\pi}\left(
\left(1-{5\over 8}{k^2\over m_b^2}\right){k\over C_F\al_sm_b}+
\left(1-2{k^2\over m_b^2}\right)\right.
\]
\be
\left.\left(
\ln\left({k\over \mu}\right)+\gamma_E
+\Psi_1\left(1-
{C_F\al_sm_b\over 2k}\right)\right)
+{11\over 16}{C_F\al_s k\over m_b}
\Psi_2\left(1-{C_F\al_sm_b\over 2k}\right)
\right)
\label{Hoangcorr}
\ee
\[
+{4\pi\over 3}{C_F\al_s\over m_b^2}
\left(1+{3\over 2}{C_A\over C_F}\right)G_C^2(0,0,k)
\]
where $\Psi_2(x)=\Psi_1'(x)$.
Note that in ref.~\cite{Hoang} the shift of the
spectrum of intermediate nonrelativistic Coulomb bound states
was treated exactly {\it i.e.} without expanding of
the energy denominators. This accounts for
a part of the higher order corrections. We, however, consistently
work in the next-to-next-to-leading order and keep only
the second order terms in eq.~\ref{Hoangcorr}.
Since this part of the corrections is relatively small
the difference between these two approaches is
really negligible for the numerical analysis of the sum rules.
The correction $\Delta G_1$ to eq.~(\ref{Hoangcorr})
due to the first iteration
of $V_1$ term of the QCD potential
has been found in ref.~\cite{KPP} where
the consistent analysis of sum rules for $b\bar b$ system in
the next-to-leading order
has been performed
$$
\Delta G_1(0,0,k)={\al_s\over 4\pi}{C_F\al_sm_b^2\over 4\pi}\left(
\sum_{m=0}^\infty F^2(m)(m+1)
\left(C_0^1+(L_k+\Psi_1(m+2))\right)C_1^1
\right.
$$
\[
-2\sum_{m=1}^\infty\sum_{n=0}^{m-1}
F(m)F(n)
{n+1\over m-n}C_1^1
+2\sum_{m=0}^\infty F(m)
\left(C_0^1+(L_k -
2\gamma_E-\Psi_1(m+1))C_1^1\right)
\]
\be
\left.+2\left(-\gamma_E L_k+
{1\over 2}L_k^2\right)C_1^1
\right)
\label{G1}
\ee
where $L_k=\ln\left({\mu \over 2k}\right)$ and
\[
F(m)={C_F\al_s m_b\over (m+1)2k}\left(m+1-{\displaystyle {C_F\al_s m_b\over 2k}}\right)^{-1}.
\]
The correction $\Delta G_2^{(2)}$ to eq.~(\ref{Hoangcorr})
due to $V_2$ part of the potential
is also known \cite{KPP}
\[
\Delta G_2^{(2)}(0,0,k)=\left({\al_s\over 4\pi}\right)^2{C_F\al_sm_b^2\over
4\pi}\left( \sum_{m=0}^\infty F^2(m)
\left((m+1)\left(C_0^2+L_k C_1^2+L_k^2 C_2^2\right)
\right.\right.
\]
$$
\left.
+(m+1)\Psi_1(m+2)\left(C_1^2+2L_k
C_2^2\right)+I(m)C_2^2\right)
$$
\be
+2\sum_{m=1}^\infty\sum_{n=0}^{m-1}
F(m)F(n)\left(-
{n+1\over m-n}\left(C_1^2 +2L_k C_2^2\right)
+J(m,n)C_2^2\right)
\label{G22}
\ee
\[
+2\sum_{m=0}^\infty F(m)
\left(C_0^2+L_k C_1^2+
(L_k^2+K(m))C_2^2-(2\gamma_E+
\Psi_1(m+1))
\left(C_1^2+2L_k C_2^2\right)\right)
\]
$$
\left.+2\left(-\gamma_E L_k+
{1\over 2}L_k^2\right)
C_1^2
+N(k)C_2^2\right)
$$
where
\[
I(m)=(m+1)\left(\Psi^2_1(m+2)-\Psi_2(m+2)+2\zeta(2)-{2\over(m+1)^2}\right)
\]
$$
-2(\Psi_1(m+1)+\gamma_E),
$$
$$
J(m,n)= 2{n+1\over m-n}\left(\Psi_1(m-n)-{1\over n+1}+2\gamma_E\right)
$$
\[
+2{m+1\over m-n}(\Psi_1(m-n+1)-\Psi_1(m+1)),
\]
$$
K(m)=2(\Psi_1(m+1)+\gamma_E)^2+\Psi_2(m+1)-\Psi_1^2(m+1)+2\gamma_E^2,
$$
\[
N(k)=\left(\gamma_E+{\pi^2\over 6}\right)L_k
-\gamma_E L_k^2+{1\over 3}
L_k^3.
\]
In this paper we complete these results by computing
the correction $\Delta G_2^{(1)}$ due to the
second iteration of $V_1$ term which of the proper
(next-to-next-to-leading) order according to counting
of smallness in nonrelativistic QCD with respect to $\al_s$ and $v$.
The result reads
\[
\Delta G_2^{(1)}(0,0,k)=\left({\al_s\over 4\pi}\right)^2
{(C_F\al_s)^2\over 4\pi}{m_b^3\over 2k}
\left( \sum_{m=0}^\infty H^3(m)(m+1)
\right.
\]
$$
\left(C_0^1+
\left(\Psi(m+2)+
L_k\right)C_1^1\right)^2
$$
\[
-2\sum_{m=1}^\infty\sum_{n=0}^{m-1}{n+1\over m-n}C_1^1
\left(H^2(m)H(n)\left(C_0^1+\left(\Psi(m+2)+
L_k\right)C_1^1\right)
\right.
\]
\be
\left.
+H(m)H^2(n)\left(C_0^1+\left(\Psi(n+2)+
L_k\right)C_1^1\right)\right)
\label{G21}
\ee
\[
+2(C_1^1)^2\left(\sum_{m=2}^\infty\sum_{l=1}^{m-1}\sum_{n=0}^{l-1}
{H(m)H(n)H(l)}{n+1\over (l-n) (m-n)}\right.
\]
$$
+\sum_{m=2}^\infty\sum_{n=1}^{m-1}\sum_{l=0}^{n-1}
{H(m)H(n)H(l)}{l+1\over (n-l)(m-n)}
$$
\[
\left.\left.
+\sum_{n=2}^\infty\sum_{m=1}^{n-1}\sum_{l=0}^{m-1}
{H(m)H(n)H(l)}{(l+1)(m+1)\over (n+1)(n-l)(n-m)}\right)\right)
\]
where
$$
H(m)=\left(m+1-{\displaystyle {C_F\al_s m_b\over 2k}}\right)^{-1}.
$$
We are going to describe the details of
this rather cumbersome
calculation
elsewhere. One remark is in order though.
Because
ultraviolet divergences in eq.~(\ref{Hoangcorr})
depend on $k$ one has
to match the calculation of these corrections to the
calculation of the Wilson coefficient $C_h(\al_s)$ (eq.~(\ref{ch}))
\cite{Hoang,Mel}.
Such matching is not necessary for the calculation of $\Delta G_1$ and
$\Delta G_2^{(i)}$
terms because their divergent parts are $k$ independent.
Thus
eqs.~(\ref{thmom},~\ref{G0}-\ref{G21}) give
the complete analytical expressions for the vacuum polarization
function of heavy quarks near the two-particle
threshold in the next-to-next-to-leading order\footnote{In
refs.~\cite{Hoang,Mel} the corrections to the
Coulomb GF due to $V_i$ ($i=1,2$)
terms of the potential were treated numerically for complex values
of energy far from the real axis.}. Eqs. (\ref{G0}-\ref{G21}) look
awkward and they can be rendered into more readable form by using
$\Psi$ functions for expressing some of the sums entering the formulae.
However for direct numerical analysis of sum rules for $b\bar b$
system this form is most suitable with respect to applicability of
efficient numerical algorithms of a symbolic system.
Obtained formulae are applied to
the analysis of the
$\Upsilon$ system for extraction of the
$b$-quark pole mass $m_b$ and the coupling
constant $\al_s$.
The sum rules are formulated in the literature \cite{KPP,Vol2}
and we will use the latest version \cite{KPP} with correct large $n$
behavior.
The moments ${\cal M}_n$
\[
{\cal M}_n=
\left.{12\pi^2\over n!}(4m_b^2)^n{d^n\over ds^n}\Pi(s)\right|_{s=0}=
(4m_b^2)^n\int_0^\infty{R(s)ds\over s^{n+1}}
\]
of the spectral density
$R(s)=12\pi {\rm Im}\Pi(s+i\epsilon)$
are compared with experimental ones
\[
{\cal M}_n^{exp} =
%(4m_b^2)^n\int_0^\infty{R(s)ds\over s^{n+1}}=
{(4m_b^2)^n\over Q_b^2}\int_0^\infty{R_b(s)ds\over s^{n+1}}
\]
under the assumption of quark-hadron duality.
The experimental moments ${\cal M}_n^{exp}$ are generated by
the function $R_b(s)$ which is the
normalizad cross section
$R_b(s) = \sigma(e^+e^-\rightarrow {\rm hadrons}_{\,b\bar b})/
\sigma(e^+e^-\rightarrow \mu^+\mu^-)$.
Here $Q_b=-1/3$ is the $b$-quark electric charge.
Numerical values
are obtained basically by saturating
the experimental moments with the contribution
of
the first six $\Upsilon$ resonances (see \cite{KPP} for details).
Their leptonic widths
$\Gamma_{k}$ and masses $M_{k}$ $(k=1\ldots 6)$
are known with good accuracy \cite{PDG}
\[
{\cal M}_n^{exp}={(4m_b^2)^n\over Q_b^2}
\left({9\pi\over \al_{QED}^2(m_b)}
\sum_{k=1}^6{\Gamma_{k}\over M_{k}^{2n+1}}
+\int_{s_0}^\infty\!{\rm d}s {R_b(s)\over s^{n+1}}\right).
\]
The rest of the spectrum beyond the resonance
region for energies larger than
$s_0\approx (11.2~{\rm GeV})^2$ (continuum contribution)
lies far from threshold and is safely
approximated by
the ordinary PT expression for the theoretical
spectral density, so there $R_b(s)\approx R(s)$.
The influence of the continuum on high moments
is almost negligible numerically
and in any case under strict control\footnote{The expressions
for the first few moments of the spectral density
are now available
in ordinary perturbation theory
with $\al_s^2$ accuracy \cite{Chet}, however, they cannot
be used in theoretical formulas for sum rules directly
because the spectrum is well known
experimentally only for energies close
to threshold due to existence of sharp resonances while the
contribution of the continuum to these low moments is large in
comparison with the resonance contribution.}.
Electromagnetic coupling constant is renormalized to the
energy of order
of $m_b$ with the result $\al_{QED}^2(m_b)=1.07 \al^2$ \cite{PDG}.
We work with moments for $10