Eur. Phys. J. C (2019) 79 :323 Page 5 of 23 323
functions, so a numerical approach must be employed. There
are many numerical methods to solve the Cornell potential
in the literature. For its simplicity and robustness, we will
use the Numerov algorithm [22,23] (also called Cowell’s
method). Further details can be found in Appendix A. We
have implemented the Numerov algorithm and the matrix
element computations in a Fortran 2008 code [24], which
contains both public and in-house routines, and is used to
predict the mass of the bound states as well as to perform fits
to data.
The potential in Eq. (1) by itself is nonetheless unable
to describe the Υ(1S) − η
b
(1S) mass splitting, as it does
not depend on the spin of the quarkonium state. Further-
more, it provides degenerate masses for the χ
bJ
multiplet
(with J ={0, 1, 2}). Hence, in order to describe with higher
accuracy the low-lying bottomonium spectrum, one has to
add 1/m
2
Q
terms to the Cornell static potential that take
into account the spin–spin, spin–orbit and tensor interac-
tions, breaking the present degeneracy [25]. Such contribu-
tions arise from the leading relativistic corrections of the
t-channel gluon and confinement interactions, giving the fol-
lowing terms [5,26]:
3,4
V
OGE
SS
(r) =
8α
Cornell
s
9m
2
Q
r
2
(S
1
· S
2
)δ(r),
V
OGE
LS
(r) =
2α
Cornell
s
m
2
Q
r
3
(L · S ),
V
CON
LS
(r) =−
σ
2m
2
Q
r
(L · S ),
V
OGE
T
(r) =
α
Cornell
s
3m
2
Q
r
3
S
12
, (13)
where S = S
1
+ S
2
is the total spin, L the relative orbital
momentum and S
12
the tensor operator of the Q Q bound
state, defined as
S
12
= 2 (S
1
·ˆr)(S
2
·ˆr) − (S
1
· S
2
), (14)
with ˆr = r/r. Analytical expressions can be found for the
spin–spin, spin–orbital and tensor operators,
S
1
· S
2
=
1
4
[2 s (s + 1) − 3] ,
(15a)
L · S=
1
2
[ j ( j + 1) − (+ 1) − s (s + 1)],
(15b)
S
12
=
2 [2 (+ 1) s (s + 1) − 3 L · S−6 L · S
2
]
(2 − 1)(2 + 3)
,
(15c)
3
Here the LS term from confinement is obtained by the exchange of a
scalar particle.
4
We denote with the superscript “OGE” the terms arising from gluon
exchanges and with “CON” those coming from the confinement inter-
action.
with ( j,,s) the quantum numbers of the Q Q state. The
specific values for the quantum numbers considered in this
article are given in Table 1.
Given the assumed large mass of the heavy quark, the
above terms are expected to be small. Therefore, their con-
tribution to the bottomonium mass will be incorporated to
the model using first-order perturbation theory. For consis-
tency, if perturbation theory is to be used to second order, one
needs to include the 1/m
4
Q
potential. In practice one needs
to take matrix elements of the operators in Eq. (13). To that
end we use the angular decomposition in Eq. (3) and write
the three-dimensional integration in spherical coordinates
d
3
r = r
2
dr dΩ such that the angular integrations are carried
out with ordinary angular momentum algebra [see Eq. (15)]
and the radial integral becomes a one-dimensional matrix ele-
ment of the reduced wave function u
n
(r), given that the Jaco-
bian factor r
2
cancels when writing |R
n
(r)|
2
=|u
n
(r)|
2
/r
2
as in Eq. (6).
3 Fitting the Cornell model to experimental data
In this section we present the results from χ
2
fitsofthe
Cornell model parameters to charmonium and bottomonium
experimental data. We will restrict ourselves to n
p
≤ 2, since
for higher values of the principal quantum number one needs
to include string-breaking effects, as can be seen in Fig. 4.
The case of bottomonium will serve as a proof of concept for
our calibration, that is, it will show that the three parameters
of the model can be determined from fits to the 8 lowest-lying
bottomonium bound states. This is crucial since in QCD we
can only reliably predict these within perturbation theory, as
argued in Ref. [27]. In many phenomenological applications
of the quark model approach, in which more sophisticated
Hamiltonians are used, the parameters of the potential are
not obtained by a full fledged χ
2
minimization, but sim-
ply adjusted by a rough comparison to data plus physically
motivated priors. This procedure is in many cases justified,
since there is a large amount of data that the model aims
to describe, but its accuracy might vary widely for different
observables. Given that it is very hard to add a theoretical
covariance matrix to the χ
2
(the matrix does not provide a
method to quantify its “modeling” error), such a fit could
lead to biased model parameter values. For the simpler situa-
tion of the naïve Cornell model, which seeks to describe only
the low-lying quarkonium spectrum, we show that the fit is
indeed possible, and we shall see that no bias is observed.
The minimization of the χ
2
is carried out using the For-
tran 77 package MINUIT [28]. We have checked that the
algorithm is very effective in finding the minimum, but due
to the strong correlation between the various fit parameters
[see Eq. (17)] it does not estimate the covariance matrix cor-
rectly. To solve this problem we compute a three-dimensional
123