# Double charmonium production at B-factories within light cone formalism.

V.V. Braguta<EMAIL_ADDRESS>Institute for High Energy Physics, Protvino,
Russia

###### Abstract

This paper is devoted to the study of the processes $e^{+}e^{-}\to
J/\Psi\eta_{c},J/\Psi\eta_{c}^{\prime},\psi^{\prime}\eta_{c},\psi^{\prime}\eta_{c}^{\prime}$
within light cone formalism. It is shown that if one disregards the
contribution of higher fock states, the twist-3 distribution amplitudes needed
in the calculation can be unambiguously determined from the twist-2
distribution amplitudes and equations of motion. Using models of the twist-2
distribution amplitudes the cross sections of the processes under study have
been calculated. The results of the calculation are in agreement with Belle
and BaBar experiments. It is also shown that relativistic and radiative
corrections to the cross sections play crucial role in the achievement of the
agreement between the theory and experiments. The comparison of the results of
this paper with the results obtained in other papers has been carried out. In
particular, it is shown that the results of papers where relativistic and
radiative corrections were calculated within NRQCD are overestimated by a
factor of $\sim 1.5$.

###### pacs:

12.38.-t, 12.38.Bx, 13.66.Bc,

##  I Introduction

Double charmonium production at B-factories is very interesting problem from
theoretical point of view. At the beginning, the interest to this problem was
caused by large discrepancy between theoretical predictions for the cross
section of the process $e^{+}e^{-}\to J/\Psi\eta_{c}$ Braaten:2002fi ; Liu:1 ;
Liu:2  and it’s first experimental measurement by Belle collaboration
Abe:2002rb . Lately, large discrepancy between theory and experiment was found
by Belle Abe:2004ww  and BaBar Aubert:2005tj  collaborations for the processes
$e^{+}e^{-}\to
J/\Psi\eta_{c}^{\prime},\psi^{\prime}\eta_{c},\psi^{\prime}\eta_{c}^{\prime},J/\Psi\chi_{c0},\psi^{\prime}\chi_{c0}$.

A number of attempts were made to explain this discrepancy. For instance, the
authors of papers Bodwin:2002fk ; Bodwin:2002kk ; Luchinsky:2003yh  studied
the possibility of admixture of $e^{+}e^{-}\to J/\psi J/\psi$ events in the
$e^{+}e^{-}\to J/\Psi\eta_{c}$. Another attempt was to attribute the
discrepancy to large radiative corrections Zhang:2005ch ; Gong:2007db ;
Zhang:2008gp . One more very popular approach to the problem Ma:2004qf ;
Bondar:2004sv ; Braguta:2005gw ; Braguta:2005kr ; Braguta:2006nf ;
Ebert:2006xq ; Choi:2007ze ; Berezhnoy:2007sp ; Ebert:2008kj  consisted in
taking into account internal motion of the quark-antiquark pair in charmonium.
Some other interesting points of view on this problem were proposed in papers
Lee:2003db ; Zhang:2008ab ; Berezhnoy:2006mz .

Among many approaches to the resolution of this challenging problem one should
mention nonrelativistic QCD (NRQCD) Bodwin:1994jh . Contrary to many other
models of quarkonium production this approach allows one to improve
systematically the accuracy of the calculation. Thus at the leading order
approximation Braaten:2002fi ; Liu:1 ; Liu:2  there is very large discrepancy
between NRQCD predictions for the cross section of the process $e^{+}e^{-}\to
J/\Psi\eta_{c}$ and the experiments. However, after the inclusion of the one
loop radiative corrections Zhang:2005ch ; Gong:2007db  this discrepancy
becomes smaller. At the next step, after taking into account radiative and
relativistic corrections simultaneously He:2007te ; Bodwin:2007ga , the
problem can be resolved. From this example one can draw a conclusion, that
relativistic and radiative corrections play very important role in the
processes with charmonia production. Probably, in a similar way one can
resolve many other puzzles with quarkonia production Bodwin .

It should be noted that within NRQCD the discrepancy between the theory and
experiments was resolved only for the process $e^{+}e^{-}\to J/\Psi\eta_{c}$.
At the same time the discrepancies were observed also for the processes
$e^{+}e^{-}\to
J/\Psi\eta_{c}^{\prime},\psi^{\prime}\eta_{c},\psi^{\prime}\eta_{c}^{\prime},J/\Psi\chi_{c0},\psi^{\prime}\chi_{c0}$.
It is not clear is it possible to apply the same approach for these processes,
since excited charmonia states are rather relativistic and the application of
NRQCD to the production of such states is questionable.

Another systematic approach to the study of hard exclusive processes is light
cone formalism (LC) Lepage:1980fj ; Chernyak:1983ej . Within this approach the
amplitude of hard exclusive process can be separated into two parts. The first
part is partons production at very small distances, which can be treated
within perturbative QCD. The second part is the hadronization of the partons
at larger distances. This part contains information about nonperturbative
dynamic of the strong interactions. For hard exclusive processes it can be
parameterized by process independent distribution amplitudes (DA), which can
be considered as hadrons’ wave functions at light like separation between the
partons in the hadron. It should be noted that within LC one does not assume
that the mesons are nonrelativistic. This approach can equally well be applied
to the production of light and heavy mesons, if the DAs of the produced meson
are known. For this reason, one can hope that within this approach one can
study the production of excited charmonia states.

The first attempts to describe the experimental results obtained at Belle and
BaBar collaborations within LC were done in papers Ma:2004qf ; Bondar:2004sv ;
Braguta:2005kr ; Braguta:2006nf . There are two very important problems common
for all these papers. The first one is that in these papers the
renormalization group evolution of the DAs was disregarded. What is very
important since the evolution of the DAs takes into account very important
part of radiative corrections – the leading logarithmic radiative corrections.
The second problem is poor knowledge of charmonia DAs, which are the key
ingredient of any calculation done within LC.

Recently, the leading twist DAs of the $S$-wave charmonia mesons have become
the objects of intensive study Choi:2007ze ; Bodwin:2006dm ; Ma:2006hc ;
Braguta:2006wr ; Braguta:2007fh ; Braguta:2007tq ; Feldmann:2007id ;
Bell:2008er ; Hwang:2008qi . Knowledge about these DAs allowed one to build
some models for the $S$-wave charmonia DAs, that can be used in practical
calculations. In this paper LC will be applied to the study the processes
$e^{+}e^{-}\to
J/\Psi\eta_{c},J/\Psi\eta_{c}^{\prime},\psi^{\prime}\eta_{c},\psi^{\prime}\eta_{c}^{\prime}$.
It will be shown that with the models of DAs proposed in papers Braguta:2006wr
; Braguta:2007fh ; Braguta:2007tq  LC predictions are in agreement with the
results obtained at Belle and BaBar experiments.

This paper is organized as follows. In the next sections a brief description
of LC is given. In section III the formula for the amplitude of the process
$e^{+}e^{-}\to VP$, where $V$ and $P$ are vector and pseudoscalar mesons, is
derived. In section IV this formula is used to calculate the cross sections of
the processes $e^{+}e^{-}\to
J/\Psi\eta_{c},J/\Psi\eta_{c}^{\prime},\psi^{\prime}\eta_{c},\psi^{\prime}\eta_{c}^{\prime}$.
The results of the calculation are discussed in section V. Finally, in the
last section the results of this paper are summarized.

##  II Brief description of light cone formalism.

In this section a brief description of light cone formalism (LC) will be
given. As an example, let us consider hard exclusive process with single meson
production. And let us assume that this meson is a pseudoscalar meson $P$. The
presence of high energy scale $E_{h}$, which is of order of the characteristic
energy of the hard exclusive process, allows one to apply factorization
theorem for the amplitude of the process $T$

$\displaystyle T=\sum_{n}C_{n}\times\langle P|O_{n}|0\rangle,$ (1)

where the coefficient $C_{n}$ describes partons production at small distances,
the matrix element $\langle P|O_{n}|0\rangle$ describes hadronization of the
partons which takes place at large distances. The sum is taken over all
possible operators $O_{n}$. For instance, the operators
$\bar{Q}\gamma_{\mu}\gamma_{5}Q,\bar{Q}\gamma_{\mu}\gamma_{5}{\overset{\leftrightarrow}{D}_{\mu_{1}}}{\overset{\leftrightarrow}{D}}_{\mu_{2}}Q,\bar{Q}\sigma_{\mu\nu}\gamma_{5}G_{\alpha\beta}Q$
are few examples of the operator $O_{n}$. Actually, there are infinite number
of the operators $O_{n}$ that contribute to the pseudoscalar meson production.

The cross section of hard exclusive process can be expanded in inverse powers
of the high energy scale $E_{h}$

$\displaystyle\sigma=\frac{\sigma_{0}}{E_{h}^{n}}+\frac{\sigma_{1}}{E_{h}^{n+1}}+...$
(2)

To determine if some operator contributes to a given term in $1/E_{h}$
expansion one uses the concept of the twist of this operator Braun:2003rp .
Thus only the leading twist – the twist-2 operators
$\bar{Q}\hat{z}\gamma_{5}Q,\bar{Q}\hat{z}\gamma_{5}({z\overset{\leftrightarrow}{D}})Q,\bar{Q}\hat{z}\gamma_{5}({z\overset{\leftrightarrow}{D}})^{2}Q,...$111$z$
here is lightlike fourvector $z^{2}=0$. contribute to the leading term in
expansion (2). From this one sees that already at the leading order
approximation infinite number of operators contribute to the cross section.
Nevertheless, it is possible to cope with infinite number of contributions if
one parameterizes all the twist-2 operators by the moments of some function
$\phi(x)$

$\displaystyle\langle
P(q)|\bar{Q}\hat{z}\gamma_{5}({-i\overset{\leftrightarrow}{D}}_{\mu_{1}})...({-i\overset{\leftrightarrow}{D}}_{\mu_{n}})Q|0\rangle\times
z^{\mu_{1}}...z^{\mu_{n}}=if_{P}(qz)^{n+1}\int_{0}^{1}dx\phi(x)(2x-1)^{n},$
(3)

where $q$ is the momentum of the pseudoscalar meson $P$, $f_{P}$ is the
constant which is defined as $\langle
P(q)|\bar{Q}\gamma_{\mu}\gamma_{5}Q|0\rangle=if_{P}q_{\mu}$, $x$ is the
fraction of momentum of the whole meson $P$ carried by quark. The function
$\phi(x)$ is called the leading twist distribution amplitude (DA). One can
think of it as about the Fourier transform of the wave function of the meson
$P$ with lightlike distance between quarks. Using the definition of this
function, factorization theorem (1) can be rewritten as

$\displaystyle T=\int_{0}^{1}dxH(x)\phi(x),$ (4)

where $H(x)$ is the hard part of the amplitude, which describes small distance
effects. This part of the amplitude can be calculated within perturbative QCD.
As it was noted, the leading twist DA parameterizes infinite set of the
twist-2 operators. This part of the amplitude describes hadronization of
quark-antiquark pair at large distances and parameterizes nonperturbative
effects in the amplitude. It should be noted that formula (4) resums the
contributions of all the twist-2 operators. If the meson $P$ is a
nonrelativistic meson, formula (4) resums relativistic corrections to the
amplitude $T$.

Now let us consider radiative corrections to formula (4). The presence of two
different energy scales, which are strongly separated $E_{h}\gg M_{P}$, gives
rise to the appearance of large logarithm $\log{E_{h}^{2}/M_{P}^{2}}$. This
logarithm enhances the role of radiative corrections. The main contribution to
amplitude (4) comes from the leading logarithmic radiative corrections
$\sim(\alpha_{s}\log{E_{h}^{2}/M_{P}^{2}})^{n}$. It turns out that these
corrections can be taken into the account in formula (4) as follows
Chernyak:1983ej ; Lepage:1980fj

$\displaystyle T=\int_{0}^{1}dxH(x,\mu)\phi(x,\mu).$ (5)

To resum the leading logarithmic corrections coming from all loops the scale
$\mu$ should be taken of order of $\sim E_{h}$. The hard part of the amplitude
$H(x,\mu)$ should be calculated at the tree level approximation. At this level
$H(x,\mu)$ depends on the renormalization scale $\mu$ only through the running
of the strong coupling constant $\alpha_{s}(\mu)$. The rest of the leading
logarithms are resummed in the DA $\phi(x,\mu)$ using renormalization group
method (see Appendix B). It should be stressed that formula (5) exactly resums
the leading logarithmic radiative corrections.

Commonly, to study the production of nonrelativistic mesons one uses effective
theory NRQCD Bodwin:1994jh . NRQCD deals with three energy scales $m_{Q}\gg
m_{Q}v\gg m_{Q}v^{2}$, where $m_{Q}$ is the mass of the heavy quark $Q$, $v\ll
1$ is relative velocity of quark antiquark pair. In the process of hard
nonrelativistic meson production there appears one additional energy scale
$E_{h}$ which is much greater than all scales $m_{Q},m_{Q}v,m_{Q}v^{2}$.
Evidently, it is not possible to apply NRQCD at this scale. From the effective
theory perspective, first, this large energy scale must be integrated out. And
this is done through the taking into account renormalization group evolution
of the DA $\phi(x,\mu)$.

This paper is devoted to the study of the process $e^{+}e^{-}\to VP$, where
$V$ and $P$ are vector and pseudoscalar mesons. Essential feature of this
process is that for it the leading order contribution in $1/E_{h}$ expansion
is zero Bondar:2004sv . So, we deal with the next-to-leading twist process.
For this process formula (4) remains valid. The only difference is that now we
have contributions coming from different twist-2 and twist-3 DAs (see next
section). As in the case of the leading twist process, formula (4) resums
relativistic corrections to the amplitude. However, if we consider the leading
logarithmic radiative corrections to the amplitude, formula (5) is incorrect.
One can expect that some leading logarithms are lost. It is only possible to
state that formula (5) resums the leading logarithms which appear in the
amplitude due to the running of the $\alpha_{s}$ and DAs. Below this
approximation will be used.

There is one common feature of all next-to-leading twist process. It is
connected to the following fact: along with two particles twist-3 operators
(for instance $\bar{Q}\gamma_{5}Q$) there appears operators of the type
$\bar{Q}\gamma_{\mu}\gamma_{5}G_{\alpha\beta}Q$. Evidently, these operators
describe higher fock state $|\bar{Q}Qg>$ contribution to the amplitude of the
process. NRQCD predicts that for nonrelativistic mesons such states are
suppressed by higher powers of relative velocity of quark-aniquark pair inside
the meson Bodwin:1994jh . For this reason, in this paper the contribution of
such states will be disregarded.

##  III The amplitude of the process: $e^{+}e^{-}\to VP$.

In this section the amplitude of the process $e^{+}e^{-}\to VP$, where
$V=J/\Psi,\psi^{\prime}$ and $P=\eta_{c},\eta_{c}^{\prime}$ will be
considered. Two diagrams that give contribution to the amplitude of this
process are shown in Fig 1. The other two can be obtained from the depicted
ones by the charge conjugation. The amplitude of the process involved can be
written in the following form:

$\displaystyle
M=-4\pi\alpha\frac{\bar{u}(k_{1})\gamma_{\mu}u(k_{2})}{s}\langle
V(p_{1},\lambda)P(p_{2})|J_{\mu}^{em}|0\rangle,$ (6)

where $\alpha$ is the electromagnetic coupling constant,
$\bar{u}(k_{1}),u(k_{2})$ are the electron and positron bispinors, $\sqrt{s}$
is the invariant mass of $e^{+}e^{-}$ system, $J_{\mu}^{em}$ is the
electromagnetic current. The matrix element $\langle
V(p_{1},\lambda)P(p_{2})|J_{\mu}^{em}|0\rangle$ can be parameterized by the
only formfactor $F(s)$:

$\displaystyle\langle
V(p_{1},\lambda)P(p_{2})|J_{\mu}^{em}|0\rangle=iq_{c}~{}F(s)~{}e_{\mu\nu\rho\sigma}\epsilon_{\lambda}^{\nu}p_{1}^{\rho}p_{2}^{\sigma},$
(7)

where $q_{c}$ is the charge of $c$ quark, $\epsilon_{\lambda}^{\nu}$ is the
polarization vector of the meson $V(p_{1},\lambda)$. The cross section of the
process under consideration can be written as follows

$\displaystyle\sigma(e^{+}e^{-}\to\ VP)$ $\displaystyle=$
$\displaystyle\frac{\pi\alpha^{2}q_{c}^{2}}{6}\biggl{(}\frac{2|\mathbf{p}|}{\sqrt{s}}\biggr{)}^{3}|F(s)|^{2}.$
(8)

In the last formula $\mathbf{p}$ is the momentum of the meson $V$ in the
center mass frame of the final mesons.

To calculate the formfactor $F(s)$ LC will be applied. As it was noted in the
previous section, within LC the formfactor is a series in inverse powers of
the characteristic energy of the process $\sqrt{s}$. The leading order
contribution to the formfactor $F(s)$ in the $1/s$ expansion was obtained in
paper Bondar:2004sv . In derivation of this expression the authors of paper
Bondar:2004sv  disregarded the mass difference of the final mesons. For the
mesons with different masses the expression for the formfactor $F(s)$, which
was derived in paper Braguta:2005kr , can be written as follows

$\displaystyle|F(s)|$ $\displaystyle=$
$\displaystyle\frac{32\pi}{9}\left|\frac{f_{V}f_{P}M_{P}M_{V}}{q_{0}^{4}}\right|\,I_{0}\,,$
(9) $\displaystyle I_{0}$ $\displaystyle=$
$\displaystyle\int^{1}_{0}dx\int^{1}_{0}dy~{}\alpha_{s}(\mu)\left\\{\frac{M_{P}}{M_{V}^{2}}\frac{Z_{t}Z_{p}V_{T}(x,\mu)P_{P}(y,\mu)}{d(x,y)\,s(x)}-\frac{1}{M_{P}}\frac{\overline{M}_{c}^{2}}{{M_{V}}^{2}}\,\frac{Z_{m}Z_{t}V_{T}(x,\mu)P_{A}(y,\mu)}{d(x,y)\,s(x)}+\right.$
$\displaystyle+$
$\displaystyle\frac{1}{2M_{P}}\frac{V_{L}(x,\mu)\,P_{A}(y,\mu)}{d(x,y)}+\frac{1}{2M_{P}}\frac{(1-2y)}{s(y)}\frac{V_{\perp}(x,\mu)\,P_{A}(y,\mu)}{d(x,y)}+$
$\displaystyle+$
$\displaystyle\left.\frac{1}{8}\biggl{(}1-Z_{t}Z_{m}\frac{4{\overline{M}}_{c}^{2}}{{M_{V}}^{2}}\biggr{)}\frac{1}{M}_{P}\,\frac{(1+y)V_{A}(x,\mu)P_{A}(y,\mu)}{d^{2}(x,y)}\right\\}.$
(11)

Where $q_{0}^{2}\simeq(s-M_{V}^{2}-M_{P}^{2})$, $x$ and $y$ are the fractions
of momenta carried by quark in the meson V and by quark in the meson P
correspondingly, $P_{A},P_{P},V_{T},V_{L},V_{\perp},V_{A}$ are the DAs defined
in Appendix B, $M_{V},M_{P}$ are the masses of the vector and pseudoscalar
mesons correspondingly,
$\overline{M}_{c}=M^{\overline{MS}}_{c}(\mu=M^{\overline{MS}}_{c})$, $Z_{t}$
and $Z_{p}$ are the renormalization factors of the local tensor and
pseudoscalar currents, the dimensionless propagators $d(x,y),s(x),s(y)$ are
defined as follows:

$\displaystyle d(x,y)$ $\displaystyle=$
$\displaystyle\frac{k^{2}}{q_{0}^{2}}=\left(x+\frac{\delta}{y}\right)\left(y+\frac{\delta}{x}\right),~{}s(x)=\frac{k_{x}^{2}}{q_{0}^{2}}-\delta=\left(x+\frac{\delta}{y(1-y)}\right),~{}s(y)=\frac{k_{y}^{2}}{q_{0}^{2}}-\delta=\left(y+\frac{\delta}{x(1-x)}\right),$
(12) $\displaystyle Z_{p}$ $\displaystyle=$
$\displaystyle\left[\frac{\alpha_{s}(\mu)}{\alpha_{s}({\overline{M}}_{c}^{2})}\right]^{\frac{-3c_{F}}{b_{o}}},\quad
Z_{t}=\left[\frac{\alpha_{s}(\mu)}{\alpha_{s}({\overline{M}}_{c}^{2})}\right]^{\frac{c_{F}}{b_{o}}},\quad
Z_{m}=\left[\frac{\alpha_{s}(\mu)}{\alpha_{s}({\overline{M}}_{c}^{2})}\right]^{\frac{3c_{F}}{b_{o}}},\quad\delta=\Biggl{(}Z_{m}\frac{{\overline{M}}_{c}}{q_{0}}\Biggr{)}^{2}\,,$

where $c_{F}=4/3,\,b_{o}=25/3$, the definitions of the fourvectors $k,k_{x}$
and $k_{y}$ are shown in Fig. 1.

Figure 1: The diagrams that contribute to the process $e^{+}e^{-}\to
V(p_{1},\lambda)P(p_{2})$ at the leading order approximation in the strong
coupling constant.

It should be noted here that if one ignores renormalization group running,
formula (9) is the analog of formula (4). So, in this case, it takes into
account infinite series of the relativistic corrections to double charmonium
production. If one takes into account the leading logarithmic corrections due
to the DAs and the running of the strong coupling constant, the scale $\mu$
should be taken of order of the characteristic energy of the process
$\sim\sqrt{s}$.

Before we proceed with the calculation of double charmonium production, the
expression for the formfactor $F(s)$ will be modified as follows:

1\. Formula (9) was obtained at the first nonvanishing approximation in the
$1/s$ expansion. This implies that all parameters in this formula must be
taken at the same level of accuracy. For instance, in the expression
$q_{0}^{2}\simeq(s-M_{V}^{2}-M_{P}^{2})$ last two terms are beyond the
accuracy of calculation and these terms must be omitted. Below the following
approximation will be used $q_{0}^{2}=s$.

2\. Now let us consider expressions (12) for the propagators
$d(x,y),s(x),s(y)$. It is seen that all these expressions contain the terms
proportional to $\sim\delta$. These terms are of NLO approximation, so,
according to the previous item, they must be omitted. Thus the propagators can
be written as

$\displaystyle d(x,y)=xy,\quad s(x)=x,\quad s(y)=y.$ (13)

However, if one substitutes these expressions to (III) and takes into account
the renormalization group evolution of the DAs, the divergence at the end
point region ($x,y\sim 0$) appears. It should be noted that the motion of
quark-antiquark pair in the end point regions $x,y\sim 0,1$ is relativistic.
So, if we consider the production of a nonrelativistic meson without the
evolution of the DAs, the end point regions in the DA are strongly suppressed.
For this reason the expression for the formfactor $F(s)$ is free from the
divergence. However, if the evolution of the DAs is taken into account, the
divergence in the $F(s)$ appears. In this paper this problem will be solved as
follows. According to the definition of $x$

$\displaystyle x=\frac{E^{c}+p^{c}_{z}}{E^{M}+p^{M}_{z}},$ (14)

where $E^{c},p^{c}_{z}$ are the energy and $z$ component of the momenum of
$c$-quark, $E^{M},p^{M}_{z}$ are the energy and $z$ component of the momenum
of meson $M$. If the energy and momenum of $c$-quark is much greater than it’s
mass, one can use the propagators in form (13), since corrections to this
approximation are suppressed. Let us now consider the kinematic region where
$c$-quark is approximately at rest and $x$ can be estimated as

$\displaystyle
x\sim\frac{M^{*}_{c}}{E^{M}+p^{M}_{z}}\simeq\frac{M_{c}^{*}}{\sqrt{s}}=x_{min}.$
(15)

$M_{c}^{*}$ here is the pole mass of $c$ quark. In this region the accuracy of
propagators (13) is not sufficient and one must take into account the
corrections. It is clear that in any approach the corrections to the
propagators regularize the whole expression for the formfactor $F(s)$. This
fact can be seen as follows: for any kinematical region of double quark and
double antiquark production shown in Fig. 1 the squares of momenta of the
gluon and quark propagators cannot be smaller than $(2M_{c}^{*})^{2}$ and
$(M_{c}^{*}+M_{V,P})^{2}$ correspondingly. So, there is no divergence in the
exact expressions for the propagators.

In this paper this effect will be taken into account as follows: the
propagators will be taken in form (13), but the integration in the expression
for the formfactor $F(s)$ will be done in the region $x,y\in(x_{min},~{}1)$.
In other words to get rid of the singularity the cut off parameters $x_{min}$
is introduced.

In principle, the calculation with propagators (12) is also possible. As it
was noted already, in this case the terms proportional to the $\delta$ play
role of the regulator of expression (III). However, one can expect that the
calculation done in this manner is less accurate and the result must be
smaller than it is. To understand this one can apply the idea of duality of
NRQCD and LC descriptions of the hard exclusive nonrelativistic mesons
production: if LC expression for the formfactor $F(s)$ is expanded in relative
velocities of quark antiquark pairs of the mesons $V$ and $P$, one will get
NRQCD result for the formfactor. 222 It should be noted that these is no
strict proof of this statement. However, one can expect that this statement is
indeed true since the amplitude in NRQCD and LC can be expended in series of
equivalent operators. Assuming that both theories can describe experiment one
can expect that these expansions in both theories coincide.

At the leading order approximation of NRQCD this idea can be reproduced if one
takes infinitely narrow approximation for the DAs
($P_{i}(x)=V_{j}(x)\sim\delta(x-1/2),~{}i=A,P,~{}~{}j=T,L,{\perp},A$ ) and
uses the following values of the parameters in formula (9)
$M_{V}=M_{P}=2M_{c}^{*},{\overline{M}}_{c}=M_{c}^{*},\mu=M_{c}^{*},f_{V,P,T}^{2}=\langle
O_{1}\rangle_{S}/M_{c}^{*}$. If this procedure is applied to (9) with
propagators (13), we will exactly reproduce the expression for the formfactor
$F(s)$ obtained within the leading order of NRQCD Braaten:2002fi

$\displaystyle F_{NRQCD}=\frac{2^{9}\pi\alpha_{s}}{9s^{2}}\langle
O_{1}\rangle_{S},$ (16)

where $\langle O_{1}\rangle_{S}$ is the NRQCD matrix element Braaten:2002fi .
Note that at the leading order approximation of NRQCD the formfactor $F(s)$
scales exactly as $1/s^{2}$. So, within NRQCD $1/s^{3}$ terms appear due to
the relativistic corrections, which are suppressed as $v^{2}$ ($v^{2}$ is the
characteristic velocity in charmonium). In LC $1/s^{3}$ terms can appear only
due to the power corrections of the leading order result. So, applying the
idea of duality of NRQCD and LC one can state that the power $1/s$ corrections
to formula (9) with propagators (13) are of order of $\sim
v^{2}(M_{c}^{*})^{2}/s$. If we further apply the same procedure but with
propagators (12), the formfactor $F(s)$ will be different from the leading
order NRQCD prediction

$\displaystyle
F(s)=F_{NRQCD}(1-14\delta+O(\delta^{2}))=F_{NRQCD}\biggl{(}1-14\frac{(M_{c}^{*})^{2}}{s}+O\biggl{(}\frac{1}{s^{2}}\biggr{)}\biggr{)}.$
(17)

To get the agreement with NRQCD prediction for the formfactor, one must expect
rather large power correction to this result ($14\delta\times F_{NRQCD}$) in
LC so that to cancel the second term in the last equation. This correction
appears at next-to-leading order approximation in $1/s$ expansion. So, the
leading order approximation of LC with propagators (12) underestimates real
result. It is not difficult to estimate the size of this effect using formula
(17). For $\sqrt{s}=10.6$ GeV and $m_{c}^{*}=1.4$ GeV the cross section
calculated with propagators (12) is smaller than that with propagators (13) by
$\sim 50\%$. Numerical calculation confirms this estimation.

3\. To calculate the cross section of the processes $e^{+}e^{-}\to
J/\Psi\eta_{c},J/\Psi\eta_{c}^{\prime},\psi^{\prime}\eta_{c},\psi^{\prime}\eta_{c}^{\prime}$
one needs to know the following constants:

$\displaystyle\langle J/\Psi(p,\epsilon)|\bar{C}\gamma_{\alpha}C|0\rangle$
$\displaystyle=$ $\displaystyle
f_{V1}M_{J/\Psi}\epsilon_{\alpha}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\langle\psi^{\prime}(p,\epsilon)|\bar{C}\gamma_{\alpha}C|0\rangle=f_{V2}M_{\psi^{\prime}}\epsilon_{\alpha},$
(18) $\displaystyle\langle
J/\Psi(p,\epsilon)|\bar{C}\sigma_{\alpha\beta}C|0\rangle_{\mu}$
$\displaystyle=$ $\displaystyle
if_{T1}(\mu)(p_{\alpha}\epsilon_{\beta}-p_{\beta}\epsilon_{\alpha})~{}~{}~{}~{}\langle\psi^{\prime}(p,\epsilon)|\bar{C}\sigma_{\alpha\beta}C|0\rangle_{\mu}=if_{T2}(\mu)(p_{\alpha}\epsilon_{\beta}-p_{\beta}\epsilon_{\alpha}),$
$\displaystyle\langle\eta_{c}(p,\epsilon)|\bar{C}\gamma_{\alpha}\gamma_{5}C|0\rangle$
$\displaystyle=$ $\displaystyle
if_{P1}p_{\alpha}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\langle\eta_{c}^{\prime}(p,\epsilon)|\bar{C}\gamma_{\alpha}\gamma_{5}C|0\rangle=if_{P2}p_{\alpha}.$

It should be noted that the operator ${\bar{C}}\sigma_{\alpha\beta}C$ is not
renormalization group invariant. For this reason the constants $f_{T1}$ and
$f_{T2}$ depend on scale as

$\displaystyle
f_{T}(\mu)=\biggl{(}\frac{\alpha_{s}(\mu)}{\alpha_{s}(\mu_{0})}\biggr{)}^{\frac{c_{F}}{b_{0}}}f_{T}(\mu_{0}),~{}~{}~{}i=1,2.$
(19)

In the derivation of formula (9) it was assumed that the tensor and vector
constants $f_{Vi},f_{Ti}(\overline{M}_{c})$ are connected to each other as
$f_{Ti}(\overline{M}_{c})/f_{Vi}=2\,{M}_{c}/M_{Vi},~{}i=1,2$ Bondar:2004sv ;
Braguta:2005kr . At the leading order approximation in relative velocity and
strong coupling constant these relations are correct. However, they are
violated due to radiative and relativistic corrections especially for the
excited mesons. For this reason, below $f_{V},f_{T}$ will be treated as
independent constants.

Introducing the modifications described above, the expression for the
formfactor $F(s)$ can be rewritten as follows

$\displaystyle|F(s)|$ $\displaystyle=$
$\displaystyle\frac{32\pi}{9}\left|\frac{f_{V}f_{P}M_{P}M_{V}}{s^{2}}\right|\,I_{0}\,,$
(20) $\displaystyle
I_{0}=\int^{1}_{x_{min}}dx\int^{1}_{y_{min}}dy~{}\alpha_{s}(\mu)\left\\{\frac{\tilde{f}_{t}}{M_{V}}\frac{M_{P}}{2{\overline{M}}_{c}}\frac{Z_{p}V_{T}(x,\mu)P_{P}(y,\mu)}{x^{2}\,y}-\frac{\tilde{f}_{t}}{M_{P}}\frac{\overline{M}_{c}}{2M_{V}}\,\frac{Z_{m}V_{T}(x,\mu)P_{A}(y,\mu)}{x^{2}\,y}+\right.$
$\displaystyle+$
$\displaystyle\frac{1}{2M_{P}}\frac{V_{L}(x,\mu)\,P_{A}(y,\mu)}{x\,y}+\frac{(1-2y)}{2M_{P}}\frac{V_{\perp}(x,\mu)\,P_{A}(y,\mu)}{x\,y^{2}}+\left.\frac{1}{8}\biggl{(}1-{\tilde{f}_{t}}Z_{m}\frac{2{\overline{M}}_{c}}{{M_{V}}}\biggr{)}\frac{1}{M}_{P}\,\frac{(1+y)V_{A}(x,\mu)P_{A}(y,\mu)}{x^{2}\,y^{2}}\right\\},$

where $\tilde{f}_{t}=f_{T}(\mu)/f_{V}$. Now we are ready to proceed with the
calculation.

##  IV Numerical results.

To calculate the cross sections of the processes $e^{+}e^{-}\to
J/\Psi\eta_{c},J/\Psi\eta_{c}^{\prime},\psi^{\prime}\eta_{c},\psi^{\prime}\eta_{c}^{\prime}$
the following values of input parameters will be used:

1. The strong coupling constant $\alpha_{s}(\mu)$ will be taken at the one loop approximation

$\displaystyle\alpha_{s}(\mu)=\frac{4\pi}{\beta_{0}\log(\mu^{2}/\Lambda^{2})},$
(22)

with $\Lambda=0.2$ GeV, $\beta_{0}=25/3$.

2. For the $\overline{MS}$ mass and the pole mass of $c$-quark the values ${\overline{M}}_{c}=1.2$ GeV and $M_{c}^{*}=1.4$ GeV will be used.

3. As it was noted in the previous section, to calculate the cross sections one needs constants (18). The constants $f_{V1}$ and $f_{V2}$ can be determined directly from the experiment. The constants $f_{T1}$ and $f_{T2}$ were calculated within NRQCD in paper Braguta:2007ge . NRQCD formalism can also be used to determine the values of the last two constants $f_{P1}$ and $f_{P2}$ (see Appendix A). So, the calculation will be done with the following values of these constants

$\displaystyle f_{V1}^{2}$ $\displaystyle=$ $\displaystyle 0.173\pm
0.004~{}\mbox{GeV}^{2},\qquad~{}~{}~{}~{}~{}~{}f_{V2}^{2}=0.092\pm
0.002~{}\mbox{GeV}^{2},$ (23) $\displaystyle f_{T1}^{2}(M_{J/\Psi})$
$\displaystyle=$ $\displaystyle 0.144\pm 0.016~{}\mbox{GeV}^{2},\quad
f_{T2}^{2}(M_{J/\Psi})=0.068\pm 0.022~{}\mbox{GeV}^{2},$ $\displaystyle
f_{P1}^{2}$ $\displaystyle=$ $\displaystyle 0.139\pm
0.048~{}\mbox{GeV}^{2},\quad~{}~{}~{}~{}~{}~{}~{}~{}~{}f_{P2}^{2}=0.068\pm
0.040~{}\mbox{GeV}^{2}.$

4. The last inputs to formula (20) are the DAs. The models of the leading twist DAs will be taken from papers Braguta:2006wr ; Braguta:2007fh ; Braguta:2007tq . To build the models for the twist-3 DAs one can apply equations of motion. This procedure is described in detail in Appendix B.

$\begin{array}[]{|c|c|c|c|c|c|c|c|c|}\hline\cr
H_{1}H_{2}&\lx@intercol\hfil\sigma_{Exp}\times Br_{H_{2}\to
charged>2}(\mbox{fb})\hfil\lx@intercol\vrule\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\sigma_{LO~{}NRQCD}(\mbox{fb})\hfil\lx@intercol\vrule\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\sigma_{NRQCD}(\mbox{fb})\hfil\lx@intercol\vrule\lx@intercol\vrule\lx@intercol&\sigma_{poten~{}model}(\mbox{fb})&\sigma_{LC}(\mbox{fb})\\\
\hline\cr&\mbox{Belle}\mbox{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Abe:2004ww}{\@@citephrase{(}}{\@@citephrase{)}}}}&\mbox{BaBar}\mbox{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Aubert:2005tj}{\@@citephrase{(}}{\@@citephrase{)}}}}&\mbox{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Braaten:2002fi}{\@@citephrase{(}}{\@@citephrase{)}}}}&~{}~{}\mbox{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Liu:1}{\@@citephrase{(}}{\@@citephrase{)}}}}{}{}&\mbox{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{He:2007te}{\@@citephrase{(}}{\@@citephrase{)}}}}&\mbox{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Bodwin:2007ga}{\@@citephrase{(}}{\@@citephrase{)}}}}&\mbox{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Ebert:2008kj}{\@@citephrase{(}}{\@@citephrase{)}}}}&\\\
\hline\cr\psi(1S)\eta_{c}(1S)&25.6\pm 2.8\pm 3.4&17.6\pm
2.8^{+1.5}_{-2.1}&3.78\pm 1.26&5.5&20.4&17.6^{+10.7}_{-8.3}&22.2\pm
1.1&14.4^{+11.2}_{-9.8}\\\ \hline\cr\psi(2S)\eta_{c}(1S)&16.3\pm 4.6\pm
3.9&-&1.57\pm 0.52&3.7&-&-&15.3\pm 0.8&10.4^{+9.2}_{-7.8}\\\
\hline\cr\psi(1S)\eta_{c}(2S)&16.5\pm 3.0\pm 2.4&16.4\pm
3.7^{+2.4}_{-3.0}&1.57\pm 0.52&3.7&-&-&16.4\pm 0.8&13.0^{+12.2}_{-11.0}\\\
\hline\cr\psi(2S)\eta_{c}(2S)&16.0\pm 5.1\pm 3.8&-&0.65\pm 0.22&2.5&-&-&9.6\pm
0.5&9.0^{+9.7}_{-8.5}\\\ \hline\cr\end{array}$ Table 1:  The second and third
columns contain experimental results measured at Belle and Babar experiments.
The $Br_{H_{2}\to charged>2}$ means the branching ratio of the decay of the
hadron $H_{2}$ into two charged particles. In the fourth and fifth columns the
results of the leading order NRQCD obtained in papers Liu:1 ; Braaten:2002fi
are shown. The NRQCD results obtained with inclusion of radiative and
relativistic corrections He:2007te ; Bodwin:2007ga  are shown in columns six
and seven. Potential model predictions Ebert:2008kj  for the cross sections
are presented in column eight. Last column contains the values of the cross
sections obtained in this paper.

There are different sources of uncertainty to the results obtained in this
paper. The most important uncertainties can be divided into the following
groups:

1. The uncertainty in the models of the distribution amplitudes $\phi_{i}(x,\mu)$, which can be modeled by the variation of the parameters of these models (39). The calculation shows that for the processes $e^{+}e^{-}\to J/\Psi\eta_{c},J/\Psi\eta_{c}^{\prime},\psi^{\prime}\eta_{c},\psi^{\prime}\eta_{c}^{\prime}$. these uncertainties are not greater than $\sim 5\%,~{}12\%,~{}28\%,~{}40\%$ correspondingly. It is seen that these uncertainties are not very large. This fact results from the property found in papers Braguta:2006wr : the evolution improves any model of DA.

2. The uncertainty due to the radiative corrections. In the approach applied in this paper the leading logarithmic radiative corrections due to the evolution of the DAs and strong coupling constant were resummed. As it was noted above for the leading twist processes the resummation of the leading logarithms in the DAs and strong coupling constant is equivalent to the resummation of the leading logarithms in the whole amplitude. Unfortunately, this is not valid for the next-to-leading twist processes, for which some of the leading logarithms are lost. For this reason the radiative corrections should be estimated as $\sim\alpha_{s}(\sqrt{s}/2)\log(s/4/(M_{c}^{*})^{2})\sim 50\%$.

3. The uncertainty due to the power corrections. This uncertainty is determined by the next-to-leading order contribution in the $1/s$ expansion. As it was noted above due to the application of the propagators in form (13) one can hope that these corrections are suppressed by the square of relative velocity of quark antiquark pair in the meson. In the calculation these corrections will be estimated as $\sim 4v^{2}_{\psi^{\prime}}M_{\psi^{\prime}}^{2}/s\sim 25\%$.

4. The uncertainty due to the regularization procedure. As it was noted above, to get rid of the divergence in the formfactor $F(s)$ the cut of parameter was introduced. Evidently, our results depend on the value of the cut of parameter $x_{min}\cdot\sqrt{s}$, which is of order of the mass of $c$-quark (see formula (15) ). To estimate this source of uncertainty the cut off parameter is varied in the region $x_{min}\cdot\sqrt{s}=1.0-1.6$ GeV. The calculation shows that at $x_{min}\cdot\sqrt{s}=1.0$ GeV the cross sections are increased by $\sim 50\%$ and at $x_{min}\cdot\sqrt{s}=1.6$ GeV the cross sections are decreased by $\sim 20\%$. It should be noted that this source of uncertainty is closely connected with the uncertainty due to the radiative corrections. However, to understand this in detail one needs the theory which takes into account all the leading logarithmic corrections.

5. The uncertainty in the values of constants (23). It should be noted that this source of uncertainty is very important especially for the production of the excited states $\psi^{\prime},\eta_{c}^{\prime}$. Thus, for the processes $e^{+}e^{-}\to J/\Psi\eta_{c},\psi^{\prime}\eta_{c},J/\Psi\eta_{c}^{\prime},\psi^{\prime}\eta_{c}^{\prime}$ the errors due to the uncertainties in the values of constants (23) are $\sim 35\%,40\%,60\%,65\%$ correspondingly.

Adding all the uncertainties in quadrature one gets the total errors of the
calculations.

The results of the calculation are presented in Table I. The second and third
columns contain experimental results measured at Belle and Babar experiments.
In the fourth and fifth columns the results of the leading order of NRQCD
approach obtained in papers Liu:1 ; Braaten:2002fi  are shown. The NRQCD
results obtained with inclusion of radiative and relativistic corrections
He:2007te ; Bodwin:2007ga  are shown in columns six and seven. Potential model
predictions Ebert:2008kj  for the cross sections are presented in column
eight. Last column contains the values of the cross sections obtained in this
paper.

##  V Discussion.

It is seen from Tab. I that within the accuracy of the calculation the results
of this paper are in agreement with Belle and BaBar experiments. It is also
seen that the uncertainty of the calculation is rather large. There are two
very important sources of uncertainty. The first one is the theoretical
problem with taking into account of all leading logarithmic radiative
corrections to the amplitude of the next-to-leading twist processes. It can be
estimated as $\sim 50-70\%$ of the cross sections. This source of uncertainty
can be reduced if the theory, which takes into account all leading logarithmic
corrections, is created. The second very important source of uncertainty is
poor knowledge of constants (18). For some processes this uncertainty can
reach $60\%$. The values of these constants used in this paper can be
considered as the first estimation of their real values. So, theoretical and
experimental study of these constants can greatly improve the accuracy of any
predictions done within LC.

$\begin{array}[]{|c|c|c|c|c|}\hline\cr
H_{1}H_{2}&\sigma_{LO~{}NRQCD}(\mbox{fb})&\sigma_{rel~{}corr}(\mbox{fb})&\sigma_{tot}(\mbox{fb})\\\
\hline\cr\psi(1S)\eta_{c}(1S)&1.9&6.3&14.4\\\
\hline\cr\psi(2S)\eta_{c}(1S)&1.0&6.2&13.0\\\
\hline\cr\psi(1S)\eta_{c}(2S)&1.0&7.8&10.4\\\
\hline\cr\psi(2S)\eta_{c}(2S)&0.53&7.2&9.0\\\ \hline\cr\end{array}$ Table 2:
The second column contains the values of the cross sections obtained at the
leading NRQCD approximation. The third column contains the cross sections
calculated in the following approximation: all relativistic corrections are
resummed, but the leading logarithmic radiative corrections are not taken into
account. The last column represents the results obtained if the relativistic
and leading logarithmic radiative corrections to the amplitude are taken into
account simultaneously.

Next, let us discuss the results obtained in other papers and compare them
with the results of this paper. First, let us consider the results of the
leading order NRQCD predictions Braaten:2002fi ; Liu:1  shown in columns four
and five of Tab. I. These results are approximately by an order of magnitude
smaller than the cross sections measured at the experiments. At the same time
LC predictions are in reasonable agreement with the experiments. This facts
lead to the question: why LC predictions are much greater than the leading
order NRQCD predictions? The answer to this question can be given within LC.
As it was noted in section III, LC can reproduce the leading order NRQCD
results. To do this the renormalization group evolution of the constants and
DAs will be disregarded and all parameters will be taken at the leading NRQCD
approximation: the constants $f_{Ti},f_{Pi}$ are equal to the constant
$f_{Vi},i=1,2$, which can be determined from the leptonic decay width, the
mass $\overline{M}_{c}=M_{c}^{*}=1.4$ GeV, $M_{V}=M_{P}=2M_{c}^{*}$, all DAs
are taken at the infinitely narrow approximation $\sim\delta(x-1/2)$. Thus one
gets the results shown in the second column of Tab. II. At the second step all
parameters used in the calculation are taken at their central values, but
without renormalization group evolution. At this step infinite series of the
relativistic corrections are resummed, but the leading logarithmic radiative
corrections are not taken into account (see section II). The results are shown
in the third column of Tab. II. At the last step, renormalization group
evolution is taken into account and the results are presented in the last
column of Tab. II. Within LC this means that the relativistic and leading
logarithmic radiative corrections are taken into account simultaneously.

From Tab. II one sees that the relativistic and leading logarithmic radiative
corrections taken into account simultaneously dramatically enhance the leading
NRQCD predictions and bring the agreement with Belle and BaBar experiments.
Very important conclusion which can be drawn from this result is that in hard
exclusive processes relativistic and leading logarithmic radiative corrections
play very important role and the consideration of such processes at the
leading NRQCD approximation is unreliable.

Looking to the results of Tab. II one can also draw a conclusion that
relativistic corrections alone cannot describe the experimental results. This
conclusion is in agreement with the results of papers Braaten:2002fi ;
Braguta:2005gw ; Ebert:2006xq ; Berezhnoy:2007sp . In these papers the cross
section of the process $e^{+}e^{-}\to J/\Psi\eta_{c}$ is enhanced due to the
relativistic corrections by a factor of $\sim 2-3$ what is in agreement with
the results obtained in this paper. At the same time this conclusion is in
disagreement with the results of paper Ebert:2008kj  (see Tab.I column eight),
where the authors tried to attribute the disagreement between theory and
experiment only to the relativistic corrections, which were calculated within
potential model. From the results shown in Tab.II one sees that this
approximation is realistic only for the production of excited states
$e^{+}e^{-}\to\psi^{\prime}\eta_{c}^{\prime}$. In this case the relativistic
corrections are much more important than the leading logarithmic radiative
corrections.

The authors of papers Bondar:2004sv ; Braguta:2005kr  took into account only
the part of the leading logarithmic radiative corrections which appears due to
the evolution of the constant $f_{T}$ and the running mass of $c$-quark. The
evolution of the DAs was disregarded. The calculation done in this paper shows
that this approximation can be applied only to the DAs of the $2S$ charmonia
mesons, for which the evolution is not very important ( see paper
Braguta:2007tq ). As it was shown in papers Braguta:2006wr ; Braguta:2007fh
the evolution of the $1S$ state charmonia DAs is very important. To compensate
the effect of the evolution of the DAs the authors of paper Bondar:2004sv
proposed rather wide model of the DAs with so called relativistic tail333It
should be noted that relativistic tail of the DAs used in this paper is absent
at the scale $\mu\sim M_{c}^{*}$ and appears due to the evolution of the DAs
at larger scales. (see paper Bodwin:2006dm ). For instance, the characteristic
velocity of the $1S$ charmonia (see formula (32)) with this DA can be
estimated as

$\displaystyle\langle v^{2}\rangle_{1S}\sim 3\langle\xi^{2}\rangle_{1S}=0.39,$
(24)

which is much larger than $\langle v^{2}\rangle_{1S}=0.2-0.3$ calculated in
papers Bodwin:2006dn ; Braguta:2006wr ; Braguta:2007fh ; Choi:2007ze ;
Bodwin:2007fz .

At the end of this section let us consider how the disagreement between theory
and experiment can be resolved within NRQCD. The authors of papers He:2007te ;
Bodwin:2007ga  took into account the relativistic and one loop radiative
corrections and got the following values of the cross section of the process
$e^{+}e^{-}\to J/\Psi\eta_{c}$

$\displaystyle\sigma(e^{+}e^{-}\to J/\Psi\eta_{c})$ $\displaystyle=$
$\displaystyle
20.4~{}\mbox{fb}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mbox{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{He:2007te}{\@@citephrase{(}}{\@@citephrase{)}}}},$
$\displaystyle\sigma(e^{+}e^{-}\to J/\Psi\eta_{c})$ $\displaystyle=$
$\displaystyle
17.6^{+10.7}_{-8.3}~{}\mbox{fb}~{}~{}~{}~{}\mbox{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Bodwin:2007ga}{\@@citephrase{(}}{\@@citephrase{)}}}}.$
(25)

So, similarly to this paper, the authors of papers He:2007te ; Bodwin:2007ga
resolved the disagreement through the taking into account of the relativistic
and radiative corrections.

Central values (25) are larger than the central values of the cross section
obtained in this paper. Below it will be shown that central values (25) are
overestimated. A simple way to understand why the cross section was
overestimated is to consider paper Gong:2007db  which has similar problem. In
this paper explicit expression for the one loop radiative corrections to the
cross section of the process $e^{+}e^{-}\to J/\Psi\eta_{c}$ were calculated.
The result of this paper can be written in the form

$\displaystyle\sigma=\sigma^{0}\times\biggl{(}1+\frac{\alpha_{s}(\mu)}{\pi}K\biggr{)},$
(26)

where $\sigma^{0}$ is the cross section at the leading order approximation,
the expression for the factor $K$ can be found in Gong:2007db . The cross
section $\sigma^{0}$ is proportional to
$|R_{J/\Psi}(0)|^{2}\times|R_{\eta_{c}}(0)|^{2}$, where
$R_{J/\Psi}(r),R_{\eta_{c}}(r)$ are the radial wave functions of the $J/\Psi$
and $\eta_{c}$ mesons. It should be noted here that the cross section $\sigma$
is very sensitive to the values of the wave functions at the origin
$R_{J/\Psi}(0),R_{\eta_{c}}(0)$, so it is very important how these parameters
were calculated. In the calculation the authors of Gong:2007db  took
$|R_{J/\Psi}(0)|=|R_{\eta_{c}}(0)|$ and the value of the $|R_{J/\Psi}(0)|$ was
taken from the leptonic decay width $\Gamma_{ee}$

$\displaystyle|R_{J/\Psi}(0)|^{2}=\frac{1}{1-\frac{16}{3}\frac{\alpha_{s}}{\pi}}\frac{M_{J/\Psi}^{2}\Gamma_{ee}}{4\alpha^{2}q_{c}^{2}}.$
(27)

Now few comments are in order. First, this expression is taken at the next-to-
leading order approximation in the strong coupling constant, what means that
some part of the one loop radiative corrections is put to the $\sigma^{0}$.
Second, if formula (27) is expanded in the $\alpha_{s}$, one will get infinite
series in the strong coupling constant. So, the application of formula (27) in
the calculation of the wave function at the origin is equivalent to the
statement that one knows all infinite series of the radiative corrections to
the wave function. Evidently, this is not correct. To be consistent with the
one loop approximation applied in paper Gong:2007db , one should expand
expression (27) in the strong coupling constant and than leave only the first
term in this expansion

$\displaystyle|R_{J/\Psi}(0)|^{2}=\biggl{(}1+\frac{16}{3}\frac{\alpha_{s}}{\pi}\biggr{)}\frac{M_{J/\Psi}^{2}\Gamma_{ee}}{4\alpha^{2}q_{c}^{2}}.$
(28)

It is not difficult to see that the $|R_{J/\Psi}(0)|^{2}$ calculated from
formula (27) is greater than that calculated form formula (28) by a factor
$1/(1-(16\alpha_{s}/3\pi)^{2})$. The cross section calculated using formula
(27) is greater than the cross section calculated using formula (28) by a
factor $1/(1-(16\alpha_{s}/3\pi)^{2})^{2}\sim 1.5$ for $\alpha_{s}=0.25$. So,
the values of the wave functions at the origin calculated in paper Gong:2007db
were overestimated, what led to the overestimation of the cross section by a
factor of $\sim 1.5$.

Similar overestimation of the wave functions at the origin and, as the result,
overestimation of the cross section takes place in papers He:2007te ;
Bodwin:2007ga . For instance, in paper Bodwin:2007ga  the authors used the
value of the NRQCD matrix element $\langle O_{1}\rangle_{J/\Psi}$, which is
analog of the wave function at the origin, calculated in paper Bodwin:2007fz
using the formula

$\displaystyle\langle O_{1}\rangle_{J/\Psi}=\frac{1}{(1-f(\langle
v^{2}\rangle_{J/\Psi})-\frac{8}{3}\frac{\alpha_{s}}{\pi})^{2}}\frac{3M_{J/\Psi}^{2}\Gamma_{ee}}{8\pi\alpha^{2}q_{c}^{2}},$
(29)

where $\langle v^{2}\rangle_{J/\Psi}$ is defined in (32),
$f(x)=x/(3(1+x+\sqrt{1+x}))$. Similar expression can be written for the
$\langle O_{1}\rangle_{\eta_{c}}$. It is clear that the application of formula
(29) is equivalent to the statement that one knows not only full $\alpha_{s}$
corrections in front of the leading NRQCD operator $\langle
O_{1}\rangle_{J/\Psi}$ but also full $\alpha_{s}$ corrections in front of all
operators that control relativistic corrections $\langle
v^{n}\rangle_{J/\Psi}$ to the amplitude of the leptonic decay. Evidently, this
is not correct. To make formula (29) more consistent with the approach applied
in paper Bodwin:2007ga , it should be modified as follows

$\displaystyle\langle
O_{1}\rangle_{J/\Psi}=\frac{3M_{J/\Psi}^{2}\Gamma_{ee}}{8\pi\alpha^{2}q_{c}^{2}}\biggl{(}\frac{1}{(1-f(\langle
v^{2}\rangle_{J/\Psi}))^{2}}+\frac{16}{3}\frac{\alpha_{s}}{\pi}\biggr{)}.$
(30)

In last formula all problems mentioned above are resolved. Similarly, one can
improve the expression for the operator $\langle O_{1}\rangle_{\eta_{c}}$.
With numerical values of paper Bodwin:2007fz  the factor of the overestimation
is $\sim 1.5$. Taking into account this factor, the cross section $17.6$ fb
obtained in paper Bodwin:2007ga  should be changed to $11.7$ fb.

Similar analysis can be done for the result of paper He:2007te . The
calculation shows that the cross section is overestimated by a factor of $\sim
1.5$. Taking into account this factor the central value of the cross sections
now is $\sim 13.6$ fb.

##  VI Conclusion.

This paper is devoted to the study of the processes $e^{+}e^{-}\to
J/\Psi\eta_{c},J/\Psi\eta_{c}^{\prime},\psi^{\prime}\eta_{c},\psi^{\prime}\eta_{c}^{\prime}$
within Light Cone Formalism (LC). The amplitude for these processes was first
derived in paper Bondar:2004sv . In the present paper the amplitude was
modified so that to archive better accuracy of the calculation and to resolve
some questions raised in Bodwin:2006dm .

To calculate the cross sections of the processes under study one needs the
twist-2 and twist-3 distribution amplitudes (DA). The models of the twist-2
DAs of the $S$-wave charmonia were proposed in papers Braguta:2006wr ;
Braguta:2007fh ; Braguta:2007tq . To get the twist-3 DAs, equations of motion
were applied. It turns out that if one ignores the contribution arising from
higher fock states, the twist-3 DAs can be unambiguously determined.

Using the models of the DAs, the cross sections of the processes
$e^{+}e^{-}\to
J/\Psi\eta_{c},J/\Psi\eta_{c}^{\prime},\psi^{\prime}\eta_{c},\psi^{\prime}\eta_{c}^{\prime}$
were calculated. Within the error of the calculation the results of this study
are in agreement with Belle and BaBar experiments. In addition, the question -
why LC predictions are much greater than the leading order NRQCD predictions -
was studied. Numerical results of the calculation shows that large
disagreement between LC and the leading NRQCD predictions can be attributed to
large contribution of relativistic and radiative corrections. From this
results one can draw a conclusion that in hard exclusive processes
relativistic and radiative corrections play very important role and the
consideration of such processes at the leading NRQCD approximation is
unreliable.

The results of this paper are in agreement with recent NRQCD study of the
process $e^{+}e^{-}\to J/\Psi\eta_{c}$ He:2007te ; Bodwin:2007ga  where the
authors took into account relativistic and one loop radiative corrections.
However, in the present paper it was shown that the results of these papers
are overestimated by a factor 1.5.

###### Acknowledgements.

The author thanks A.K. Likhoded, A.V. Luchinsky for useful discussion. This
work was partially supported by Russian Foundation of Basic Research under
grant 07-02-00417, CRDF grant Y3-P-11-05 and president grant MK-2996.2007.2.

##  Appendix A  Calculation of the constants $f_{P1},f_{P2}$.

To calculate the values of the constants $f_{Pi},~{}i=1,2$ one can apply NRQCD
formalism. At the NLO approximation of NRQCD the constants $f_{Pi}$ can be
written as follows Bodwin:1994jh ; Braaten:1998au

$\displaystyle f_{Pi}^{2}M_{P_{i}}^{2}$ $\displaystyle=$ $\displaystyle\langle
P_{i}|\varphi^{+}\chi|0\rangle\langle
0|\chi^{+}\varphi|P_{i}\rangle\times\biggl{(}1-4\frac{\alpha_{s}}{\pi}-{\langle
v^{2}\rangle_{i}}\biggr{)},$ (31)

where

$\displaystyle\langle v^{2}\rangle_{i}=-\frac{1}{m_{c}^{2}}\frac{\langle
0|\chi^{+}({\overset{\leftrightarrow}{\bf
D}})^{2}\varphi|P_{i}\rangle}{\langle 0|\chi^{+}\varphi|P_{i}\rangle},$ (32)

$P_{i}$ here is the $\eta_{c}$ meson if $i=1$ and the $\eta_{c}^{\prime}$
meson if $i=2$. To continue the calculation one also needs NRQCD expression
for the decay width $\Gamma\left[\eta_{c}\to\gamma\gamma\right]$:

$\displaystyle\Gamma\left[\eta_{c}\to\gamma\gamma\right]=\frac{4\pi
q_{c}^{2}\alpha^{2}}{M_{\eta_{c}}^{3}}\langle\eta_{c}|\varphi^{+}\chi|0\rangle\langle
0|\chi^{+}\varphi|\eta_{c}\rangle\times\biggl{(}1+\frac{\pi^{2}-20}{3}\frac{\alpha_{s}}{\pi}-\frac{\langle
v^{2}\rangle_{\eta_{c}}}{3}\biggr{)}$ (33)

To determine the constant $f_{P1}$ let us express this constant through the
width $\Gamma\left[\eta_{c}\to\gamma\gamma\right]$444In the derivation of (34)
the expression $M_{\eta_{c}}=2M_{c}^{*}+M_{c}^{*}\langle
v^{2}\rangle_{\eta_{c}}$ was used.

$\displaystyle f_{P1}^{2}=\frac{M_{c}^{*}}{2\pi
q_{c}^{2}\alpha^{2}}\Gamma\left[\eta_{c}\to\gamma\gamma\right]\biggl{(}1-\frac{\pi^{2}-8}{3}\frac{\alpha_{s}}{\pi}-\frac{1}{6}\langle
v^{2}\rangle_{\eta_{c}}\biggr{)},$ (34)

where $M_{c}^{*}$ is the mole mass of $c$-quark. The value of the constant
$f_{P1}$ will be calculated with the following set of parameters:
$\Gamma\left[\eta_{c}\to\gamma\gamma\right]=7.2\pm 0.7\pm 2.0$ Yao:2006px ,
$\alpha_{s}=0.25$, $\langle v^{2}\rangle_{\eta_{c}}=0.25$ Bodwin:2006dn ,
$\langle v^{2}\rangle_{\psi^{\prime}}=0.54$ Braguta:2007tq , $M_{c}^{*}=1.4\pm
0.2$ GeV. To estimate the error of the calculation one should take into
account that within NRQCD the constant is double series in the relativistic
and radiative corrections. At the NNLO approximation one has the relativistic
corrections $\sim\langle v^{2}\rangle^{2}$, the radiative corrections to the
short distance coefficient of the operator $\langle
0|\chi^{+}(\vec{\sigma}\vec{\epsilon})\varphi|V_{i}(\epsilon)\rangle$
$\sim\alpha_{s}^{2}$ and the radiative corrections to the short distance
coefficient of the operator $\langle
0|\chi^{+}(\vec{\sigma}\vec{\epsilon})({\overset{\leftrightarrow}{\bf
D}})^{2}\varphi|V_{i}(\epsilon)\rangle$ that can be estimated as
$\sim\alpha_{s}\langle v^{2}\rangle$. In addition, there is an experimental
uncertainty in the measurement of $\Gamma\left[\eta_{c}\to\gamma\gamma\right]$
and the uncertainty in the $m_{c}$. Adding all these uncertainties in
quadrature one can estimate the error of the calculation. Thus one gets

$\displaystyle f_{P1}^{2}=0.139\pm 0.048~{}~{}\mbox{GeV}^{2}.$ (35)

Unfortunately, one cannot apply formula (34) to get the value of the $f_{P2}$.
Since today only the product
$\Gamma\left[\eta_{c}^{\prime}\to\gamma\gamma\right]\times
Br(\eta_{c}^{\prime}\to K^{0}_{S}K^{\pm}\pi^{0})$ has been measured Yao:2006px
and there is no model independent way to determine
$\Gamma\left[\eta_{c}^{\prime}\to\gamma\gamma\right]$.

To estimate the value the $f_{P2}$ one could use the value of the constants
$f_{V2}$ of $f_{T2}(M_{J/\Psi})$, which are equal up to the relativistic
correction and radiative corrections. In this paper the constant
$f_{T2}(M_{J/\Psi})$ will be taken to get the value of $f_{P2}$

$\displaystyle f^{2}_{P2}=0.068\pm 0.040~{}~{}\mbox{GeV}^{2}.$ (36)

##  Appendix B Models for the distribution amplitudes.

In this section models and renormalization group evolution of the DAs needed
in the calculation will be considered. The DAs of the vector meson $V$ are
defined as follows Bondar:2004sv

$\displaystyle\langle
V(p,\epsilon)|\bar{Q}_{\alpha}(z)Q_{\beta}(-z)|0\rangle_{\mu}=\frac{f_{V}M_{V}}{4}\int_{0}^{1}dxe^{i(pz)(2x-1)}\biggl{\\{}\biggl{(}\hat{\epsilon}-\hat{p}\frac{({\epsilon}z)}{(pz)}\biggr{)}V_{\perp}(x,\mu)+\hat{p}\frac{({\epsilon}z)}{(pz)}V_{L}(x,\mu)+$
$\displaystyle\frac{f_{T}(\mu)}{f_{V}}\frac{1}{M_{V}}~{}\sigma_{\mu\nu}{\epsilon}^{\mu}p^{\nu}~{}V_{T}(x,\mu)+\frac{1}{8}\biggl{(}1-\frac{f_{T}(\mu)}{f_{V}}\frac{2M_{c}(\mu)}{{M_{V}}}\biggr{)}~{}e_{\mu\nu\sigma\rho}\gamma^{\mu}\gamma_{5}{\epsilon}^{\nu}p^{\sigma}z^{\rho}~{}V_{A}(x,\mu)\biggr{\\}}_{\beta\alpha}.$
(37)

The DAs of the pseudoscalar meson are defined as

$\displaystyle\langle
P(p)|\bar{Q}_{\alpha}(z)Q_{\beta}(-z)|0\rangle_{\mu}=i\frac{f_{P}M_{V}}{4}\int_{0}^{1}dxe^{i(pz)(2x-1)}\biggl{\\{}\frac{\hat{p}\gamma_{5}}{M_{P}}P_{A}(x,\mu)-\frac{M_{P}}{2M_{c}(\mu)}\gamma_{5}P_{P}(x,\mu)\biggr{\\}}_{\beta\alpha}.$
(38)

where $x$ is the fraction of momentum of the meson $V$ carried by quark, the
constants $f_{T}(\mu),f_{V}$ are defined in equations (18), $M_{c}(\mu)$ is
the running mass of $c$-quark. Expressions (37, 38) are defined at the scale
$\mu$.

The models of the leading twist DAs $V_{L}(x,\mu),V_{T}(x,\mu),P_{A}(x,\mu)$
were proposed in papers Braguta:2006wr ; Braguta:2007fh ; Braguta:2007tq .
According to these models the functions
$P_{A}(x,\mu\sim\overline{M}_{c}),V_{T}(x,\mu\sim\overline{M}_{c}),V_{L}(x,\mu\sim\overline{M}_{c})$
are equal to the function $\phi_{1S}(x,\mu\sim\overline{M}_{c})$ for the $1S$
state mesons and to the function $\phi_{2S}(x,\mu\sim\overline{M}_{c})$ for
the $2S$ state mesons, which have the form Chernyak:1983ej

$\displaystyle\phi_{1S}(x,\mu\sim\overline{M}_{c})$ $\displaystyle\sim$
$\displaystyle\varphi_{as}(x)\mbox{ Exp}\biggl{[}-\frac{b}{4x(1-x)}\biggr{]},$
$\displaystyle\phi_{2S}(x,\mu\sim\overline{M}_{c})$ $\displaystyle\sim$
$\displaystyle\varphi_{as}(x)(a+(2x-1)^{2})\mbox{
Exp}\biggl{[}-\frac{b}{4x(1-x)}\biggr{]},$ (39)

where $\varphi_{as}(x)=6x(1-x)$ is the asymptotic function. For the $1S$
charmonium states the constant $b$ can vary within the interval $3.8\pm 0.7$.
For the $2S$ charmonium states the constants $a$ and $b$ can vary within the
intervals $0.03^{+0.32}_{-0.03}$ and $2.5^{+3.2}_{-0.8}$ correspondingly. The
renormalization group evolution of the leading twist DAs is well known
Chernyak:1983ej  and it can be written in the form

$\displaystyle\phi(x,\mu)=6x(1-x)\biggl{[}1+\sum_{n=2,4..}\biggl{(}\frac{\alpha_{s}(\mu)}{\alpha_{s}(\mu_{0})}\biggr{)}^{\epsilon_{n}/{b_{0}}}a_{n}(\mu_{0})C_{n}^{3/2}(2x-1)\biggr{]},$
(40)

where $a_{n}(\mu_{0})$ is the coefficient of the expansion in Gegenbauer
polynomials $C_{n}^{3/2}(z)$ at scale $\mu_{0}$, the constant $b_{0}=25/3$,
the anomalous dimensions $\epsilon_{n}$ for the functions $V_{L},P_{A}$ are
defined as

$\displaystyle\epsilon_{n}$ $\displaystyle=$
$\displaystyle\frac{4}{3}\biggl{(}1-\frac{2}{(n+1)(n+2)}+4\sum_{j=2}^{n+1}\frac{1}{j}\biggr{)},$
(41)

the anomalous dimension $\epsilon_{n}$ for the function $V_{T}$

$\displaystyle\epsilon_{n}$ $\displaystyle=$
$\displaystyle\frac{4}{3}\biggl{(}4\sum_{j=2}^{n+1}\frac{1}{j}\biggr{)}.$ (42)

Further let us consider the twist-3 DA $P_{P}(x,\mu)$. It turns out the if one
ignores higher fock states it is possible to connect the twist-3 DA
$P_{P}(x,\mu)$ to the twist-2 DA $P_{A}(x,\mu)$ using equations of motion
Ball:1998je

$\displaystyle\langle\xi^{n}\rangle_{P}=\delta_{n0}+\frac{n-1}{n+1}\biggl{(}\langle\xi^{n-2}\rangle_{P}-r(\mu)\langle\xi^{n-2}\rangle_{A}\biggr{)},$
(43)

where $\xi=2x-1$, $\langle\xi^{n}\rangle_{P,A}$ are the moments of the DAs
$P_{P}$ and $P_{A}$ correspondingly, $r(\mu)=4M_{c}(\mu)^{2}/M_{\eta}^{2}$. To
solve these equations one can expand the $P_{P}$ in a series of Gegenbauer
polynomials $G^{1/2}_{n}(z)$ Braun:1989iv

$\displaystyle
P_{P}(x,\mu)=\biggl{[}1+\sum_{n=2,4..}b_{n}(\mu)C_{n}^{1/2}(2x-1)\biggr{]}.$
(44)

Substituting expressions for the DAs $P_{P}$ and $P_{A}$ (40) and (44) to
equations of motion (43), one can solve these equations recursively and relate
the coefficients $b_{n}(\mu)$ to the coefficients $a_{n}(\mu)$. For instance,
for the first tree coefficients one can get the following formulas

$\displaystyle
b_{2}(\mu)=-\frac{5}{2}r(\mu),~{}~{}~{}b_{4}(\mu)=-\frac{27}{20}r(\mu)-\frac{81}{10}r~{}a_{2}(\mu),~{}~{}~{}b_{6}(\mu)=-\frac{13}{14}r(\mu)-\frac{39}{7}r(\mu)a_{2}(\mu)-\frac{195}{14}r(\mu)a_{4}(\mu).$
(45)

It is clear that relativistic motion in the nonrelativistic system must be
strongly suppressed. The behavior of DA in the end point region $x\sim 0,1$ is
determined by the relativistic motion. So, the DA of nonrelativistic system
must be strongly suppressed in the end point region. As it was shown in paper
Braguta:2006wr  for the leading twist DAs such suppression can be achieved if
there is fine tuning of the coefficients $a_{n}(\mu)$ at the scale
$\mu\sim\overline{M}_{c}$. Similarly, to get the suppression in the end point
region for the function $P_{P}$ one requires the fine tuning in the
coefficients $b_{n}(\overline{M}_{c})$, which can be achieved through the fine
tuning of the constants $r(\overline{M}_{c})$ and $a_{n}(\overline{M}_{c})$.
If we put $r(\overline{M}_{c})=4\overline{M}_{c}^{2}/M_{\eta}^{2}$ the fine
tuning for the higher moments will be broken what leads to large relativistic
motion in the end point region. To avoid this problem it will be assumed that
$r(\overline{M}_{c})$ is a free parameter, which will be adjusted through the
requirement that the moments of the $P_{P}(x,\mu\sim\overline{M}_{c})$ must be
equal to the moments of the $P_{A}(x,\mu\sim\overline{M}_{c})$ to the leading
order in relative velocity expansion. Using the last statement and (43) one
can obtain the expression for the $r(\overline{M}_{c})$

$\displaystyle
r(\overline{M}_{c})=1-3\langle\xi^{2}\rangle_{A}+9(\langle\xi^{2}\rangle_{A})^{2}-5\langle\xi^{4}\rangle_{A}-27(\langle\xi^{2}\rangle_{A})^{3}+30\langle\xi^{2}\rangle_{A}\langle\xi^{4}\rangle_{A}-7\langle\xi^{6}\rangle_{A}+O(v^{8}).$
(46)

All moments in this expression are taken at scale $\mu\sim\overline{M}_{c}$.
The exact expression for the $r(\overline{M}_{c})$ can be written in the
following form

$\displaystyle\frac{2}{r(\overline{M}_{c})}=\int_{-1}^{1}d\xi\frac{1+\xi^{2}}{(1-\xi^{2})^{2}}P_{A}\biggl{(}\frac{1+\xi}{2},\mu\sim\overline{M}_{c}\biggr{)},$
(47)

where $\xi=2x-1$.

The same approach can be applied to the twist-3 DAs $V_{A},V_{\perp}$. These
functions can be expanded in a series of Gegenbauer polynomials Ball:1998sk :

$\displaystyle
V_{\perp}(x,\mu)=\biggl{[}1+\sum_{n=2,4..}c_{n}(\mu)C_{n}^{1/2}(2x-1)\biggr{]},$
$\displaystyle
V_{A}(x,\mu)=6x(1-x)\biggl{[}1+\sum_{n=2,4..}d_{n}(\mu)C_{n}^{3/2}(2x-1)\biggr{]}.$
(48)

The coefficients $c_{n}(\mu)$ and $d_{n}(\mu)$ can be related to the
coefficients $a_{n}(\mu)$ of the functions $V_{L},V_{T}$ through the equations
of motion Ball:1998sk

$\displaystyle(n+1)\langle\xi^{n}\rangle_{\perp}=\langle\xi^{n}\rangle_{L}+\frac{n(n-1)}{2}(1-\delta(\mu))\langle\xi^{n-2}\rangle_{A},$
(49)
$\displaystyle\frac{1}{2}(n+2)(1-\delta(\mu))\langle\xi^{n}\rangle_{A}=\langle\xi^{n}\rangle_{\perp}-\delta(\mu)\langle\xi^{n}\rangle_{T},$

where $\langle\xi^{n}\rangle_{L,T,\perp,A}$ are the moments of the DAs
$V_{L},V_{T},V_{\perp},V_{A}$,
$\delta(\mu)=2f_{T}(\mu)/f_{V}M_{c}(\mu)/M_{V}$. Similarly, to the constant
$r(\overline{M}_{c})$ the value of the constant $\delta(\overline{M}_{c})$
must be adjusted. We adjust it through the requirement that the moments of the
DA $V_{\perp}$ must be equal to the moments of the $V_{L}$ to the leading
order approximation in relative velocity expansion. Thus one can get

$\displaystyle\delta(\overline{M}_{c})=1-\biggl{(}2\langle\xi^{2}\rangle_{L}-6(\langle\xi^{2}\rangle_{L})^{2}+4\langle\xi^{4}\rangle_{L}-22\langle\xi^{2}\rangle_{L}\langle\xi^{4}\rangle_{L}+18(\langle\xi^{2}\rangle_{L})^{3}+6\langle\xi^{6}\rangle_{L}\biggr{)}+O(v^{8}).$
(50)

It should be noted that with this value of the $\delta(\overline{M}_{c})$ the
moments of the function $V_{A}$ are equal to the moments of the $V_{L}$ to the
leading order in relative velocity. The exact expression for the
$\delta(\overline{M}_{c})$ is

$\displaystyle\delta(\overline{M}_{c})=\frac{\int_{-1}^{1}\frac{d\xi}{1-\xi^{2}}V_{L}\biggl{(}\frac{1+\xi}{2},\mu\sim\overline{M}_{c}\biggr{)}}{\int_{-1}^{1}{d\xi}\frac{1+\xi^{2}}{(1-\xi^{2})^{2}}V_{T}\biggl{(}\frac{1+\xi}{2},\mu\sim\overline{M}_{c}\biggr{)}}$
(51)

## References

  * (1) E. Braaten and J. Lee, Phys. Rev. D 67, 054007 (2003) [arXiv:hep-ph/0211085];
  * (2) K. Y. Liu, Z. G. He and K. T. Chao, Phys. Lett. B 557, 45 (2003) [arXiv:hep-ph/0211181];
  * (3) K. Y. Liu, Z. G. He and K. T. Chao, arXiv:hep-ph/0408141.
  * (4) K. Abe et al. [Belle Collaboration], Phys. Rev. Lett. 89, 142001 (2002), [arXiv:hep-ex/0205104].
  * (5) K. Abe et al. [Belle Collaboration], Phys. Rev. D 70, 071102 (2004), hep-ex/0407009.
  * (6) B. Aubert [BABAR Collaboration], hep-ex/0506062.
  * (7) G. T. Bodwin, J. Lee and E. Braaten, Phys. Rev. Lett. 90, 162001 (2003) [arXiv:hep-ph/0212181].
  * (8) G. T. Bodwin, J. Lee and E. Braaten, Phys. Rev. D 67, 054023 (2003) [Erratum-ibid. D 72, 099904 (2005)] [arXiv:hep-ph/0212352].
  * (9) A. V. Luchinsky, arXiv:hep-ph/0301190.
  * (10) Y. J. Zhang, Y. j. Gao and K. T. Chao, Phys. Rev. Lett. 96, 092001 (2006) [arXiv:hep-ph/0506076].
  * (11) B. Gong and J. X. Wang, Phys. Rev. D 77, 054028 (2008) [arXiv:0712.4220 [hep-ph]].
  * (12) Y. J. Zhang, Y. Q. Ma and K. T. Chao, arXiv:0802.3655 [hep-ph].
  * (13) J. P. Ma and Z. G.Si, Phys. Rev. D 70, 074007 (2004), [arXiv:hep-ph/0405111].
  * (14) A. E. Bondar and V. L. Chernyak, Phys. Lett. B 612, 215 (2005) [arXiv:hep-ph/0412335].
  * (15) V. V. Braguta, A. K. Likhoded and A. V. Luchinsky, Phys. Rev. D 72, 094018 (2005) [arXiv:hep-ph/0506009].
  * (16) V. V. Braguta, A. K. Likhoded and A. V. Luchinsky, Phys. Rev. D 72, 074019 (2005) [arXiv:hep-ph/0507275].
  * (17) V. V. Braguta, A. K. Likhoded and A. V. Luchinsky, Phys. Lett. B 635, 299 (2006) [arXiv:hep-ph/0602047].
  * (18) D. Ebert and A. P. Martynenko, Phys. Rev. D 74, 054008 (2006) [arXiv:hep-ph/0605230].
  * (19) H. M. Choi and C. R. Ji, Phys. Rev. D 76, 094010 (2007) [arXiv:0707.1173 [hep-ph]].
  * (20) A. V. Berezhnoy, arXiv:hep-ph/0703143.
  * (21) D. Ebert, R. N. Faustov, V. O. Galkin and A. P. Martynenko, arXiv:0803.2124 [hep-ph].
  * (22) J. Lee, J. Korean Phys. Soc. 45, S354 (2004) [arXiv:hep-ph/0312251].
  * (23) Y. J. Zhang, Q. Zhao and C. F. Qiao, arXiv:0806.3140 [hep-ph].
  * (24) A. V. Berezhnoy and A. K. Likhoded, Phys. Atom. Nucl. 70, 478 (2007) [arXiv:hep-ph/0602041].
  * (25) G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D 51, 1125 (1995) [Erratum-ibid. D 55, 5853 (1997)] [arXiv:hep-ph/9407339].
  * (26) Z. G. He, Y. Fan and K. T. Chao, Phys. Rev. D 75, 074011 (2007) [arXiv:hep-ph/0702239].
  * (27) G. T. Bodwin, J. Lee and C. Yu, Phys. Rev. D 77, 094018 (2008) [arXiv:0710.0995 [hep-ph]].
  * (28) G.T. Bodwin, Talk presented at the conference Quark Confinement and Hadron Spectrum, Mainz (Germany), 1-6 September 2008.
  * (29) G. P. Lepage and S. J. Brodsky, Phys. Rev. D 22, 2157 (1980).
  * (30) V. L. Chernyak and A. R. Zhitnitsky, Phys. Rept. 112, 173 (1984).
  * (31) G. T. Bodwin, D. Kang and J. Lee, Phys. Rev. D 74, 114028 (2006) [arXiv:hep-ph/0603185].
  * (32) J. P. Ma and Z. G. Si, Phys. Lett. B 647, 419 (2007) [arXiv:hep-ph/0608221].
  * (33) V. V. Braguta, A. K. Likhoded and A. V. Luchinsky, Phys. Lett. B 646, 80 (2007) [arXiv:hep-ph/0611021].
  * (34) V. V. Braguta, Phys. Rev. D 75, 094016 (2007) [arXiv:hep-ph/0701234].
  * (35) V. V. Braguta, arXiv:0709.3885 [hep-ph].
  * (36) T. Feldmann and G. Bell, arXiv:0711.4014 [hep-ph].
  * (37) G. Bell and T. Feldmann, JHEP 0804, 061 (2008) [arXiv:0802.2221 [hep-ph]].
  * (38) C. W. Hwang, arXiv:0811.0648 [hep-ph].
  * (39) V. M. Braun, G. P. Korchemsky and D. Mueller, Prog. Part. Nucl. Phys. 51, 311 (2003) [arXiv:hep-ph/0306057].
  * (40) V. V. Braguta, Phys. Rev. D 78, 054025 (2008) [arXiv:0712.1475 [hep-ph]].
  * (41) G. T. Bodwin, D. Kang and J. Lee, Phys. Rev. D 74, 014014 (2006) [arXiv:hep-ph/0603186].
  * (42) G. T. Bodwin, H. S. Chung, D. Kang, J. Lee and C. Yu, Phys. Rev. D 77, 094017 (2008) [arXiv:0710.0994 [hep-ph]].
  * (43) E. Braaten and Y. Q. Chen, Phys. Rev. D 57, 4236 (1998) [Erratum-ibid. D 59, 079901 (1999)] [arXiv:hep-ph/9710357].
  * (44) W. M. Yao et al. [Particle Data Group], J. Phys. G 33, 1 (2006).
  * (45) P. Ball, JHEP 9901, 010 (1999) [arXiv:hep-ph/9812375].
  * (46) V. M. Braun and I. E. Filyanov, Z. Phys. C 48, 239 (1990) [Sov. J. Nucl. Phys. 52, 126 (1990 YAFIA,52,199-213.1990)].
  * (47) P. Ball, V. M. Braun, Y. Koike and K. Tanaka, Nucl. Phys. B 529, 323 (1998) [arXiv:hep-ph/9802299].