Lattice Gluon Propagator and One-Gluon-Exchange Potential

TL;DR

This paper investigates the interquark potential using a nonperturbative gluon propagator from lattice QCD, showing the limitations of the one-gluon-exchange approximation and proposing an improved potential for quarkonium spectra.

Contribution

It introduces a nonperturbative lattice-based gluon propagator into the OGE approximation and enhances the potential with a linear term to better match quarkonium spectra.

Findings

The VLGP potential is non-confining and insufficient for the infrared sector.

Adding a linear term improves the potential's agreement with observed spectra.

Estimated interquark distances for charmonium and bottomonium states.

Abstract

We consider the interquark potential in the one-gluon-exchange (OGE) approximation, using a fully nonperturbative gluon propagator from large-volume lattice simulations. The resulting VLGP potential is non-confining, showing that the OGE approximation is not sufficient to describe the infrared sector of QCD. Nevertheless, it represents an improvement over the perturbative (Coulomb-like) potential, since it allows the description of a few low-lying bound states of charmonium and bottomonium. In order to achieve a better description of these spectra, we add to VLGP a linearly growing term. The obtained results are comparable to the corresponding ones in the Cornell-potential case. As a byproduct of our study, we estimate the interquark distance for the considered charmonium and bottomonium states.

Figures

Tables

Table 1. Table 1: Ranges of parameter values and iterative step sizes used to obtain the charmonium and bottomonium eigenenergies E 𝐸 E with the V L ​ G ​ P subscript 𝑉 𝐿 𝐺 𝑃 V_{LGP} potential. We also show the range and step of integration for the radial distance r 𝑟 r .

	Charmonium	Bottomonium	Step
E	$2.50$ to $4.50\mathrm{GeV}$	$8.50$ to $12.50\mathrm{GeV}$	$0.04\mathrm{GeV}$
$m_c$	$1.00$ to $2.25\mathrm{GeV}$		$0.001\mathrm{GeV}$
$m_b$		$4.00$ to $6.00\mathrm{GeV}$	$0.001\mathrm{GeV}$
${\alpha }_s$	$0.10$ to $1.00$		$0.01$
$r$	$200.0$ to $0.0{\mathrm{GeV}}^{-1}$		$0.05{\mathrm{GeV}}^{-1}$

Table 2. Table 2: Ranges of parameter values and iterative step sizes used to obtain the charmonium and bottomonium eigenenergies E 𝐸 E with the V L ​ G ​ P + L subscript 𝑉 𝐿 𝐺 𝑃 𝐿 V_{LGP+L} and Cornell potentials. For α s subscript 𝛼 𝑠 \alpha_{s} , taken as a fixed parameter, we list the corresponding perturbative values (see Section II ). We also show the range and step of integration for the radial distance r 𝑟 r .

	Charmonium	Bottomonium	Step
E	$2.00$ to $6.00\mathrm{GeV}$	$8.50$ to $12.50\mathrm{GeV}$	$0.04\mathrm{GeV}$
$m_c$	$1.00$ to $2.00\mathrm{GeV}$		$0.01\mathrm{GeV}$
$m_b$		$4.15$ to $4.85\mathrm{GeV}$	$0.01\mathrm{GeV}$
$\sigma$	$0.10$ to $0.50{\mathrm{GeV}}^2$		$0.01{\mathrm{GeV}}^2$
${\alpha }_s$	$0.2663$	$0.1843$
$r$	$50.0$ to $0.0{\mathrm{GeV}}^{-1}$		$0.0125{\mathrm{GeV}}^{-1}$

Table 3. Table 3: Experimental spectrum of charmonium states and spin-averaged values using two different methods (see text).

Particle Name	Mass ($\mathrm{GeV}$)	$J^{PC}$	$l$	MAV1 ($\mathrm{GeV}$)	MAV2 ($\mathrm{GeV}$)
${\eta }_c\left(1S\right)$	$2.9836\left(7\right)$	$0^{-+}$	0	$\mathrm{3.068\hspace{0.17em}59}\left(18\right)$	$3.040\left(57\right)$
$J/\psi \left(1S\right)$	$\mathrm{3.096\hspace{0.17em}916}\left(11\right)$	$1^{--}$
${\chi }_{c0}\left(1P\right)$	$\mathrm{3.414\hspace{0.17em}75}\left(31\right)$	$0^{++}$	1	$\mathrm{3.525\hspace{0.17em}32}\left(11\right)$	$3.485\left(70\right)$
${\chi }_{c1}\left(1P\right)$	$\mathrm{3.510\hspace{0.17em}66}\left(7\right)$	$1^{++}$
$h_c\left(1P\right)$	$\mathrm{3.525\hspace{0.17em}38}\left(11\right)$	$1^{+-}$
${\chi }_{c2}\left(1P\right)$	$\mathrm{3.556\hspace{0.17em}20}\left(9\right)$	$2^{++}$
${\eta }_c\left(2S\right)$	$3.6394\left(13\right)$	$0^{-+}$	0	$\mathrm{3.674\hspace{0.17em}43}\left(34\right)$	$3.663\left(23\right)$
$\psi \left(2S\right)$	$\mathrm{3.686\hspace{0.17em}109}\left(14\right)$	$1^{--}$
$\psi \left(3770\right)$	$\mathrm{3.773\hspace{0.17em}15}\left(33\right)$	$1^{--}$	0 or 2	$\mathrm{3.773\hspace{0.17em}15}\left(33\right)$	$\mathrm{3.773\hspace{0.17em}15}\left(33\right)$
$X\left(3872\right)$	$\mathrm{3.871\hspace{0.17em}69}\left(17\right)$	$1^{++}$	1	$\mathrm{3.871\hspace{0.17em}69}\left(17\right)$	$\mathrm{3.871\hspace{0.17em}69}\left(17\right)$
${\chi }_{c0}\left(2P\right)$	$3.9184\left(19\right)$	$0^{++}$	1	$3.9257\left(25\right)$	$3.9228\left(44\right)$
${\chi }_{c2}\left(2P\right)$	$3.9272\left(26\right)$	$2^{++}$
$\psi \left(4040\right)$	$4.039\left(1\right)$	$1^{--}$	0 or 2	$4.039\left(1\right)$	$4.039\left(1\right)$
$\psi \left(4160\right)$	$4.191\left(5\right)$	$1^{--}$	0 or 2	$4.191\left(5\right)$	$4.191\left(5\right)$
$X\left(4260\right)$	$4.251\left(9\right)$	$1^{--}$	0 or 2	$4.251\left(9\right)$	$4.251\left(9\right)$
$X\left(4360\right)$	$4.361\left(13\right)$	$1^{--}$	0 or 2	$4.361\left(13\right)$	$4.361\left(13\right)$
$\psi \left(4415\right)$	$4.421\left(4\right)$	$1^{--}$	0 or 2	$4.421\left(4\right)$	$4.421\left(4\right)$
$X\left(4660\right)$	$4.664\left(12\right)$	$1^{--}$	0 or 2	$4.664\left(12\right)$	$4.664\left(12\right)$

Table 4. Table 4: Experimental spectrum of bottomonium states and spin-averaged values using two different methods (see text). We include the unconfirmed state η b ​ ( 1 ​ S ) subscript 𝜂 𝑏 1 𝑆 \eta_{b}(1S) .

Particle Name	Mass ($\mathrm{GeV}$)	$J^{PC}$	$l$	MAV1 ($\mathrm{GeV}$)	MAV2 ($\mathrm{GeV}$)
${\eta }_b\left(1S\right)$	$9.3980\left(32\right)$	$0^{-+}$	0	$9.4447\left(10\right)$	$9.429\left(31\right)$
${\rm Y}\left(1S\right)$	$\mathrm{9.460\hspace{0.17em}30}\left(26\right)$	$1^{--}$
${\chi }_{b0}\left(1P\right)$	$\mathrm{9.859\hspace{0.17em}44}\left(42\right)$	$0^{++}$	1	$\mathrm{9.899\hspace{0.17em}73}\left(46\right)$	$9.886\left(26\right)$
${\chi }_{b1}\left(1P\right)$	$\mathrm{9.892\hspace{0.17em}78}\left(26\right)$	$1^{++}$
$h_b\left(1P\right)$	$9.8993\left(10\right)$	$1^{+-}$
${\chi }_{b2}\left(1P\right)$	$\mathrm{9.912\hspace{0.17em}21}\left(26\right)$	$2^{++}$
${\rm Y}\left(2S\right)$	$\mathrm{10.023\hspace{0.17em}26}\left(31\right)$	$1^{--}$	0	$\mathrm{10.023\hspace{0.17em}26}\left(31\right)$	$\mathrm{10.023\hspace{0.17em}26}\left(31\right)$
${\rm Y}\left(1D\right)$	$10.1637\left(14\right)$	$2^{--}$	2	$10.1637\left(14\right)$	$10.1637\left(14\right)$
${\chi }_{b0}\left(2P\right)$	$10.2325\left(4\right)$	$0^{++}$	1	$\mathrm{10.260\hspace{0.17em}24}\left(46\right)$	$10.251\left(18\right)$
${\chi }_{b1}\left(2P\right)$	$\mathrm{10.255\hspace{0.17em}46}\left(22\right)$	$1^{++}$
${\chi }_{b2}\left(2P\right)$	$\mathrm{10.268\hspace{0.17em}650}\left(22\right)$	$2^{++}$
${\rm Y}\left(3S\right)$	$10.3552\left(2\right)$	$1^{--}$	0	$10.3552\left(2\right)$	$10.3552\left(2\right)$
${\chi }_b\left(3P\right)$	$10.534\left(9\right)$	${?}^{?+}$	1	$10.534\left(9\right)$	$10.534\left(9\right)$
${\rm Y}\left(4S\right)$	$10.5794\left(12\right)$	$1^{--}$	0	$10.5794\left(12\right)$	$10.5794\left(12\right)$
${\rm Y}\left(10860\right)$	$10.876\left(11\right)$	$1^{--}$	0 or 2	$10.876\left(11\right)$	$10.876\left(11\right)$
${\rm Y}\left(11020\right)$	$11.019\left(8\right)$	$1^{--}$	0 or 2	$11.019\left(8\right)$	$11.019\left(8\right)$

Table 5. Table 5: Results for the charmonium eigenstates using the V L ​ G ​ P subscript 𝑉 𝐿 𝐺 𝑃 V_{LGP} potential. We leave α s subscript 𝛼 𝑠 \alpha_{s} and m c subscript 𝑚 𝑐 m_{c} as free parameters, obtaining α s = 0.95 ± 0.02 subscript 𝛼 𝑠 plus-or-minus 0.95 0.02 \alpha_{s}=0.95\pm 0.02 and m c = 2.064 − 0.009 + 0.010   GeV subscript 𝑚 𝑐 times superscript subscript 2.064 0.009 0.010 gigaelectronvolt m_{c}=$2.064_{-0.009}^{+0.010}\text{\,}\mathrm{GeV}$ . The states 1S, 2S and 1P were used as inputs in the fits. Long dashes represent states that have not been observed experimentally.

Charmonium Spectrum
$V_{LGP}$
State	Mass	Deviation from average
	($\mathrm{GeV}$)	spin state ($\mathrm{GeV}$)
1S	$3.05\left(11\right)$	$0.01$
1P	$3.64\left(13\right)$	$-0.02$
2S	$3.64\left(13\right)$	$0.16$
1D	$3.64\left(13\right)$	—oii
2P	$4.10\left(13\right)$	$0.18$
3S	$4.12\left(15\right)$	—oii
2D	$4.13\left(18\right)$	—oii
3P	$4.13\left(24\right)$	—oii
4S	$4.13\left(140\right)$	—oii

Table 6. Table 6: Results for the bottomonium eigenstates using the V L ​ G ​ P subscript 𝑉 𝐿 𝐺 𝑃 V_{LGP} potential. We leave α s subscript 𝛼 𝑠 \alpha_{s} and m b subscript 𝑚 𝑏 m_{b} as free parameters, obtaining α s = 0.513 − 0.010 + 0.009 subscript 𝛼 𝑠 subscript superscript 0.513 0.009 0.010 \alpha_{s}=0.513^{+0.009}_{-0.010} and m b = 5.10947 − 0.00014 + 0.00016   GeV subscript 𝑚 𝑏 times subscript superscript 5.10947 0.00016 0.00014 gigaelectronvolt m_{b}=$5.10947^{+0.00016}_{-0.00014}\text{\,}\mathrm{GeV}$ . The states 1S, 2S and 1P were used as inputs in the fits. Long dashes represent states that have not been observed experimentally.

Bottomonium Spectrum
$V_{LGP}$
State	Mass	Deviation from average
	($\mathrm{GeV}$)	spin state ($\mathrm{GeV}$)
1S	$9.43\left(30\right)$	$0.00$
1P	$10.02\left(32\right)$	$0.14$
2S	$10.02\left(32\right)$	$0.00$
1D	$10.02\left(32\right)$	$-0.14$
2P	$10.13\left(32\right)$	$-0.12$
3S	$10.16\left(32\right)$	$-0.19$
3P	$10.22\left(32\right)$	$-0.31$
4S	$10.22\left(32\right)$	$-0.36$
2D	$10.22\left(32\right)$	—oii

Table 7. Table 7: Quark masses and string tension obtained from our preferred fit. These parameters are used to obtain the spectrum in Tables 8 and 9 .

$V_{LGP+L}$	Cornell Potential	Quark Mass
		in Ref. Agashe:2014kda
$m_c=1.16\left(3\right)\mathrm{GeV}$	$m_c={1.11}_{-0.02}^{+0.08}\mathrm{GeV}$	$m_c=1.275\left(25\right)\mathrm{GeV}$
$m_b={4.61}_{-0.01}^{+0.02}\mathrm{GeV}$	$m_b={4.58}_{-0.01}^{+0.04}\mathrm{GeV}$	$m_b\left(\overline{MS}\right)=4.18\left(3\right)\mathrm{GeV}$
$\sigma ={0.23\left(1\right)\mathrm{GeV}}^2$	$\sigma ={{0.26}_{-0.03}^{+0.01}\mathrm{GeV}}^2$	$m_b\left(1S\right)=4.66\left(3\right)\mathrm{GeV}$
${\chi }^2=6.20$	${\chi }^2=12.13$

Table 8. Table 8: Results for the charmonium eigenstates using V L ​ G ​ P + L subscript 𝑉 𝐿 𝐺 𝑃 𝐿 V_{LGP+L} and the Cornell potentials in a constrained fit (see text). Long dashes represent states that have not been observed experimentally.

Charmonium Spectrum
$V_{LGP+L}$
State	Mass	Deviation from average
	($\mathrm{GeV}$)	spin state ($\mathrm{GeV}$)
1S	$2.96\left(11\right)$	$-0.10$
1P	$3.46\left(12\right)$	$-0.07$
2S	$3.69\left(13\right)$	$0.01$
1D	$3.81\left(14\right)$	—oii
2P	$4.02\left(14\right)$	$0.09$
3S	$4.22\left(15\right)$	—oii
2D	$4.31\left(15\right)$	—oii
3P	$4.50\left(16\right)$	—oii
4S	$4.69\left(17\right)$	—oii
Cornell Potential
State	Mass	Deviation from average
	($\mathrm{GeV}$)	spin state ($\mathrm{GeV}$)
1S	$2.93\left(17\right)$	$-0.14$
1P	$3.42\left(19\right)$	$-0.11$
2S	$3.69\left(20\right)$	$0.01$
1D	$3.80\left(21\right)$	—oii
2P	$4.04\left(22\right)$	$0.12$
3S	$4.28\left(24\right)$	—oii
2D	$4.36\left(24\right)$	—oii
3P	$4.58\left(26\right)$	—oii
4S	$4.79\left(27\right)$	—oii

Table 9. Table 9: Results for the bottomonium eigenstates using V L ​ G ​ P + L subscript 𝑉 𝐿 𝐺 𝑃 𝐿 V_{LGP+L} and the Cornell potentials in a constrained fit (see text). Long dashes represent states that have not been observed experimentally.

Bottomonium Spectrum
$V_{LGP+L}$
State	Mass	Deviation from average
	($\mathrm{GeV}$)	spin state ($\mathrm{GeV}$)
1S	$9.47\left(30\right)$	$0.04$
1P	$9.86\left(31\right)$	$-0.03$
2S	$10.00\left(32\right)$	$-0.01$
1D	$10.11\left(160\right)$	$-0.05$
2P	$10.24\left(33\right)$	$-0.01$
3S	$10.37\left(33\right)$	$0.01$
2D	$10.44\left(140\right)$	—oii
3P	$10.56\left(34\right)$	$0.03$
4S	$10.67\left(34\right)$	$0.10$
3D	$10.73\left(220\right)$	—oii
4P	$10.84\left(35\right)$	—oii
Cornell Potential
State	Mass	Deviation from average
	($\mathrm{GeV}$)	spin state ($\mathrm{GeV}$)
1S	$9.49\left(31\right)$	$0.06$
1P	$9.84\left(33\right)$	$-0.04$
2S	$10.00\left(33\right)$	$-0.01$
1D	$10.10\left(34\right)$	$-0.06$
2P	$10.25\left(34\right)$	$0.00$
3S	$10.39\left(35\right)$	$0.03$
2D	$10.45\left(35\right)$	—oii
3P	$10.59\left(36\right)$	$0.05$
4S	$10.72\left(37\right)$	$0.14$
3D	$10.77\left(240\right)$	—oii
4P	$10.89\left(38\right)$	—oii

Table 10. Table 10: Typical interquark distances for charmonium. Errors are expected to be negligible.

Charmonium
State	$V_{LGP+L}$	Cornell Potential
	distance ($\mathrm{fm}$)	distance ($\mathrm{fm}$)
1S	$0.40$	$0.42$
1P	$0.64$	$0.64$
2S	$0.79$	$0.78$
1D	$0.84$	$0.82$
2P	$0.97$	$0.95$
3S	$1.10$	$1.07$
2D	$1.13$	$1.09$
3P	$1.25$	$1.21$
4S	$1.36$	$1.33$

Table 11. Table 11: Typical interquark distances for bottomonium. Errors are expected to be negligible.

Bottomonium
State	$V_{LGP+L}$	Cornell Potential
	distance ($\mathrm{fm}$)	distance ($\mathrm{fm}$)
1S	$0.22$	$0.24$
1P	$0.38$	$0.38$
2S	$0.47$	$0.47$
1D	$0.51$	$0.50$
2P	$0.59$	$0.58$
3S	$0.67$	$0.65$
2D	$0.69$	$0.67$
3P	$0.77$	$0.75$
4S	$0.84$	$0.81$
3D	$0.85$	$0.81$
4P	$0.93$	$0.90$

Equations

$$V(r)=-\frac{4}{3}\frac{{\alpha }_s}{r}+\sigma r,$$

$$S_{fi}\equiv ⟨f|i⟩={\delta }_{fi}+i(2\pi {)}^4{\delta }^{(4)}(Q-P)T_{fi},$$

$$T_{fi}=\frac{1}{(2\pi {)}^6}\frac{m^2}{\sqrt{E_{p_1}{E}_{p_2}{E}_{q_1}{E}_{q_2}}}\left(t_{\text{exch}}+t_{\text{annihil}}\right),$$

$$t_{\text{exch}}=$$

$$\times \overline{v}(p_2,{\sigma }_2){\gamma }^{\nu }v(q_2,{\tau }_2)$$

$$t_{\text{annihil}}=$$

$$\times \overline{u}(q_1,{\tau }_1){\gamma }^{\nu }v(q_2,{\tau }_2).$$

$$u(p,1/2)$$

$$v(p,1/2)$$

$$T_{fi}=\frac{1}{(2\pi {)}^6}\left(t_{\text{exch}}+t_{\text{annihil}}\right).$$

$$t_{\text{exch}}$$

$$=e^2{P}_{00}(k){\delta }_{{\sigma }_1{\tau }_1}{\delta }_{{\sigma }_2{\tau }_2},$$

$$k=p_1-q_1=\left(\begin{array}{cc}0,& \vec{k}\end{array}\right).$$

$$k=\left(\begin{array}{cc}2m,& 0\end{array}\right).$$

$$T_{fi}=\frac{1}{(2\pi {)}^6}\frac{e^2}{{\vec{k}}^2}.$$

$$V(\vec{r})$$

$$=-\frac{e^2}{4\pi r}.$$

$$t_{\text{exch}}$$

$$\times \overline{v}(p_2,{\sigma }_2)c_{2,i}^{†}{\lambda }^b{\gamma }^{\nu }{c}_{2,f}v(q_2,{\tau }_2)$$

$$t_{\text{annihil}}$$

$$\times \overline{u}(q_1,{\tau }_1)c_{1,f}^{†}{\lambda }^b{\gamma }^{\nu }{c}_{2,f}v(q_2,{\tau }_2),$$

$$c_{1,f}^{†}{\lambda }^a{c}_{1,i}{c}_{2,i}^{†}{\lambda }^a{c}_{2,f}=\frac{1}{3}\mathrm{Tr}{\lambda }^a{\lambda }^a=\frac{{\delta }^{aa}}{6}=\frac{4}{3}$$

$$c_{2,i}^{†}{\lambda }^a{c}_{1,i}{c}_{1,f}^{†}{\lambda }^a{c}_{2,f}=\frac{1}{3}(\mathrm{Tr}{\lambda }^a)(\mathrm{Tr}{\lambda }^a)=0.$$

$$P_{\mu \nu }^{ab}=-\frac{g_{\mu \nu }{\delta }^{ab}}{k^2},$$

$$V(r)=-\frac{4}{3}\frac{g_0^2}{4\pi r}=-\frac{4}{3}\frac{{\alpha }_s}{r}.$$

$$V(r)=-\frac{4}{3}\frac{{\alpha }_s}{r}+\sigma r,$$

$$P_{\mu \nu }^{ab}(k)=\frac{C(s+k^2)}{t^2+u^2{k}^2+k^4}\left({\delta }_{\mu \nu }-\frac{k_{\mu }{k}_{\nu }}{k^2}\right){\delta }^{ab},$$

$$\displaystyle\begin{aligned} C&=$0.784$,\;&s&=$2.508\text{\,}{\mathrm{GeV}}^{2}$\,,\\ t&=$0.720\text{\,}{\mathrm{GeV}}^{2}$,\;&u&=$0.768\text{\,}\mathrm{GeV}$\,,\end{aligned}$$

$$P_{00}^{ab}\big{(}\vec{k}\,\big{)}=\frac{C\,\big{(}s+\vec{k}^{2}\big{)}}{t^{2}+u^{2}\vec{k}^{2}+\vec{k}^{4}}\,\delta^{ab}\,.$$

$$T_{fi}=\frac{4}{3}\,\frac{g_{0}^{2}}{(2\pi)^{6}}\,\frac{C\,\big{(}s+\vec{k}^{2}\big{)}}{t^{2}+u^{2}\vec{k}^{2}+\vec{k}^{4}}\,.$$

$$\begin{gathered} V(\vec{r})=-\frac{4}{3}\frac{{\alpha }_s}{r}\frac{C}{2\pi i} \\ \times {\int }_{-\infty }^{\infty }\frac{\left(s+k^2\right)\left(e^{ikr}-e^{-ikr}\right)}{t^2+u^2{k}^2+k^4}kdk, \end{gathered}$$

$$k_{m,n}$$

$$\theta$$

$$\begin{gathered} \mathrm{Res}\left[\frac{\left(s+k^2\right)k{e}^{\pm ikr}}{t^2+u^2{k}^2+k^4},k_{m,n}\right]= \\ \frac{1}{2}\frac{\left(s+k_{m,n}^2\right)e^{\pm i{k}_{m,n}r}}{u^2+2k_{m,n}^2}. \end{gathered}$$

Full text

††preprint: LPT-Orsay-17-13

I Introduction

A reliable description of heavy quarkonia states is of great interest for our understanding of nonperturbative aspects of QCD Brambilla:2014jmp and is expected to be important in guiding the search for physics beyond the standard model Love:2008ys . A fortuitous advantage in the study of such states is that, due to the large mass of the heavy quarks, various approximations may be adopted. For example, an expansion in inverse powers of the heavy-quark mass $m$ is performed in potential nonrelativistic QCD (pNRQCD) Brambilla:2004jw , and lattice simulations (especially for bottomonium systems) are applied to effective actions obtained by an expansion in powers of the heavy-quark velocity $v/c$. Similarly, in the relativistic quark model with the quasipotential approach, radiative corrections may be included and treated perturbatively in the case of heavy quarkonia Ebert:2011jc . This possibility of exploring different scales of the problem separately is also helpful in methods more directly based on QCD, such as studies of Dyson-Schwinger and Bethe-Salpeter equations Fischer:2014cfa .

An early but still successful approach to describe heavy quarkonia is given by nonrelativistic potential models, to which relativistic corrections may also be added Radford:2007vd .111Note that these corrections may be computed from lattice data for the Wilson loop Koma:2006si ; Koma:2008zza . The idea is to view confinement as an “a priori” property of QCD, modeling the interquark potential to incorporate some known features of the interaction at both ends of the energy scale, i.e. at small and large distances. The simplest such model, the Cornell — or Coulomb-plus-linear — potential Eichten:1978tg ; Eichten:1979ms ; Eichten:2002qv , is obtained by supplementing the high-energy (perturbative) part of the potential with an explicit confining term. Hence, the resulting expression is a sum of two terms: the first one comes from the quark-antiquark interaction in the one-gluon-exchange (OGE) approximation using a tree-level gluon propagator, and the second one is a linearly rising potential. We have

[TABLE]

where ${\alpha_{s}}$ is the strong coupling constant and $\sigma$ is the string tension. The first term may be associated with scattering of the quark-antiquark pair inside the meson and is analogous to the Coulomb potential in the QED case. The second term corresponds to linear confinement as observed from the strong-coupling expansion of the Wilson loop in lattice gauge theory with static quarks. The Cornell potential provides a spin-independent description of the interquark potential for heavy quarks, with parameters determined by fitting a few known states (see e.g. Chung:2008sm ) or by comparison with lattice simulations. For a recent determination of these parameters, see Ref. Kawanai:2011xb . The numerical procedure for obtaining the mass spectrum for the Cornell potential, as well as for other commonly used potentials is reviewed in detail in Lucha:1991vn .

More generally, the static interquark potential may be defined conveniently in terms of the Wilson loop, or it may be obtained (perturbatively) by taking the nonrelativistic limit in the Bethe-Salpeter equation describing the bound state of two heavy fermions. This yields a Schrödinger equation, to which a linear term is added a posteriori. It would be interesting, nevertheless, to have a better insight about confinement as an emergent property of the interquark interaction induced by the gluon propagator, rather than as a built-in feature. Of course, at the hadronic scale, the full gluon propagator in QCD is very different from the perturbative one and it should contain information about confinement. In order to use this nonperturbative information we propose to substitute the free gluon propagator in the OGE term of the potential, as described above, by a fully nonperturbative one, obtained from lattice simulations. We want to check if this replacement leads to an improved description of the spectra, possibly without the need to include the linearly rising term explicitly. To this end, we use the data generated in studies of the SU(2) gluon propagator in Landau gauge on very large lattices (up to $128^{4}$), reported in Cucchieri:2007md ; Cucchieri:2007rg . We note that lattice data for propagators in the $SU(2)$ and $SU(3)$ cases have essentially the same behavior apart from a global constant Cucchieri:2007zm , which can be fixed by choosing a specific multiplicative renormalization condition, as done in the momentum-subtraction scheme. Of course, to include all QCD effects in the analysis, one should consider a gluon propagator obtained from unquenched SU(3) simulations. On the other hand, such simulations have been done Bowman:2004jm ; Ilgenfritz:2006he ; Kamleh:2007ud ; Silva:2010vx only for rather small physical volumes up to now, and with associated unquenching effects that seem to be modest, at the quantitative rather than qualitative level. Moreover, we are interested in the origin of the linearly confining term of the static interquark potential, which should already show up in the pure-Yang-Mills sector of the theory, which is confining, without the need to include unquenching effects. Thus, we choose to use our SU(2) lattice data Cucchieri:2007md ; Cucchieri:2007rg ; Cucchieri:2007zm , for which data with good accuracy and well controlled finite-volume effects are available.

We organize this paper in the following way. In Section II we review the procedure for obtaining the Coulomb potential in QED as the nonrelativistic limit of $e^{-}e^{+}$ scattering (at tree level) and the analogous calculation in the heavy-quark case. We then follow the same procedure using the lattice gluon propagator to obtain a nonperturbatively corrected OGE potential, i.e. we use directly the fit obtained in Ref. Cucchieri:2011ig and perform the Fourier transform analytically to get the potential. The result is compared to the perturbative (Coulomb-like) potential in Fig. 2. In Section III we describe the numerical method for obtaining the mass spectra associated with a given interquark potential in the nonrelativistic approximation. We also outline our choices for the interquark potentials, the fitting parameters, and the experimental data used for input and comparison. Our results for the spectra and interquark distances are reported in Section IV and our conclusions in Section V. Preliminary versions of our study have been presented in Refs. Serenone:2012yta ; Serenone:2014ota ; Serenone:2015qra . We note again that our aim is to gain a qualitative understanding of the interplay between perturbative and nonperturbative features of the interquark potential. Our approach is similar in spirit to the one in Refs. Gonzalez:2011zc ; Vento:2012wp ; Gonzalez:2012hx , but our conclusions are somewhat different.

II Potential from Lattice Propagator

Let us first review how the Coulomb potential is obtained in the nonrelativistic limit of QED from the application of Feynman rules to the electron-positron system. The scattering-matrix $S_{fi}$, from which the interaction potential may be obtained, is given by

[TABLE]

where $Q$ and $P$ correspond respectively to the final and initial total momentum and $T_{fi}$ is the scattering amplitude. The two tree-level Feynman diagrams contributing to $T_{fi}$ (see Fig. 1) correspond to the $t$ and $s$ channels, respectively coming from scattering with one photon exchange and to annihilation and creation of an $e^{-}e^{+}$ pair. We get

[TABLE]

where

[TABLE]

and

[TABLE]

We follow the notation in bjorken1965relativistic ; Lucha:1991vn ; Lucha:1995zv : $p_{i}$ denotes the momentum of the incoming particles and $q_{i}$ of the outgoing ones. The particles’ initial and final spins are respectively $\sigma_{i}$ and $\tau_{i}$. We represent the photon propagator by a function $P_{\mu\nu}(k)$ of the photon momentum $k$.

We then make the nonrelativistic approximation, i.e. we impose the kinetic energy of the system to be much smaller than its rest energy ($|\vec{p}|\ll m\cong E$). The four-component state vectors become

[TABLE]

and

[TABLE]

In this approximation, $T_{fi}$ can be written as

[TABLE]

To compute the exchange term, we adopt the Dirac representation for the gamma matrices and the center-of-momentum frame, obtaining

[TABLE]

with222We are using the metric $g_{\mu\nu}=\operatorname{diag}(-1,1,1,1)$, which will be more convenient when we consider Wick rotations later.

[TABLE]

For the annihilation term, note that conservation of momentum at the vertices implies that

[TABLE]

For QED, the Feynman-gauge propagator is given by the expression $P_{\mu\nu}(k)=-g_{\mu\nu}/k^{2}$. As seen in Eqs. (6) and (7), the spinors are momentum-independent in the nonrelativistic approximation and, while $t_{\rm exch}$ is proportional to $1/\vec{k}^{2}$, we see that $t_{\rm annihil}$ will be proportional to $1/4m^{2}$. Thus, we can neglect annihilation effects and the scattering amplitude is given by

[TABLE]

The potential can then be obtained as an inverse Fourier transform, which leads to the Coulomb potential

[TABLE]

For QCD, we replace the photon by the gluon and the electron-positron pair by a quark-antiquark pair. The scattering amplitude will continue to be expressed as a sum of the two terms, now given by

[TABLE]

and

[TABLE]

where $c_{(1,2),(i,f)}$ are three-component color vectors and $\lambda^{a}$ are the Gell-Mann matrices.

Let us note that, with respect to the QED case, the terms $t_{\rm exch}$ and $t_{\text{annihil}}$, which contribute to the scattering amplitude $T_{fi}$ in Eq. (3), now have multiplicative (Casimir) factors, coming from the sum over colors. This sum is obtained assuming that the incoming/outgoing quarks and antiquarks have equal probability of being in a given color state and imposing a color-diagonal gluon propagator. Then, these factors are given respectively by

[TABLE]

and

[TABLE]

Therefore, annihilation effects do not contribute, independently of the nonrelativistic approximation. If we now assume a free (i.e. tree-level) gluon propagator

[TABLE]

we obtain a Coulomb-like interquark potential

[TABLE]

Notice that the above potential is non-confining. This could have been expected, since we have performed a purely perturbative calculation, while confinement is a nonperturbative phenomenon. The confinement property can then be obtained by addition of a linear term, as described in Section I, leading to the Cornell, or Coulomb-plus-linear, potential Eichten:1978tg ; Eichten:1979ms ; Eichten:2002qv

[TABLE]

which describes surprisingly well the states of charmonium and bottomonium.

As mentioned in Section I, we substitute the free propagator by a fully nonperturbative one in the OGE term. More precisely, we use the propagator333The energy scale used to convert $s$, $t$, and $u$ from lattice to physical units was set using the value $\sqrt{\sigma}=0.44$ GeV for the string tension.

[TABLE]

obtained from fits of lattice data for a pure $SU(2)$ gauge theory in Landau gauge Cucchieri:2011ig . Note that the above parameters correspond to a value $1/k^{2}$ at 2 $\mathrm{GeV}$. Here we choose to normalize the propagator to $1/k^{2}$ at $k\to\infty$, i.e. we adopt $C=1$.

We now follow the same procedure as in the QED case. From Eq. (9) we notice that, in the nonrelativistic approximation, only the component $P_{00}(0,\vec{k})$ survives in the $t_{\text{exch}}$ term and thus the term $k_{\mu}k_{\nu}/k^{2}$ vanishes [see Eq. (10)]. Lastly, in order to convert the propagator in Eq. (21), which was evaluated in Euclidean space, to Minkowski space, we undo the Wick rotation, taking $\delta_{\mu\nu}\rightarrow-g_{\mu\nu}$. We obtain444Let us recall that the propagator is a gauge-dependent quantity. A gauge-independent potential obtained from the (Coulomb-gauge) propagator is discussed in Popovici:2010fy .

[TABLE]

This leads us to the following scattering amplitude

[TABLE]

The potential is obtained, as was done in the QED case [see Eq. (13)], as a Fourier transform of the scattering amplitude $T_{fi}$. We use spherical coordinates for $\vec{k}$ and set $\vec{r}=r\,\hat{z}$. The angular integration is then trivial, resulting in555For the evaluation of this integral only, we will denote $\left\lvert\vec{k}\right\rvert=k$.

[TABLE]

where $\alpha_{s}=g_{0}^{2}/4\pi$ [see Eq. (19)]. The integral in Eq. (24) can be solved using residue calculations. The four poles in the integrand are symmetrically distributed in the four quadrants of the complex plane. We index these poles in the following way

[TABLE]

where

[TABLE]

The associated contour integral is performed by considering its two terms separately: for the term with $e^{ikr}$ (respectively with $e^{-ikr}$) we close the contour above (respectively below). The residues are given by

[TABLE]

The result is simplified by noticing that $\,k_{1,0}=-k_{0,0}$, $\,k_{0,1}=-k_{1,1}$ and $\,k_{1,1}=k_{0,0}^{*}$. The obtained potential, which we call the lattice-gluon-propagator potential $V_{LGP}$, is then

[TABLE]

where $\Re$ indicates the real part.

We note that the only difference with respect to the perturbative (Coulomb-like) case [see Eq. (19)] is given by the expression within brackets. In order to get a quantitative comparison between the two results, we only need to set the value of the strong coupling constant $\alpha_{s}$ in Eqs. (19) and (28). To this end, we evaluate $\alpha_{s}$ —at the energy scale of the mass of the 1S quarkonia states [respectively $J/\psi$ and $\Upsilon$(1S) in the charmonium and bottomonium cases]— by using the four-loop formula and the $\Lambda_{\mbox{QCD}}$ values in Ref. (Agashe:2014kda, , Section 9). This yields $\alpha_{s}\approx 0.2663$ for the charmonium and $\alpha_{s}\approx 0.1843$ for the bottomonium. The resulting potentials are compared (for the charmonium case) in Fig. 2.

We see that the two curves are clearly different, with $V_{LGP}(r)$ rising above zero at around the hadronic scale (i.e. for $r\approx 1$ fm). For larger distances, the curve drops and it can be observed that the potential $V_{LGP}$ is also non-confining. Thus, since (tree-level) perturbation theory was applied, the property of confinement was lost, even though the used propagator was obtained nonperturbatively. Nevertheless, one may hope to describe the first few bound states of the spectrum solely using $V_{LGP}$. This is done in Section IV. We also consider the addition of a linearly rising term $\sigma\,r$ to the potential in order to model confinement, as done for the Cornell-potential case. In this case, the resulting expression is the lattice-gluon-propagator-plus-linear potential

[TABLE]

Note that the nonrelativistic approximation removes any spin dependence from the interactions. This means that, in our description, states with different spin values will be degenerate.

III Numerical Method

Let us consider a central (nonrelativistic) potential describing the interaction between two particles. Since we are dealing with a two-particle system, we can write the Hamiltonian in terms of relative coordinates and use separation of variables in the resulting partial differential equation to isolate the angular part of the wave function, given by the spherical harmonics. Lastly, we perform the usual substitution of variables in the radial wave function $R(r)=f(r)/r$ to obtain the ordinary differential equation (ODE) for $f(r)$

[TABLE]

where $\mu$ is the reduced mass

[TABLE]

$l$ is the quantum number associated with the angular momentum and $m$ is the mass of the heavy (charm or bottom) quark. We use units such that $c=\hbar=1$. Notice as well the addition of the rest mass of the particles, which will allow us to compare the eigenvalue directly with the mass values given in Ref. Agashe:2014kda .

The above ODE has to be solved with proper boundary-value conditions. The first condition is that $f(0)=0$. This comes from the requirement that $R(0)$ be non-singular. A second condition is that $f(r\to\infty)=0$ and comes from the fact that $R(r)$ is normalized, i.e.

[TABLE]

In the limit of large $r$, the potential is dominated by the linearly rising term and the ODE becomes

[TABLE]

The general solution of this equation is the linear combination of the Airy functions $Ai(\rho)$ and $Bi(\rho)$ abramowitz1964handbook , where $\,\rho=(2\mu\sigma)^{1/3}\,r$. However, the Airy function of the second kind $\operatorname{Bi}(\rho)$ diverges at large $\rho$ and therefore it does not obey the boundary condition at infinity. For $\rho>0$, the Airy function of the first kind can be written as

[TABLE]

where $\operatorname{K}_{\nu}(x)$ is the modified Bessel function of the second kind. One can try the Ansatz $f(\rho)=g(\rho)\operatorname{Ai}(\rho)$ and use the property $\operatorname{K}_{\nu}^{\prime}(x)=\nu\operatorname{K}_{\nu}(x)/x-K_{\nu+1}$ in the ODE in Eq. (30) to obtain a second-order ODE with coefficients in terms of $\operatorname{Ai}(\rho)$ and $\operatorname{K}_{4/3}(x)$. However, by expressing these functions as a power series in $\rho$, one clearly sees that an analytic solution would be challenging, even for the simpler case of the Cornell potential. We therefore seek a numerical solution of the problem.

To this end, we use the so-called shooting method press2007numerical . It consists in picking trial values in a discretized range for the eigenenergies, integrating the ODE for each of these values to obtain the corresponding wave function, and choosing the energies for which the wave function obeys the boundary conditions approximately. We use the backward second-order Runge-Kutta method to integrate the wave function, starting from a maximum value $r_{max}$ for the radial coordinate until the origin, in steps of $dr$ (we adapt the method as presented in Ref. press2007numerical by adopting a negative integration step). We choose $r_{max}$ sufficiently large so that we can use $f(r_{max})=\operatorname{Ai}(r_{max})$ and $f^{\prime}(r_{max})=\operatorname{Ai}^{\prime}(r_{max})$ as initial conditions. In practice, the wave function will not obey the boundary conditions exactly since the proposed energy is unlikely to be an exact eigenenergy. Nevertheless, we may count the number of nodes of the wave function: each time we observe an increase in the number of nodes when compared with the previously proposed energy, the desired eigenenergy will be between the two proposed values. We further refine our method by adapting the bisection method to search for the eigenenergy in this interval, thus allowing the use of a coarse grid without loss of precision.

We first test the above method for the Cornell potential [see Eq. (20)], with parameters fixed to $\sigma=$1\text{,}{\mathrm{GeV}}^{2}, $2\mu=$1\text{\,}\mathrm{GeV} and $4\alpha_{s}/3=1$. These are the values used in Ref. 2014arXiv1411.2023H , which adopts a different approach (the asymptotic iteration method) for solving the problem. We find agreement with their values up to the 4th and in some cases even 5th decimal place. Similarly, Ref. Chung:2008sm uses yet another numerical method to compute the eigenenergies for a different set of parameters, allowing comparison with our results. In this case we find agreement up to the 3rd decimal places. We must consider that, in this comparison, our parameters are close to but not identical to the ones used in Ref. Chung:2008sm , which might explain the slightly worse agreement than in the comparison with Ref. 2014arXiv1411.2023H .

III.1 Fitting Procedure

In general, we consider an expression for the potential with free parameters, to be fitted to a few experimental values. To find the best fit, we set up a grid of values for these parameters. Then, we compute the eigenenergies for each proposed set of parameters and select the one that best describes the observed spectrum. As a criterion for choosing the optimal parameters we consider the minimization of $\chi^{2}$ in the description of a few input values from experiment, i.e. we pick the set of parameters minimizing

[TABLE]

where $\sigma_{i}$ is the experimental error associated with the energy $E_{i,\text{experimental}}$ and the $E_{i}$’s are the eigenenergies computed numerically. In order to establish a confidence level for our parameters, we use the method described in detail in Ref. press2007numerical , which consists in determining the region in parameter space for which $\chi^{2}/d.o.f.$ increases by less than one unit with respect to its minimum value, for each of the parameters separately. In cases for which the obtained confidence level is asymmetric, we adopt the larger value as the error.

The above prescription indirectly allows us to establish confidence limits for the eigenenergies, by the so-called Monte Carlo method press2007numerical . More precisely, for each parameter, we draw $N=1000$ random numbers following a Gaussian distribution, centered at the optimal value of the fitted parameter and with standard deviation given by the symmetrized error, and evaluate the spectrum for each (generated) synthetic set of parameter values. The corresponding set of eigenenergies is then used to estimate the confidence limits for the bound-state masses.

Notice that the procedure described here can be applied to any central potential. In the next section, we perform several calculations using this method, considering two approximately nonrelativistic systems: charmonium and bottomonium. Of course, since bottomonium states are heavier in comparison with their kinetic energy, we expect to obtain better results in this case than for charmonium.

Let us now outline our choices for the interquark potentials, the fitting parameters, and the experimental data used for input and comparison. We start from the $V_{LGP}$ potential [see Eq. (28)] obtained purely from the lattice gluon propagator. In this case, as free parameters in the fits, we take the strong coupling constant $\alpha_{s}$ and the mass $m$ of the heavy quark. A motivation for including $m$ as a free parameter is that quark masses are not observable directly and depend on the renormalization scheme. The ranges of parameter values ($m$ and $\alpha_{s}$) and corresponding step sizes used to find the eigenenergies in the case of the $V_{LGP}$ potential are given in Table 1. We also list the range and step of integration for the radial distance $r$.

Next, we consider the $V_{LGP+L}$ and Cornell potentials, which share the same parameters [see Eqs. (29) and (20)]. Here one has the string tension $\sigma$ as a possible additional parameter. In order to have a fair comparison, our calculations are done with two free parameters for the three potentials separately. Namely, for the $V_{LGP+L}$ and Cornell potentials, we choose to leave $\sigma$ (which is of nonperturbative nature) and $m$ free and to fix $\alpha_{s}$ to its perturbative value (at the appropriate energy scale). Furthermore, we perform a combined (constrained) fit of charmonium and bottomonium results, leaving as free parameters the two heavy quark masses and $\sigma$. This is our preferred fit. The ranges of parameter values (heavy-quark mass and string tension $\sigma$) and corresponding step sizes used to find the eigenenergies in the case of the $V_{LGP+L}$ and Cornell potentials, as well as the calculated (fixed) values of $\alpha_{s}$, are given in Table 2. We also list the range and step of integration for the radial distance $r$.

Regarding the choice of experimental data for bound-state masses, we recall that spin interactions are not considered in our approach. This implies energy values with high degeneracy (in comparison with the experimental data) and we thus average over states with different spin. A possible averaging procedure, used in Ref. Olsson:1994cv , is to take the degeneracy of each state as a weight. The spin-averaged mass of the states with principal quantum number $n$ and in the $X$-wave state ($X=S,P,D\dots$) is then given by

[TABLE]

where $m_{i}(nX)$ is the mass of each of the $N_{l}$ states with the same angular momentum $l$, and $g_{i}$ is the degeneracy of the state. The uncertainty associated with the above average may be estimated by propagation of errors, taking the width of the resonance peak666We recall that bound states are identified by plotting a histogram of number of particles (cross-section) detected in a collision versus the energy of the collision. When a resonance is found, it is associated to a bound state. as the uncertainty in each mass $m_{i}(nX)$. We refer to this averaging procedure as “MAV1”. A second possibility to average over different spins is to imagine that, if the experiments were not very precise, we would not see several narrow nondegenerate states, but broad degenerate ones, i.e. a low-precision experiment would see the peaks merged. We thus take the spin-averaged mass from the midpoint between the state with lowest energy and the one with highest energy. The error is estimated as half of the distance between these two states. We refer to this method as “MAV2”.

The results corresponding to the two averaging procedures described above are reported in Tables 3 and 4 respectively for charmonium and bottomonium states. As experimental data, we choose to include only the states present in the meson summary table of Ref. Agashe:2014kda that are regarded as established particles. Also, we omit charged states from our tables, since quarkonia states must be neutral. As inputs in the fits, we use the states $1S$, $1P$ and $2S$ (of charmonium and bottomonium).

By fitting (independently) the charmonium and bottomonium spectra using MAV1, for the $V_{LGP+L}$ and Cornell potentials, we find very large values of $\chi^{2}/d.o.f.$, varying from $1.9\times 10^{2}/d.o.f.$ (for bottomonium and $V_{LGP+L}$) to $23\times 10^{3}/d.o.f.$ (for charmonium and Cornell potential). Instead, the procedure MAV2 gives acceptable values for $\chi^{2}/d.o.f.$ in the charmonium case, and larger values for bottomonium (still, orders of magnitude smaller than with MAV1). These $\chi^{2}/d.o.f.$ values improve if an unconfirmed state of bottomonium is included [namely, the $\eta_{b}(1S)$]. We choose MAV2 as our preferred method.

IV Results

We now follow the procedure described in Section III and obtain, for a given potential, higher eigenenergies of the spectrum from fits to a few low-lying states. A natural first attempt is to consider the $V_{LGP}$ potential in Eq. (28) for the charmonium and bottomonium spectra. As explained in the previous section, we do this by leaving the strong coupling constant $\alpha_{s}$ and the mass of the heavy quark as free parameters, taking the 1S, 2S and 1P energy states as inputs in the fits. Spin averages are done using the MAV2 method. The corresponding results are presented in Tables 5 and 6. Note that we also show the difference between each evaluated mass and the corresponding experimental value. Long dashes represent states that have not been observed experimentally. As can be seen, although the potential is non-confining (see discussion in Section II), the existence of a few lowest states is qualitatively reproduced in the spectrum. This is in agreement with the study in Gonzalez:2012hx for charmonium states using an equivalent approach.

However, it is clear that the pure OGE potential $V_{LGP}$ is not enough to model the spectrum beyond its lowest states, or even to provide a quantitative description of these states. In fact, the spacings between energy levels are not compatible with the experimental values, both for the higher (estimated) states and for the lower ones used as inputs. Moreover (see Tables 5 and 6), we find that the energy states “saturate” around a maximum value. As for the fit parameters, we obtain $\alpha_{s}=0.95\pm 0.02$ and $m_{c}=$2.064_{-0.009}^{+0.010}\text{,}\mathrm{GeV}, for the charmonium, and $\alpha_{s}=0.513^{+0.009}_{-0.010}$ and $m_{b}=$5.10947^{+0.00016}_{-0.00014}\text{\,}\mathrm{GeV}, for the bottomonium. We note that, while the bottom quark mass $m_{b}$ is not very far from the experimental one (see the third column in Table 7) the charm quark mass $m_{c}$ is almost twice the experimental datum (see again Table 7). At the same time, in both cases, the value obtained for $\alpha_{s}$ is quite far from the perturbatively estimated one (see Table 2).

In the remainder of this section, we thus use the potential $V_{LGP+L}$ [see Eq. (29)], obtained by the addition of a linearly growing term to $V_{LGP}$, as well as the Cornell potential [see Eq. (20)] to generate the spectra, and perform a comparison of the results with the experimental data. The idea is to combine the feature of an improved description of the short-distance behavior of the system, as found above using the $V_{LGP}$ potential, with the imposition of a linear behavior at large distances, which should help in obtaining the higher energy states. We will carry out the spectrum calculation —as explained at the end of Section III— using the same set of parameters for the two potentials. As above, we consider the MAV2 averaging method (see Tables 3 and 4), including the unconfirmed $\eta_{b}(1S)$ state.

The data obtained in the independent fits of charmonium and bottomonium spectra are then used to set up a constrained fit, i.e. one with a common value for the string tension $\sigma$ of the two systems. Notice that, for this constrained fit, we have three free parameters ($m_{c}$, $m_{b}$ and $\sigma$) and six states as inputs in the fit (the states 1S, 2S and 1P of charmonium and bottomonium), resulting in three degrees of freedom. The results of this fit using the $V_{LGP+L}$ and Cornell potentials are shown in Table 7. The corresponding spectra are reported in Tables 8 and 9. A visual representation of the spectra is provided in Figs. 3 and 4. Let us remark that the obtained value for the string tension $\sigma$ in our preferred fit (see Table 7) is rather close to the input value used to set the scale for the lattice gluon propagator, $\sigma\approx 0.194$ GeV2 (see Footnote 3), providing a nice consistency check.

For the charmonium spectrum, we obtain smaller errors and (nevertheless) a slightly better agreement with the spin-averaged experimental values in the $V_{LGP+L}$ case than in the Cornell-potential one (see Table 8). In the bottomonium case the results obtained with the two confining potential are comparable (see Table 9). Also, the central value for the string tension in the $V_{LGP+L}$ is slightly closer to the one used to set the energy scale for the lattice propagator.

The fact that the calculated spectra are very similar in the $V_{LGP+L}$ and Cornell potential cases can be understood if we note that, although the pure-OGE potentials were visibly different (see Fig. 2), the inclusion of the linear term brings the two potentials closer, as shown in Fig. 5.

An advantage of our approach is that we have direct access to the radial wave function $f(r)$. We plot, as an example, the wave functions777The wave functions obtained using our code are not normalized. We interpolate the data and normalize $f(r)$ a posteriori. for the 1S state for both potentials in the charmonium and bottomonium cases in Fig. 6. Thus, we see that the similarity between the two potentials (and the obtained spectra) is present for the wave functions as well. Also, note that the wave function is more extended for the charmonium states, as expected.

This direct access to the wave function can be of interest in other applications, such as effective field theories, for which one needs information on the typical distance between the quarks Brambilla:2014jmp . We estimate this quantity by computing

[TABLE]

Some of these typical distances are presented in Tables 10 and 11.

Finally, we could also estimate decay widths, which are proportional to $|R(0)|^{2}$. Notice, however, that this calculation would require a more strict control of the numerical integration in the region near the origin, since the function $R(r)=f(r)/r$ typically shows a divergence for $r\to 0$. This is beyond the scope of the present work.

V Conclusions

We briefly reviewed the potential-model approach for determining the spectrum of quarkonia and discussed the simplest such approach, the Cornell potential. We then modified the procedure for obtaining the OGE potential, by replacing the free gluon propagator with one obtained using lattice simulations. The resulting $V_{LGP}$ potential is different from the Coulomb-like potential, but is still non-confining. Inspection of Fig. 2 shows that, up to the hadronic scale, the potential rises above zero, with a trend to rise further. This is no longer true for larger values of $r$, for which the potential is damped. In fact, in order to obtain a confining (linear) potential, the gluon propagator should show a strong divergence, of $1/k^{4}$, in the infrared limit, as proven in Ref. West:1982bt . Also, an oscillating behavior — due to the complex poles of the lattice propagator Cucchieri:2011ig — is observed. We solve the associated Schrödinger equation numerically and compare our results with the spin-averaged spectrum in Tables 3 and 4. The spectrum obtained from this potential shows the interesting qualitative feature of approximately reproducing a few low-lying eigenstates. This confirms our expectation that the short-distance behavior of the potential is improved by using the fully nonperturbative gluon propagator instead of the tree-level perturbative one. A quantitative description of the spectrum including higher states is, however, not possible.

We therefore add a linear term to $V_{LGP}$, obtaining the $V_{LGP+L}$ potential in Eq. (29). We then compute the eigenenergies for the $V_{LGP+L}$ and Cornell potentials, both for charmonium and for bottomonium states. The spectra obtained using $V_{LGP+L}$ show a slight improvement over the Cornell potential, but no qualitative differences are observed. In particular, the resulting potentials are rather similar, as seen in Fig. 5.

We were also able to obtain the wave functions for all the states, which allows us to estimate the corresponding interquark distances. Let us note that the wave functions are remarkably similar for the $V_{LGP+L}$ and Cornell potentials (see Fig. 6), even though the potentials are not identical (see Fig. 5). This might suggest that the wave function is somewhat insensitive to details of the potential. In fact, a visual comparison between our wave functions and the one presented in (Kawanai:2013aca, , Fig. 5) (corresponding to a different parametrization of the Cornell potential) shows that they are also essentially identical.

Let us mention that a study using a similar method was carried out in Refs. Gonzalez:2011zc ; Gonzalez:2012hx to propose a potential for heavy-quarkonium states. In that case, the gluon propagator was taken from a study of Schwinger-Dyson equations Aguilar:2008xm . This propagator is in qualitative agreement with the lattice results we use. The main difference with respect to our study is that these authors do not include the linear term in the potential, but consider an additive contribution888Let us recall that a constant term in the interquark potential can also be related to the infrared divergence of the Fourier integral of a “confining” gluon propagator $1/k^{4}$ Lucha:1991vn . to the OGE potential, in such a way that the zero of the proposed potential coincides with the Cornell one. This corresponds to fixing the (infinite) self-energy of the static sources Necco:2003jf , which, however, is not present when considering only the OGE diagram at tree level. The spectrum obtained in Gonzalez:2012hx is in general agreement with the expected values.

We stress again that the aim of this work was to gain a qualitative understanding of the interplay between perturbative and nonperturbative features of the interquark potential. As verified in our study, even though the full nonperturbative gluon propagator was used, the potential is non-confining, i.e. confinement is washed away by the use of the (tree-level) perturbative approximation for the interaction. Nevertheless, the resulting potential (with the addition of a linear term) provides a slightly better description of the spectra, with the same number of fit parameters as the Cornell potential. Our preferred fit is done considering simultaneously the charmonium and bottomonium spectra, leaving as free parameters the two heavy quark masses and the string tension $\sigma$. We remark that leaving a known quantity as a free parameter allows a further check of the model’s consistency. Of course, it would be interesting to check if the inclusion of relativistic corrections, as done in Refs. Radford:2007vd ; Olsson:1994cv , would allow a more accurate description of the spectra.

Acknowledgements.

The authors thank B. Blossier and F. Navarra for useful comments. W.S. thanks the Brazilian funding agencies FAPESP and CNPq for financial support. A.C. and T.M. thank CNPq for partial support.