%Title: $M(B^*_c)M(B_c)$ Splitting from Nonrelativistic Renormalization Group
%Author: A.A. Penin, A. Pineda, V.A. Smirnov, M. Steinhauser
%Published: Phys.Lett. B593 (2004) 124134.
%hepph/0403080
\documentclass[12pt]{article}
%\usepackage{axodraw,a4wide,epsfig}
\usepackage{a4wide,epsfig}
\voffset0cm
\hoffset0cm
\oddsidemargin0cm
\evensidemargin0cm
\topmargin0cm
\textwidth16.5cm
\textheight22cm
\setlength{\arraycolsep}{0.5mm}
\newcommand{\logminus}{\ln\left(2z^{\beta_0}\right)}
\newcommand{\betnull}{\beta_0}
\newcommand{\mrM}{\mu_r}
\newcommand{\arccot}{\mathop{\mbox{arccot}}\nolimits}
\newcommand{\bfm}[1]{\mbox{\boldmath$#1$}}
\newcommand{\bff}[1]{\mbox{\scriptsize\boldmath${#1}$}}
\def\als{\alpha_{\rm s}}
\newcommand{\cf}{C_F}
\renewcommand{\textfraction}{0}
\renewcommand{\topfraction}{1}
\renewcommand{\bottomfraction}{1}
\newcommand{\nn}{\nonumber}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\renewcommand{\textfraction}{0}
\renewcommand{\topfraction}{1}
\renewcommand{\bottomfraction}{1}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\title{\vskip3cm{\baselineskip14pt
\centerline{\normalsize DESY 04042\hfill}
\centerline{\normalsize TTP0406 \hfill}
\centerline{\normalsize UBECMPF0405 \hfill}
%\centerline{\normalsize hepph/0403080\hfill}
%\centerline{\normalsize January 2004\hfill}
}
\vskip1.5cm
$M(B^*_c)M(B_c)$ Splitting from Nonrelativistic Renormalization Group}
\author{A.A. Penin$^{a,b}$, A. Pineda$^c$,
V.A. Smirnov$^{d,e}$, M. Steinhauser$^e$\\
{\small\it $^a$ Institut f{\"u}r Theoretische Teilchenphysik,
Universit{\"a}t Karlsruhe, 76128 Karlsruhe, Germany}\\
{\small\it $^b$ Institute for Nuclear Research, Russian Academy of Sciences,
%60th October Anniversary Prospect 7a,
117312 Moscow, Russia}\\
{\small\it $^c$ Dept. d'Estructura i Constituents de la Mat\`eria,
U. Barcelona,
%Diagonal 647,
E08028 Barcelona,
%Catalonia,
Spain}\\
{\small\it $^d$ Institute for Nuclear Physics, Moscow State
University, 119992 Moscow, Russia}\\
{\small\it $^e$ II. Institut f\"ur Theoretische Physik, Universit\"at Hamburg,
%Luruper Chaussee 149,
22761 Hamburg, Germany}}
\date{}
\maketitle
\thispagestyle{empty}
\begin{abstract}
We compute the hyperfine splitting in a heavy quarkonium composed of different
flavors in nexttoleading logarithmic approximation using the nonrelativistic
renormalization group. We predict the mass difference of the vector and
pseudoscalar charmbottom mesons to be $M(B^*_c)M(B_c)=65 \pm 24\,{(\rm
th)}\,{}^{+19}_{16}\,(\delta\alpha_s)$~MeV.
\medskip
\noindent
PACS numbers: 12.38.Bx, 14.65.Fy, 14.65.Ha
\end{abstract}
%\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
The recently discovered charmbottom heavy quarkonium completes the well
investigated charmonium and bottomonium families and offers a new perspective
in the study of the nonrelativistic dynamics of the strong interactions. The
first experimental observation of about twenty events interpreted as the
decays of the $B_c$ meson by CDF collaboration \cite{Abe} does not match the
precision of the spin one charmonium and bottomonium measurements. The
statistics, however, is expected to increase significantly in future
experiments at Tevatron and LHC greatly improving the accuracy of the
data. Note that only the pseudoscalar (spin singlet) state has been observed
so far while the vector (spin triplet) meson $B_c^*$ is still to be
discovered. This distinguishes $c\bar b$ quarkonium from the $b\bar b$ system,
where it is the pseudoscalar $\eta_b$ meson, which asks for experimental
detection.
From the theoretical point of view, the charmbottom mesons are ``in between''
the approximately Coulomb $b\bar b$ mesons and the $c\bar c$
mesons. Therefore, a simultaneous analysis of all three quarkonia could shed
new light on the balance between the perturbative and nonperturbative effects
and further check whether a perturbative analysis provides a reliable starting
point for them. Moreover, since the nonperturbative effects in the $c\bar b$
system are suppressed with respect to the $c\bar c$ meson, the former could be
a cleaner place to determine the charm quark mass (provided the experimental
accuracy is good enough). Another point to be stressed is that, though the
leading order dynamics of the $c\bar b$ state is quite similar to the $b\bar
b$ and $c\bar c$ one (up to the value of the reduced mass) the higher order
relativistic and perturbative corrections are different. Thus the comparison
of $c\bar b$ and $b\bar b$ ($c\bar c$) properties could help to establish fine
details of the nonrelativistic dynamics.
The spectrum of the charmbottom quarkonium has been subject of numerous
investigations based on potential models
\cite{Eichten:1994gt,Gershtein:1994jw}, lattice simulations \cite{Sha}, and
pNRQCD \cite{BraVai}. This last analysis computed the ground state energy
within a pure perturbative approach. We consider that this analysis further
indicates that a perturbative approach can be a good starting point for
studing the $B_c$ system.
In the present paper we focus on the hyperfine splitting (HFS) $E_{\rm hfs}$
of the $B_c$, {\it i.e.} the mass difference between the singlet and triplet
spin states $M(B^*_c)M(B_c)$. The QCD study of the heavy quarkonium HFS has a
long history \cite{BucNg,recent}. For the sameflavor quarkonium the
nexttoleading order (NLO) ${\cal O}(\alpha_s)$ correction to the ground
state HFS can be found in \cite{PenSte} in a closed analytical form. The
leading order HFS is proportional to the fourth power of the strong coupling
constant $\alpha_s(\nu)$ and thus the low order calculations suffer from
strong spurious dependence on the renormalization scale $\nu$, which
essentially limits the numerical accuracy of the approximation. Hence, the
proper fixing of the normalization scale becomes mandatory for the HFS
phenomenology. The dynamics of the nonrelativistic bound state, however, is
characterized by three well separated scales: the hard scale of the heavy
quark mass $m$, the soft scale of the bound state momentum $mv$, and the
ultrasoft scale of the bound state energy $mv^2$, where $v\propto\alpha_s$ is
the velocity of the heavy quark inside the approximately Coulomb bound state.
To make the procedure of scale fixing selfconsistent one has to resum to all
orders the large logarithms of the scale ratios. For the sameflavor case
this problem has been solved in Ref.~\cite{KPPSS} within the nonrelativistic
renormalization group (NRG) approach and the nexttoleading logarithmic (NLL)
result for HFS has been derived. The renormalization group improved result
shows better stability with respect to the scale variation. Moreover, the use
of the NRG significantly improves the agreement with the experimental value of
HFS in charmonium in comparison to the NLO computation. Below we generalize
the analysis to the differentflavor quarkonium case and apply the result to
predict the splitting $M(B^*_c)M(B_c)$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Renormalization group running of the spin dependent potential}
To derive the NRG equations necessary for the NLL analysis of the HFS, we rely
on the method based on the formulation of the nonrelativistic effective theory
\cite{CasLep} known as potential NRQCD (pNRQCD) \cite{PinSot1}. The method was
developed in Ref.~\cite{Pin} where, in particular, the leading logarithmic
(LL) result for HFS has been obtained (see also Ref.~\cite{HMS}). In pNRQCD
the HFS is generated by the spinflip potential in the effective Hamiltonian,
which in momentum space has the form
\be
\delta {\cal H}_{\rm spin} = D_{S^2,s}^{(2)}{4C_F\pi\over 3m_1m_2}\bfm{S}^2,
\qquad {\bfm S}={{\bfm \sigma}_1+{\bfm \sigma}_2\over 2}\,,
\label{pot}
\ee
where ${\bfm \sigma}_{1}$ and ${\bfm \sigma}_{2}$ are the spin operators of
the quark and antiquark with masses $m_1$ and $m_2$,
$C_F=(N_c^21)/(2N_c)$, and $D_{S^2,s}^{(2)}$ is the Wilson coefficient, which
incorporates the effects of the modes that have been integrated out. In
effective theory calculations such couplings become singular as a result of
the scale separation. The renormalization of these singularities allows one
to derive the equations of the NRG,
which describe the running of the effectivetheory couplings, {\it i.e.} their
dependence on the effectivetheory cutoffs. The solution of these equations
sums up the logarithms of the scale ratios.
In general, one should consider
the soft, potential and ultrasoft running of $D_{S^2,s}^{(2)}$
corresponding to the ultraviolet
divergences of the soft, potential, and ultrasoft regions \cite{BenSmi}. We
denote the corresponding cutoffs as $\nu_s$, $\nu_p$ and $\nu_{us}$,
respectively. $\nu_{us}$ and $\nu_{p}$ are correlated as was first realized in Ref. \cite{LMR}.
A natural relation between them is
$\nu_{us}=\nu_p^2/(2m_r)$, where $m_r=m_1m_2/(m_1+m_2)$ is the reduced mass.
The dependence on $\nu_s$ first emerges
in the LL approximation after integrating out the hard modes.
It disappears after subsequent integrating out the soft modes
giving rise to a dependence on $k$, the
threedimensional momentum transfer between the quark and antiquark.
Thus the soft running effectively stops at $\nu_s=k$.
The dependence on $\nu_p$ emerges for the first time in the NLL approximation
and cancels out in the timeindependent Schr\"odinger perturbation
theory for heavy quarkonium observables. Thus, in pNRQCD one considers
$D_{S^2,s}^{(2)}$ as a function of $k$ and $\nu_p$.
For the calculation of the spectrum it is convenient to expand this
$k$dependent potential around $k=\nu_s$
\be
D_{S^2,s}^{(2)}(k,\nu_p)=D_{S^2,s}^{(2)}(\nu_s,\nu_p)+
\ln\left({k\over\nu_s}\right)
\nu_s{d \over d\nu_s} D_{S^2,s}^{(2)}(\nu_s,\nu_p)+\ldots
\,.
\label{momdep}
\ee
The characteristic momentum for the Coulomb system
is $ \alpha_sm_r$ and for $\nu_s\sim \alpha_sm_r$
the average of $\ln\left({k/\nu_s}\right)$
over bound state wave function does not produce a large logarithm
while the derivative in $\ln \nu_s$ results in extra factor of $\alpha_s$.
Thus, for the calculation of the HFS in NLL approximation
one can take the first term on the right hand side of
Eq.~(\ref{momdep}) in the NLL approximation, the second term in the LL
approximation and neglect the higher derivative terms.
Once expanded, the potential is a function of $\nu_s$ and $\nu_p$ (we should
not forget that there is also a dependence on $m_i$, the masses of the heavy
quarks, and $\nu_h$, the matching scale of the order of the heavy quark
masses). Let us start with the discussion of the soft running. To the NLL
approximation it is determined by the following NRG equation
\begin{eqnarray}
\nu_s {d\over d\nu_s}D_{S^2,s}^{(2)}
&=&
\alpha_s
c_F(m_1)c_F(m_2)\gamma_s
\label{Dsoftrun}
\,,
\end{eqnarray}
where $c_F$ is the effective Fermi coupling,
\begin{eqnarray}
\gamma_s&=&\gamma^{(1)}_s{\alpha_s \over
\pi}+
\gamma^{(2)}_s{\alpha_s^2 \over \pi^2}+\cdots
\end{eqnarray}
is the soft anomalous dimension and $\alpha_s = \alpha_s(\nu_s)$ is
renormalized in the $\overline{\rm MS}$ scheme. The running of the
coefficient $c_F$ is known in NLL
approximation \cite{ABN}. It reads
\bea
c_F(m_i)&=&
z^{{\gamma_0\over 2}}
\left[ 1 + \frac{\alpha_s(\nu_h)}{4\pi}
\left(c_1+\frac{\gamma_0}{2}\ln\frac{\nu_h^2}{m_i^2}\right)
% \right.\nn\\&&\left.\mbox{}
+ \frac{\alpha_s(\nu_h)  \alpha_s(\nu_s)}{4\pi}\left(
\frac{\gamma_1}{2\beta_0}  \frac{\gamma_0\beta_1}{2\beta_0^2}
\right) + \dots \right]
\,,
\nonumber\\
\eea
where $z=\left(\alpha_s(\nu_s)/\alpha_s(\nu_h)\right)^{1/\beta_0}$,
$\nu_h\sim m_i$ is the hard matching scale, $c_1 = 2(C_A+C_F)$
and the one and twoloop anomalous dimensions read~\cite{ABN}
\be
\gamma_0 = 2 C_A \,, \qquad
\gamma_1 = \frac{68}{9}\,C_A^2  \frac{52}{9}\,C_A T_F\,n_l
\,.
\ee
Here $C_A=N_c$, $T_F=1/2$,
$n_l$ is the number of massless quark flavors,
and $\beta_i$ is the $(i+1)$loop
coefficient of the QCD $\beta$ function
\be
\beta_0 = \frac{11}{3}C_A\frac{4}{3} T_Fn_l\,, \qquad
\beta_1 = \frac{34}{3}C_A^2\frac{20}{3} C_A T_Fn_l  4C_FT_Fn_l\,.
\ee
The value of oneloop anomalous dimension
\begin{eqnarray}
\gamma^{(1)}_s &=& \frac{\beta_0}{2} + \frac{7}{4}C_A
\end{eqnarray}
can be extracted from the result of Ref.~\cite{Pin}.
The result for the twoloop coefficient
\begin{eqnarray}
\gamma^{(2)}_s&=&
\frac{1}{216}
\left[
{{C_A}}^2\,\left( 5  36\,{\pi }^2 \right) +
88\,{C_A}\,{n_l}\,{T_F} +
4\,{n_l}\,{T_F}\,\left( 27\,{C_F} 
40\,{n_l}\,{T_F} \right)
\right]
\,,
\label{gamma2}
\end{eqnarray}
is new. It was obtained by an explicit calculation of the subleading
singularities of the twoloop soft diagrams using the approach of
\cite{PinSot2,CMY,KPSS}. In this approach, dimensional regularization with
$D=42\varepsilon$ is used to handle the divergences, and the formal
expressions derived from the Feynman rules of the effective theory are
understood in the sense of the threshold expansion \cite{BenSmi}. Thus the
practical calculation reduces to the evaluation of the coefficients of the
quadratic and linear soft poles in $\varepsilon$. Our approach possesses two
crucial virtues: the absence of additional regulator scales and the automatic
matching of the contributions from different scales. For the reduction of the
twoloop Feynman integrals to the master ones the method of Ref.~\cite{Bai}
was used.
The solution of Eq.~(\ref{Dsoftrun}) can be written as a sum of the LL
and NLL contributions.
The LL result is already known and reads \cite{Pin}
(see also \cite{HMS})
\be
\left(D_{S^2,s}^{(2)}\right)^{LL}=\alpha_s(\nu_h)
\left[
1
+{2\beta_07C_A \over 2\beta_04C_A}\left(
z^{2C_A+\beta_0} 1 \right) \right]\,.
\label{DLL}
\ee
For the NLL term we obtain
\be
\left(\delta D_{S^2,s}^{(2)}\right)^{NLL}_{s}
=B_1\alpha_s^2(\nu_h)(z^{\gamma_0+\beta_0}1)+
B_2\alpha_s^2(\nu_h)(z^{\gamma_0+2\beta_0}1)
\,,
\label{DNLLs}
\ee
where
\bea
B_1&=& \frac{
{\beta_1}{\gamma_0} 
2 \beta_0^2
\left[c_1
+{\gamma_0 \over 2}\ln\left(\frac{\nu_h^2}{m_1m_2}\right)\right]
 \beta_0\gamma_1}{2
{{\beta_0}}^2
\left( {\beta_0}  {\gamma_0} \right)
\pi } {\gamma^{(1)}_s}
\,,
\\
B_2&=& \frac{
 {\beta_1}\,{\gamma_0}\,
{\gamma^{(1)}_s} +
{\beta_0}\,{\gamma_1}\,
{\gamma^{(1)}_s} +
{{\beta_0}}\,
\left( {\beta_1}\,{\gamma^{(1)}_s} 
4\,\beta_0\,{\gamma^{(2)}_s} \right) }{2\,
{{\beta_0}}^2\,
\left( 2\,{\beta_0}  {\gamma_0} \right)
\,\pi }
\,.
\eea
The potential running starts to contribute in NLL order.
To compute it we inspect all operators that lead to spindependent
ultraviolet divergences in the timeindependent perturbation theory
contribution with one and two potential loops \cite{Pin,KniPen2,Hil}.
They are
\begin{itemize}
\item[(i)] the ${\cal O}(v^2,\alpha_sv)$ operators
\cite{BucNg},
\item[(ii)] the tree ${\cal
O}(v^4)$ operators, some of which can be checked against the QED analysis
\cite{CMY,Pac},
\item[(iii)] the oneloop ${\cal O}(\alpha_sv^3)$ operators for which
only the Abelian parts are known \cite{CMY}, while the nonAbelian parts
are new.
\end{itemize}
In the NLL approximation, we need the LL soft and ultrasoft running of the
${\cal O}(v^2)$ and ${\cal O}(v^4)$ operators, which enter the twoloop
timeindependent perturbation theory diagrams, and the NLL soft and ultrasoft
running of the ${\cal O}(\alpha_sv)$ and ${\cal O}(\alpha_sv^3)$ operators,
which contribute at one loop. The running of the ${\cal O}(v^2,\alpha_sv)$
operators is already known within pNRQCD \cite{Pin}. The running of the other
operators is new. For some of them, it can be obtained using the
reparameterization invariance \cite{Manohar}.
We refrain from writing the corresponding system of NRG equations, which
is rather lengthy, and only present its solution, which can be cast in the form
\bea
\left(\delta D_{S^2,s}^{(2)}\right)^{NLL}_{p}&=&
\pi \alpha_s^2(\nu_h) \sum_{i=1}^{18} A_i f_i
\,,
\label{DNLLp}
\eea
where the coefficients $A_i$ and $f_i$ are given in the Appendix.
To get this result we rescale the ultrasoft cutoff to
$\nu_{us}=\nu_p^2/\nu_h$. The difference to the
previous definition is beyond the NLL accuracy.
The LL result~(\ref{DLL}) obeys the tree level matching condition
\be
\left.\left(D_{S^2,s}^{(2)}\right)^{LL}\right_{\nu=\nu_h}=
\alpha_s(\nu_h) \,,
\ee
while Eqs.~(\ref{DNLLs}) and~(\ref{DNLLp}) vanish at $\nu=\nu_h$ by
construction.
We then use the known oneloop result of the potential \cite{BucNg} to obtain the NLO
matching condition at the scale $k=\nu_s=\nu_p=\nu_h$. It reads
\bea
\left(D_{S^2,s}^{(2)}\right)_{\rm 1loop}&=&
\left[
{5 \over 9}T_Fn_l
\frac{5}{36}C_A +C_F+{7 \over 8}C_A\ln\left(\frac{\nu_h^2}{m_1m_2}\right)
\right.\nn\\&&
\left.
\frac{3}{4}\left(
C_F\frac{m_1m_2}{m_1+m_2}
+\frac{1}{2}\left(C_A2C_F\right)\frac{m_1+m_2}{m_1m_2}
\right)\ln\left(\frac{m_2}{m_1}\right)
\right]
{\alpha_s^2(\nu_h) \over \pi}
\,.
\nn\\
\label{DNLL1l}
\eea
Note that in the limit $m_1=m_2\equiv m_q$ this equation does not reproduce
the sameflavor equalmass expression
\bea
\left(D_{S^2,s}^{(2)}\right)^{q\bar q}_{\rm 1loop}&=&
\left[
{5 \over 9}T_Fn_l+{3 \over 2}(1\ln 2)T_F+{11C_A9C_F \over 18}
\right.\left.
+{7 \over 4}C_A\ln\frac{\nu_h}{m_q}
\right]
{\alpha_s^2(\nu_h) \over \pi}
\,,\nonumber
\\
\eea
because of the twogluon annihilation contribution present in the latter
case.
Thus the NLL approximation for the Wilson coefficient is given by the sum
\be
\left(D_{S^2,s}^{(2)}(\nu)\right)^{NLL}=\left(D_{S^2,s}^{(2)}(\nu)\right)^{LL}+\left(\delta
D_{S^2,s}^{(2)}(\nu)\right)^{NLL}_{s}+\left(\delta
D_{S^2,s}^{(2)}(\nu)\right)^{NLL}_{p}+\left(D_{S^2,s}^{(2)}\right)^{~}_{\rm
1loop}
\,.
\label{Dsum}
\ee
where $D_{S^2,s}^{(2)}(\nu)\equiv D_{S^2,s}^{(2)}(\nu,\nu)$ and we combine the soft and potential
running by setting $\nu_s=\nu_p=\nu$,
which is consistent at the order of interest.
From Eqs.~(\ref{pot}) and~(\ref{momdep})
we obtain the final result for the NLL spinflip potential
\be
\delta {\cal H}_{\rm spin} =\left[\left(D_{S^2,s}^{(2)}(\nu)\right)^{NLL}+
{\gamma_s^{(1)}\over\pi}\left(\alpha_s^2c_F^2\right)^{LL}\ln\left({k\over\nu}\right)\right]
{4C_F\pi\over 3m_1m_2}\bfm{S}^2
\,.
\label{finpot}
\ee
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Hyperfine splitting in NLL approximation}
We are now in the position to derive the
NLL result for the HFS. It is obtained by computing the
corrections to the energy levels with the insertion of the
potential~(\ref{finpot})
in the quantum mechanical perturbation theory.
The result for principal quantum number $n$ reads
\bea
E_{n,{\rm hfs}}^{NLL} &=&{4 \over 3} { C_F^2 \alpha_s\over n} E_n^C
\left\{
(1+2\delta \phi_n)
\left(D_{S^2,s}^{(2)}(\nu)\right)^{LL}
\right.
\nn\\&&
+\left(\ln{\left({n\nu\over \bar\nu }\right)}
+\Psi_1(n+1)+\gamma_E+{n1 \over 2n}
\right){\gamma_s^{(1)}\over\pi}\left(\alpha_s^2c_F^2\right)^{LL}
\nn\\&&
+\left.
\left(\delta D_{S^2,s}^{(2)}(\nu)\right)^{NLL}_{s}
+ \left(\delta D_{S^2,s}^{(2)}(\nu)\right)^{NLL}_{p}
+\left(D_{S^2,s}^{(2)}\right)^{NLL}_{\rm 1loop}
\right\}
\,,
\label{HFS}
\eea
where $\bar\nu=2C_F\alpha_sm_r$,
$E_n^C=  C_F^2\alpha_s^2m_r/(4n^2)$, $\Psi_n(z) =d^n \ln \Gamma
(z)/dz^n$, $\Gamma (z)$ is the
Euler $\Gamma$function, and $\gamma_E=0.577216\ldots$ is Euler's constant.
In Eq.~(\ref{HFS}) the first order correction to the Coulomb wave
function at the origin due to oneloop contribution to the static
potential reads \cite{KPP}
\be
\delta \phi_n ={\alpha_s \over \pi}
\left[{3\over 8}a_1+
{\beta_0\over 4}
\left(
3\ln{\left({n\nu\over \bar\nu }\right)}
+\Psi_1(n+1)2n\Psi_2(n)1+\gamma_E+{2\over n}
\right)
\right]
\,,
\ee
where $a_1=31C_A/920T_Fn_l/9$. Furthermore, the second line
of Eq.~(\ref{HFS}) results from the second term in square brackets
in Eq.~(\ref{finpot}) after average over the Coulomb wave
function.
By expanding the
resummed expression up to ${\cal O}(\alpha_s^2)$, we get
\bea
&&
E_{n,{\rm hfs}}^{NLL} ={4 \over 3} { C_F^2 \alpha_s^2\over n} E_n^C
\left\{
1 + {{\als \over \pi}}\,\left[ {C_F} + \frac{7\,{C_A}\,{L^n_{\alpha_s}}}{4}
+
\frac{7\,{C_A}}{8} \,\ln\left(\frac{4\,m_r^2}{ m_1 m_2 } \right)
\right.\right.\nn\\&&
+
\left( \frac{3\,{C_F}\,{m_r}}{ {m_1}  {m_2} } +
\frac{3\,{C_A}\,\left( {m_1} + {m_2} \right) }{8\,\left( {m_1}  {m_2} \right) } \right) \,
\ln\left(\frac{{m_1}}{{m_2}} \right)
+
\frac{{n_f}\,{T_F}\,
\left( 15  11\,{n} +
12\,{{n}}^2\,\Psi_2(n) \right) }{9\,{n}}
\nn
\\
&&
\left.

\frac{{C_A}\,\left( 393  41\,{n} 
126\,{\gamma_E}\,{n} 
126\,{n}\,\Psi_1(n) +
264\,{{n}}^2\,\Psi_2(n) \right) }{72\,{n}}
\right]
\nn
\\
&&
+ {{{\als^2 \over \pi^2}}}{{L^n_{\alpha_s}}}\left[ L^n_{\alpha_s}
\left( \frac{19\,{{C_A}}^2}{6} 
\frac{5\,{C_A}\,{n_f}\,{T_F}}{6} \right)
\right.\nn\\&&
+
\left( \frac{{{C_A}}^2}{6} 
\frac{11\,{C_A}\,{C_F}}{8} 
\frac{{{C_F}}^2\,\left( {{m_1}}^2 + {{m_2}}^2 \right) }{{\left( {m_1} + {m_2} \right) }^2} \right) \,{\pi }^2 
\frac{2\,{C_F}\,{n_f}\,{T_F}}{3}
\nn\\&&
+
\left( \frac{11\,{{C_A}}^2\,\left( {m_1} + {m_2} \right) }{8\,\left( {m_1}  {m_2} \right) } +
\frac{4\,{C_F}\,{n_f}\,{T_F}\,{m_r}}{{m_1}  {m_2} } +
{C_A}\,\left( \frac{11\,{C_F}\,{m_r}}{{m_1}  {m_2} }
\right.\right.\nn\\&&\left.\left.

\frac{{n_f}\,{T_F}\,\left( {m_1} + {m_2} \right)}{2\,\left( {m_1}  {m_2} \right) } \right) \right) \,
\ln\left(\frac{{m_1}}{{m_2}} \right)
+
\left( \frac{19\,{{C_A}}^2}{6} 
\frac{5\,{C_A}\,{n_f}\,{T_F}}{6} \right) \,
\ln\left(\frac{4\,{m_r^2}}{m_1 m_2}\right)
\nn\\&&

\frac{{{C_A}}^2\,\left( 1380  305\,{n} 
450\,{\gamma_E}\,{n} 
450\,{n}\,\Psi_1(n) +
924\,{{n}}^2\,\Psi_2(n) \right) }{144\,
{n}}
\nn\\&&\left.\left.
+ {C_A}\,
\left( \frac{43\,{C_F}}{12} +
\frac{{n_f}\,{T_F}\,
\left( 114  109\,{n}  18\,{\gamma_E}\,{n} 
18\,{n}\,\Psi_1(n) +
84\,{{n}}^2\,\Psi_2(n) \right) }{36\,
{n}} \right) \right]
\right\}
\label{ser}
\,,
\eea
where $\alpha_s\equiv \alpha_s(\nu)$, $\Psi_n(x)=d^n\ln\Gamma(x)/dx^n$,
$L^n_{\alpha_s}=\ln\left(C_F\alpha_s/n\right)$ and $\nu_h=2m_r$ and
$\nu=\bar{\nu}/n$ has been chosen. The ${\cal
O}(\alpha_s^2\ln^2\alpha_s)$ term is known \cite{Pin,HMS}, while the ${\cal
O}(\alpha_s^2\ln\alpha_s)$ term is new.
The equalmass case expression \cite{KPPSS}, relevant for charmonium and
bottomonium, can be deduced from Eq. (\ref{HFS}) by replacing
\be
\left(D_{S^2,s}^{(2)}\right)_{\rm 1loop}\to
\left(D_{S^2,s}^{(2)}\right)^{q\bar q}_{\rm 1loop}
\ee
and setting $m_1=m_2$. After including the onephoton annihilation
contribution, the Abelian part of the equalmass
result
reproduces the ${\cal O}(m\alpha_s^6\ln\alpha_s)$ and
${\cal O}(m\alpha_s^7\ln^2\alpha_s)$ corrections to the positronium
HFS (see {\it e.g.} \cite{Pac,CMY}).
\section{Numerical estimates and conclusions}
For the numerical estimates, we adopt the strategy of \cite{KPPSS} and replace
the onshell mass of the charm and bottom quarks by one half of the physical
masses of the ground state of bottomonium and charmonium~\cite{Hag}. In
practice, we take $m_b=4.73$ GeV and $m_c=1.5$ GeV, consistent with the
accuracy of our computation. Furthermore, we take $\alpha_s(M_Z)$ as an input
and run\footnote{For the running and decoupling of $\alpha_s$ we
use the program {\tt RunDec}~\cite{Chetyrkin:2000yt}.}
with fourloop accuracy down to the matching scale $\nu_h$ to ensure
the best precision. Below the matching scale the running of $\alpha_s$ is
used according to the logarithmic precision of the calculation in order not to
include nexttonexttoleading logarithms in our analysis. In
Fig.~\ref{fig1}, the HFS for the charmbottom quarkonium ground state is
plotted as a function of $\nu$ in the LO, NLO, LL, and NLL approximations for
the hard matching scale value $\nu_h= 1.95$~GeV. As we see, the LL curve
shows a weaker scale dependence compared to the LO one. The scale dependence
of the NLO and NLL expressions is further reduced, and, moreover, the NLL
approximation remains stable at the physically motivated scale of the inverse
Bohr radius, $C_F\alpha_s m_r\sim 0.9$~GeV, where the fixedorder expansion
breaks down. At the scale $\nu'\approx 0.85$~GeV, which is close to the inverse
Bohr radius, the NLL correction vanishes. Furthermore, at $\nu''=
0.92$~GeV, the result becomes independent of $\nu$; {\it i.e.}, the NLL curve
shows a local maximum corresponding to $E_{\rm hfs}=65$~MeV, which we take as
the central value of our estimate. The NLL curve also shows an impressive
stability with respect to the hard matching scale variation in the physical
range $m_c<\nu_h