Calibrating the Na\"ive Cornell Model with NRQCD

Building on results from a previous study on heavy meson spectroscopy, we calibrate the Cornell model employing NRQCD predictions for the lowest-lying bottomonium states up to N$^3$LO, in which the bottom mass is varied within a wide range. We find that the Cornell model mass parameter can be identified, within perturbative uncertainties, with the MSR mass at the scale $R = 1\,$GeV. This identification holds for any value of $\alpha_s$ or the bottom mass, and for all perturbative orders investigated. Furthermore we show that: a) the `string tension' parameter is independent of the bottom mass, and b) the Coulomb strength $\kappa$ of the Cornell model can be related to the QCD strong coupling constant $\alpha_s$ at a characteristic non-relativistic scale. We also show how to remove the $u=1/2$ renormalon of the static QCD potential and sum-up large logs related to the renormalon subtraction by switching to the low-scale, short-distance MSR mass, and using R-evolution. Our R-improved expression for the static potential remains independent of the heavy quark mass value and agrees with lattice QCD results for values of the radius as large as $0.8\,$fm, and with the Cornell model potential at long distances.


Introduction
The discovery of the J/ψ in 1974 [1,2] caused a revolution in hadron spectroscopy. Despite the limited development of Quantum Chromodynamics (QCD) for heavy quarkonium systems at that time, the large mass of the c and b quarks made a description in terms of non-relativistic phenomenological models justified. Among those we count the early works of Eichten [3], Godfrey [4], Stanley [5] or Bhanot [6], which showed that both light-and heavy-meson spectroscopy could be described in an unified and simple framework.
Within the quark model working hypothesis, quarks are assumed to be low-energy spectators, bounded due to flavor-independent gluonic degrees of freedom. Perturbative dynamics dominate at short distances, while at long distances non-perturbative effects become manifest. The short-distance interaction is assumed to be dominated by a single t-channel gluon exchange (s-channel gluons are power suppressed), what leads to a Coulomb interaction proportional to the strong coupling constant α s , in which the electric change is replaced by the first SU (N c ) Casimir C F = (N 2 c − 1)/(2N c ) = 4/3 for N C = 3. The longrange linear confining interaction has been confirmed by lattice QCD calculations [7,8], and can be understood under the flux-tube picture where the confinement arises from the chromoelectric potential. Recent investigations also show that it can be understood purely in terms of perturbation theory [9,10]. The simplest quark model for heavy quarkonium is based on the so-called Cornell potential, which assumes linear plus Coulomb-like interactions : V Cor (r) = σ r − κ r . (1.1) Here κ and σ are purely phenomenological constants of the model with, at first glance, no direct connection with QCD parameters, which need to be defined within a given scheme and in general depend on the renormalization scale µ. In any case, to ease a comparison with the O(α s ) expression for the QCD static potential we define κ ≡ C F α Cornell s . However, even though the quark model approach to study heavy quarkonium has some common features with a full fledged QCD approach, and therefore one expects to be connected to QCD basic principles, it is an "ad hoc" construction whose foundations do not need to be completely solid from a theoretical point of view. Indeed, to our knowledge, no work has ever studied the connection between quark potential models and non-relativistic QCD (NRQCD). Since these models are still widely used by a relatively large scientific community, establishing a connection between them and QCD appears certainly warranted.
The small velocity of the charm and bottom quarks in QQ bound states enables the use of non-relativistic effective theories within QCD to study heavy quarkonia. The use of non-relativistic approaches implies that heavy quarkonia bound states can be organized in terms of the non-relativistic quantum numbers (j, , s) and the radial excitation number n, while the hyperfine splittings are power corrections, starting at order O(m −2 Q ). For heavy quarkonia we can distinguish three well defined scales : the heavy quark mass m Q acting as the hard scale, the soft scale determined by the relative momentum of the QQ system (p ∼ m Q v with v c) and the ultrasoft scale marked by the average kinetic energy of the heavy quarkonia (E ∼ m Q v 2 ). These scales have a well-defined hierarchy (m Q p E), which allows for significant simplifications.
NRQCD [11] is obtained from QCD by integrating out the heavy quark mass m Q , and one can exploit the fact that m Q Λ QCD to perform a perturbative matching. This implies that NRQCD inherits all the light degrees of freedom from QCD. The NRQCD Lagrangian is thus expanded as a series in 1/m Q powers, factorizing the terms which contribute at the hard scale as Wilson coefficients. Even though the hard scale lays in the perturbative domain, the soft and ultrasoft scales are probed e.g. when building up QQ systems, which in some cases could make computations in terms of partonic degrees of freedom unreliable. Additionally, there is still a fundamental problem in NRQCD : it does not distinguish soft and ultrasoft scales, which complicates the power-counting.
A solution to this problem arrived with the construction of EFTs such as velocity NRQCD (vNRQCD) [12] and potential NRQCD (pNRQCD) [13,14], which describe the interactions of a non-relativistic system with ultrasoft gluons, organizing the perturbative expansions in α s and the velocity of heavy quarks systematically. Such EFTs only include degrees of freedom relevant for QQ systems near threshold, while the rest of degrees of freedom are integrated out. Indeed, pNRQCD is obtained from NRQCD by integrating out soft gluons and quarks, and potential gluons, under the assumption p Λ QCD .
The specific treatment of the remaining degrees of freedom will depend on the relation between E and the hadronic scale Λ QCD . At long distances QCD becomes strongly coupled and hadronic degrees of freedom emerge. The QQ system is assumed to be dominated by perturbative physics and non-perturbative corrections are taken into account by local or non-local condensates. Under these assumptions, the static potential of QCD for colorsinglet states can be computed, which is by itself an interesting field of research. On the one hand, it is fundamental to study the energy spectrum and, on the other, it is a good tool to explore weak and strong coupling regimes and analyze phenomena such as confinement. In fact, the static potential can be directly compared with lattice QCD simulations [15].
In this article we attempt to calibrate the simplest realization of the Cornell model against NRQCD. To that end we compare observables that can be reliably predicted both in the theory and the model, the mass of the lowest-lying QQ bound states, varying the quark mass and the strong coupling constant. This exercise is inspired by a similar analysis carried out in Ref. [16] for the top quark parameter embedded in the parton-shower Monte Carlo Pythia [17,18]. As a byproduct, we show that the Cornell potential agrees for large values of r with the QCD static potential once the latter is expressed in terms of the MSR mass and improved with all-order resummation of large renormalon-related logs via R-evolution. Our R-improved static potential also compares nicely with lattice QCD simulations from Refs. [19,20].
This article is organized as follows : In Sec. 2 we discuss the resolution of the twobody problem for the Cornell potential. A fit for the parameters of the Cornell model to experimental bottomonium and charmonium data is carried out in Sec. 3. Section 4 introduces the concepts of the MSR mass and R-evolution. A discussion of the QCD static potential, and how its leading renormalon is canceled is presented in Sec. 5. Sec. 6 contains a comparison to lattice results. NRQCD analytic results for the mass of QQ bound states are compiled in Sec. 7. The calibration procedure is explained in Sec. 8, while the results are presented in Sec. 9. Sec. 10 contains our conclusions. In Appendix A we discuss the numerical resolution of the Schrödinger equation for the Cornell potential using the Numerov method.
2 The Cornell Potential for the QQ system As already mentioned previously, we want to obtain the spectrum of the Cornell potential and compare it with NRQCD predictions to extract relations among the Cornell model parameters and QCD fundamental constants. Therefore, we need to solve a two-body problem with a central potential.
The relative wave function of the quark-antiquark pair ψ can be factorized in an angular part, expressed in terms of spherical harmonics, and the radial wave function R n (r) : Here is an integer non-negative number which represents the orbital angular momentum, and n is a natural number accounting for radial excitation. It simply accounts for the infinitely many (but denumerable) bound states that can be found for a given value of . It should not be confused with the principal quantum number n p = n + > 0, which bounds the possible values of the orbital angular momentum 0 ≤ ≤ n p − 1. The Schrödinger equation can be simplified if written in terms of the reduced wave function u n = r R n (r), yielding an ordinary differential equation for u n (r), where V Cor is given in Eq. (1.1) and m Q is the mass of the heavy quark in the Cornell model. Note that we are using c = = 1 units. Here E n is the binding energy 1 , and the mass of the quarkonium bound state is given by M n = 2 m Q + E n . In Fig. 2 we show the effective Cornell potential for n = 1, 2 and = 0, 1, for charmonium and bottomonium, using for the parameters the values obtained in Sec. 3, when fitting to experimental data. The figure also shows the approximations at large and small r, as well as the binding energies and classical radii.
The wave function has to be normalized to one, implying that the reduced radial wave function should be square integrable on the positive real axis : where we have used dΩ Y m (Ω)Y m (Ω) = δ m,m δ , , the orthogonality property of the spherical harmonics. This condition is satisfied if at large distances u(r) falls off sufficiently rapidly. The solution of the above equation requires the proper determination of two boundary conditions. On the one hand, for small r Eq. (2.2) becomes u nl (r) + ( + 1)u n (r)/r 2 = 0, which implies a power-like behavior of the reduced wave function with exponent + 1, For large r, the Coulomb and centrifugal terms are largely suppressed, and we get the following asymptotic Schrödinger equation for the radial wave function : which is independent of and can be analytically solved in terms of Airy function The energy quantization of Eq. (2.6) is achieved simply imposing the boundary condition at small r for = 0, that is, that the reduced wave function vanishes at r = 0. The Airy function has an infinite denumerable number of zeros which happen to be all negative. These can be computed numerically and we denote them by a n such that Ai(a n ) = 0 and |a n | > |a n | if n > n . With this at hand we find for the Airy potential self-energies : Therefore the exact solution of the Airy Schrödinger equation takes a very simple form : u Airy n (r) = N Airy n Ai a n + r (m Q σ) 1/3 . (2.9) Indeed Eq. (2.9) behaves linearly near r = 0, and the n-th state crosses the horizontal axis exactly n − 1 times, as expected. This behavior can be seen in Fig. 1   From the asymptotic behavior of the Airy functions we can obtain a simpler boundary condition at large r : . (2.10) For the Cornell potential in Eq. (1.1) it is not possible to find analytic expressions for the energy eigenvalues and eigenfunctions, 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 [21,22] (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 [23], 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.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. Furthermore, 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 interactions, breaking the present degeneracy [24]. Such contributions arise from the s-channel gluon exchange and the leading relativistic corrections of the t-channel gluon and confinement interactions, giving the following terms [4,25] where S = S 1 + S 2 is the total spin, L the relative orbital momentum and S 12 the tensor operator of the QQ bound state, defined as withr = r/r. Analytical expressions can be found for the spin-spin, spin-orbital and tensor operators, with (j, , s) the quantum numbers of the QQ state. The specific values for the quantum numbers considered in this article are given in Tab. 1. Given the assumed large mass of the heavy quark, the above terms are expected to be small. Therefore, their contribution to the bottomonium mass will be incorporated to the model using first-order perturbation theory. For consistency, 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. (2.11). To that end we use the angular decomposition in Eq. (2.1) 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 2 Here the LS term from confinement is obtained by the exchange of a scalar particle. momentum algebra [ see Eq. (2.13) ] and the radial integral becomes a one-dimensional matrix element of the reduced wave function u n (r), given that the Jacobian factor r 2 cancels when writing |R n (r)| 2 = |u n (r)| 2 /r 2 as in Eq. (2.4).

Fitting the Cornell Model to Experimental data
In this section we present the results from χ 2 fits of the 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. [26]. The minimization of the χ 2 is carried out using the Fortran 77 package MINUIT [27]. 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. (3.2) ] it does not estimate the covariance matrix correctly. To solve this problem we compute a three-dimensional grid of 5 3 elements centered around the minimum found by MINUIT, which is then adjusted to a 3D quadratic polynomial by a linear regression. We find that the χ 2 function clearly behaves in Gaussian way near the minimum, which permits an estimate of the covariance matrix in the Hessian approximation, that is, from the quadratic terms of our regression.

Bottomonium Fits
We take the experimental values from the PDG [28], which are collected in Table 3 of Ref. [26]. Since experimental data is extremely precise, we find χ 2 min /d.o.f. = 258, much larger than unity. As a penalty for this deficiency we rescale the uncertainties on the parameters that come out from the fit by the square root of the reduced χ 2 at its minimum. With this procedure we get : The quoted uncertainties correspond to 1 standard deviation for each parameter (i.e. it corresponds to a 68 % confidence level in each of them), as obtained from the intersection of the χ 2 function with the horizontal hyperplane χ 2 min + 1. The 68 % confidence level in the space spanned by these three parameters is obtained by combining the correlation matrix in Eq. (3.2) with the uncertainties given in Eq. (3.1) rescaled by the factor 1.878. It corresponds to the intersection of the χ 2 function with the hyperplane χ 2 min + 3.53. However, the uncertainties as given in Eq. (3.1) together with the correlation matrix are used to compute the incertitude of any function of these parameters (e.g. the masses of the bound states).
Eq. (3.2) shows that indeed there is a very strong correlation between the three parameters, very close to 100 % (anti-)correlation. We find a very strong negative correlation between m Cornell where the order of columns and rows is the same as in Eq. (3.1). These results can be easily interpreted : an increase of either m Cornell Q or σ make the mass of all bound states larger, while one has to decrease α Cornell s produce the same effect. In Fig. 3(a) we show the reduced wave functions, eigenfunctions of the static Cornell potential, for n = 1, 2 and = 0, 1, values that correspond to the states included in the fit. For a given n value, the peak is shifted to the right for larger values of . The exponential falloff starts at larger r values for higher values of n, and is roughly independent. As expected, for n = 2 the wave function crosses the x axis once. In all cases the wave function is strongly suppressed already at r = 1.5 fm, therefore justifying our choice r max = 4 fm. Taking N = 5000 corresponds to a numerical uncertainty in the resolution of the Schrödinger equation of 0.08 % and 0.0004 % for n = 0, = 0, 1 states, respectively, and 0.034 % and 0.0002 % for n = 1, = 0, 1. This choice of N inflicts an uncertainty on the extraction m Cornell b , σ, and α Cornell s from a fit to data of 0.013 %, 0.22 % and 0.017 %, respectively. 3 These are a factor of 28, 24 and 237 smaller than the fit uncertainties, and therefore negligible. The numerical uncertainty associated to r max = 4 fm is even smaller. As an additional check on the accuracy of our numerical method, we compute the = 0 states with α Cornell s = 0, for which we know the exact solution, shown in Eq. (2.8). The numeric solutions reproduce the exact result with an accuracy of 2 × 10 −6 %. At this point it is also worth mentioning that neglecting the Coulomb-like interaction (that is, setting α Cornell s = 0) overshots the exact result of the static energies by 5.5 % and 3 % for = 0 and n = 1, 2, respectively, and 0.6 % and 0.4 % for = 1 and n = 1, 2. Likewise, setting σ = 0 one can use the known solution for the hydrogen atom bound states, which can be obtained from Eq. (7.1) truncated at O(α 2 s ), which undershoots the exact results by 2.4 % and 6.1 %, for = 0 and n = 1, 2, respectively, and 5 % and 8.2 % for = 1 and n = 1, 2. It might seem that either of these approximations is good, but it only appears so because most of the contribution comes from the quark mass. If we compare binding energies then the approximations become much worse, rising up to 600 % for the Coulomb limit and 1400 % for the Airy limit. Therefore none of the two terms in the Cornell potential can be treated as a perturbation with respect to the other for bottomonium.
We have checked the robustness of our fits by removing one, two, or three points from our dataset. It tuns out that if one keeps the two states with n p = 1 and the highest mass state Υ(3 3 S 1 ), the largest variation from the default is always smaller than 10 MeV for the Cornell mass parameter, and below 2 % for either σ or α Cornell s .
In Fig. 4 we compare the predictions of the Cornell model for the 8 states that have been used in the fit, employing the best-fit values and propagating the fit uncertainties into the masses through the covariance matrix. We include experimental data with n p = 3, which have not been included in the fit, but serve as an illustration for the limitations of this simple model. States with n p > 2 can only bee described if some sort of "string breaking" is implemented [29,30]. The N 3 LO QCD predictions of Ref. [26] are also shown for illustration. Table 2. Experimental masses of the charmonium states, up to n = 2. Data from Ref. [28] Experimental data for charmonium has been taken from the PDG [28], and for the reader's convenience is collected in Tab. 2. Given the astonishing precision of some charmonium mass measurements, the reduced χ 2 at the minimum is in this case even larger than for bottomonium. We find χ 2 min /d.o.f. = 1.3 × 10 5 and, consequently, apply the same (d) Figure 2. Effective static Cornell potential for = 0 (left panels) and = 1 (right panels), for bottom (upper panels) and charm (lower panels), as defined in Eq. (2.3). The first two energy eigenvalues are shown as horizontal lines in red (ground state) and green (first excited state). The corresponding classical radius for these energies are signaled with vertical, black, dashed lines. We also show the asymptotic behavior for large r in dashed pink, and the small r dominant terms in dashed orange (Coulombic for left plot and centrifuge for right plot), and for the right plots the Coulombic plus centrifuge terms added together, in dashed gray.

Charmonium Fits
penalty as was done in Sec. 3.1. We obtain : Two comments are in order : a) The value of α Cornell s is, as expected, larger than for bottomonium. To the extent that this model parameter can be related to the QCD strong coupling at some characteristic energy, one expects the typical scale for charmonium smaller than for bottomonium, which translates into a larger α s value for the former. b) The uncertainties for the quark mass and α s are larger than for bottomonium, while remain the same for σ.   The larger errors can be explained by the huge penalty factor applied here, namely 363, which is 23 times larger than for bottomonium. Long-distance interactions matter more in charmonium, since it is a more extended system in which softer gluons are more often exchanged. Accordingly, the confining parameter parameter is fixed less ambiguously than the rest. The correlation among the three parameters is even stronger than for bottomonium, and follows the exact same pattern : Similar to what we found in bottomonium, the static binding energies are not well reproduced if any of the two terms in the potential is set to zero, but neglecting the Coulomb term is a much better approximation. Setting σ = 0 produces binding energies which are incorrect by up to 2300 %, but setting α s = 0 overshoots the exact answer by 60 % at worst. This confirms the expectation that charmonium bound states are more afflicted by long-distance effects, parametrized by σ in the Cornell model, and also explains why this parameter is determined more accurately in Eq. (3.3). In Fig. 3(b) we show the charmonium reduced wave functions. They have identical shapes to those of bottomonium, but as anticipated from dimensional arguments, the wave functions extend to larger distances as compared to bottomonium states, and have accordingly lower (shallower) peaks (valleys). However, it is still justified to use r max = 4 fm as at 2 fm the wave function is already negligibly small. We have checked that using r max = 6 fm produces a shift in the 4-th or 5-th decimal place in any of the fit parameters shown in Eq. (3.3). Since we will not calibrate the Cornell model for charmonium, we do not study the optimal choice of N , the number of steps in the Numerov method. We simply take 200 000 for which the associated numerical errors are several orders of magnitude smaller than the fit uncertainties.
We have performed a similar robustness check in our charmonium fits, and we conclude that the stability is not as good as for bottomonium. We have only removed a single point from our datasets, and always keeping the two lowest-mass states and the highest-mass state, and we find that the impact on the Cornell mass parameter can be as large as 70 MeV, while in the other two parameters can reach deviations of 6 % and 17 % for σ and α Cornell s , respectively. Since we are not going to calibrate the charm mass, these finding are only a minor concern, which again seem to signal that the non-relativistic approximation is worse for charmonium.
In Fig. 5 we again confront the predictions of the Cornell model after the fit with the 8 masses that enter our χ 2 function. In this case we refrain from showing results for n p = 3, since the DD threshold is quite low in charmonium. Similarly, we only show QCD predictions (taken from Ref. [26]) for n p = 1, since higher n p states cannot be reliably predicted in perturbation theory. It is worth noting that the Cornell model prediction for the J/Ψ is extremely precise, with an error bar which seems unnaturally small given the uncertainties quoted in Eq. (3.3). The explanation for this apparent contradiction is the large cancellations that happen between the different error sources due to the very strong correlations among the three parameters. The J/Ψ, having an experimental uncertainty 40 times smaller than the next most precise measurement, overly drives the fit and renders the correlation such that this particular mass is predicted with an uncertainty roughly 40 times smaller than for the rest of states. Such a situation does not happen for bottomonium, therefore all states are predicted with similar incertitudes.  Fig. 4 for charmonium bound states. Here QCD has been fitted to n p = 1 experimental data, for which non-perturbative effects are quite small, while the Cornell model has been fitted to states with n p ≤ 2.

The MSR mass and R-evolution
The MSR scheme (see Refs. [31][32][33] for a review), can be seen as the natural extension of the MS mass for renormalization scales below the heavy quark mass. It relies on the fact that the renormalon ambiguity in the MS-pole relation is independent of the value of the mass. The MSR mass is then directly defined from this relation, and depends on an infrared scale R. In addition, the subtraction series relating the MSR and pole masses can be expressed in terms of α s (µ), with µ the MS renormalization scale. This is essential to cancel the renormalon in other series when the pole mass expressed in terms of the MSR mass : 4 Exploiting the µ independence of the MSR mass, the d n,k coefficients can be expressed as a function of d k,0 through the following recursion relation : with β i the coefficients of the QCD beta function, defined as (4. 3) The renormalization group equation of the MSR mass is given by : with γ R n the R-anomalous dimension coefficients [33], which are renormalon free. The γ R n can be easily calculated from the first line of Eq. (4.1) using that the pole mass is R-independent 5) and the renormalon cancels between the first term and the sum. The solution of the RGE in Eq. (4.4) sums up powers of log(R 1 /R 2 ) to all orders in perturbation theory : The MSR mass can be easily matched to the MS mass at the scale R = m Q (m Q ), and then run down from there to any value of R < m Q .
The MS and MSR bottom masses (as well as the QCD static potential and the bottomonium masses) receive corrections from the finite mass of the charm quark. Here and in what follows we assume the charm quark is massless, as its mass plays no role in the calibration of the Cornell model but might add complications when scanning over the bottom quark mass over a large range. We close this section noting that if we compute the MSR mass from the MS employing the world average values m b = 4.18 ± 0.03 GeV and α

The QCD Static Potential
As in the Cornell model, to compute the energy levels of QQ states in NRQCD at leading order in 1/m Q one needs the static QCD potential. It is defined as the potential between two static quarks, that is, the color-neutral interaction between two infinitely heavy colortriplet states. It is well defined up to O(α 4 s ), where the static approximation breaks down and the potential becomes time-dependent. This feature manifests itself in dimensional regularization as an 1/ pole, that once regulated leaves behind a dependence on the "ultrasoft" scale µ us . The pole and µ us logarithms cancel in observables such as quarkonium masses when adding ultrasoft effects from wave-function renormalization. Nevertheless, to define the static potential at four loops we simply subtract the 1/ pole. In position space the perturbative contribution to the static QCD potential can be written in a compact form as follows : 5 The coefficients a i,0 are known to four loops [34][35][36][37][38][39][40][41][42][43] and a i,j>0 can be derived from the former requiring that V QCD does not depend on µ, with a recursion relation identical to Eq. (4.2) : For charmonium (bottomonium) one has n = 3(4) and the coefficients a i,0 take the following numerical values : The V us QCD (r) term in Eq. (5.1) depends on the ultrasoft factorization scale µ us and takes the following form : Following e.g. Ref. [44], for the ultrasoft factorization scale we will take the expression that takes into account the power counting of pNRQCD [13,14]. In this article we do not perform any resummation of large ultrasoft logarithms besides those that can be absorbed in the running of α s . This ultrasoft resummation has been performed at N 3 LL in Ref. [45]. In Refs. [9,10] it is shown from renormalon dominance arguments and in the framework of the operator product expansion of pNRQCD, that perturbation theory alone should be capable of describing both the Coulomb and linear behavior of the static potential, and that non-perturbative corrections start at O(Λ 3 QCD r 2 ). In the following we shall confirm that claim numerically using the MSR scheme.
The static QCD potential suffers from a factorial divergent growth at large orders, also known as a u = 1/2 renormalon, which translates into an O(Λ QCD ) ambiguity. This ambiguity happens to be r-independent, and its nature depends only on the coefficients of the QCD beta function. In Figs. 7(a) and 7(c) this bad perturbative behavior manifests itself as a (roughly constant) vertical shift of the potential in the region between 0.05 fm and 0.2 fm every time a new order is included. 6 Furthermore, none of the orders makes the QCD static potential close to the Cornell model. It is well understood [46][47][48] that the ambiguity exactly matches that of the pole mass except for a factor of −2, such that the static energy E stat (r, µ) = 2 m pole Q + V QCD (r, µ) is renormalon free. Since the pole mass ambiguity is independent of the quark mass itself, the cancellation happens irrespectively of the specific numeric value for the mass. Therefore one only needs to re-write the pole mass in terms of a short-distance scheme to make the cancellation manifest. For the cancellation to take place one also needs to express the perturbative series δm SD Q that relates the pole mass with a short-distance mass m SD Q in terms of α s (µ), as done in the second line of Eq. (4.1). There are powers of log(µ/R) in δm SD Q that may become large if µ and R are very different. Since one has to choose µ such that log(r µ e γ E ) ∼ O(1), that is µ should depend on r, renormalization schemes with a fixed value of R such as the MS, are disfavored. Following Ref. [26] we use the MSR mass and choose µ = R to simultaneously minimize logs in the potential and in δm MSR Q . Since the canonical choice µ = 1/r quickly dives into non-perturbative values, we freeze it to 1 GeV once it reaches this value. To avoid a kink in the potential we smoothly convert the 1/r behavior into a constant employing a transition function between r = 0.08 fm and 0.2 fm. This function is composed of two quadratics smoothly connected with each other and with 1/r and 1 GeV at the junction points, as has been employed for instance in Ref. [49] in the context of event shapes. A graphical representation of this piecewise form, which will be referred to as "profile function", can be seen in Fig. 6, and has the following analytical form where we use r m = (r 0 + r 1 )/2 and µ c is the constant renormalization scale at long distances, taken as 1 GeV. The six coefficients c i in f t1 (r) and f t2 (r) are obtained by imposing continuity of the functions and their derivatives at the transition points r 0 , r m and r 1 .
Consequently, we obtain For our specific choice r 0 = 0.08 fm and r 1 = 0.2 fm, the coefficients c i take the following numerical values :  [50][51][52] in the renormalon-subtracted scheme, but with a somewhat more abrupt transition between the canonical and frozen regimes.
Our next goal is to define a short-distance potential which is renormalon free and independent of the quark mass. This is a sensible criterion, since the static potential is defined in the limit of infinitely heavy quarks. It will however depend on the scale R 0 at which the renormalon is subtracted. The concept of R-evolution comes in very handy to accomplish our task. Some basic algebra brings the static energy into a convenient form : Which allows us to define a renormalon-free potential at a given scale R 0 . To avoid the occurrence large logs of R 0 /µ in δm MSR Q (R 0 , µ) we express it in terms of δm MSR Q (R, µ) and use R-evolution to sum up large logs : The term ∆ MSR , defined in Eq. (4.6), is mass independent and sums up large logs of R/R 0 to all orders. Therefore we can write the following R-improved expression for the MSR-scheme Static QCD Potential : The R-improved static potential is similar to the Renormalon Subtracted scheme used in Ref. [45], based on [53], and to the analysis of Ref. [9] based on renormalon dominance.  1). Other values should simply shift the potential vertically, but for this specific choice we observe that the static MSR potential converges very nicely towards the Cornell model for moderate values of r. When more perturbative orders are added, the agreement becomes better over larger distances. On the other hand, for high values of r, since the renormalization scale is frozen, log(rµ) becomes large, which makes perturbation theory unreliable. At small distances all orders agree very well because α s becomes very small, but disagree with the Cornell model. So we can conclude that the Cornell model and QCD agree for moderate values of r, but disagree in the ultraviolet, as the model does not incorporate logarithmic modifications due to the running of α s . A legitimate question is then if this difference in the UV can be absorbed in the definition of the quark mass. We shall answer this question in Sec. 9 of this article.
The ultrasoft potential in Eq. (5.4) is only a small contribution of the total MSR R-improved potential. At N 3 LO we find its weight is only 0.6 % at short distances, quickly becoming below the per-mil for r 1 fm. To fully asses the impact of this term one would need to do a thorough study of perturbative uncertainties through scale variation, which is beyond the scope of this article.

Comparison to Lattice QCD
In this section we perform a comparison to lattice QCD results. We start by comparing the Cornell model parameters as obtained in fits to data to specific lattice studies that determine them. Finally we compare our R-improved QCD static potential to lattice simulations. Our result for σ in Eq. 3) is worse, and this fact will be discussed next. When it comes to the α Cornell s parameter our bottomonium (charmonium) result is roughly a factor of 2 larger than what is found in these lattice analyses. Given that the static potential at short distances can be described in perturbation theory, we know that loop corrections modify the short-distance 1/r behavior. Therefore fitting a Coulomb-like function to lattice output might be meaningless, and the discrepancy should be of little concern.
Lattice simulations for the static potential use three dynamical quarks, and therefore their results should be interpreted as the charmonium potential. Hence, it is confusing that we find better agreement for the linear confining parameter in the fits to the bottomonium spectrum. However, as pointed out in Ref. [56], the charm quark is effectively decoupled for the low-lying bottomonium states. Therefore the same static potential should be used for charmonium and bottomonium, as long as some (small) charm mass corrections are included in the latter. Therefore one could think that also the same Cornell potential should be used for charmonium and bottomonium systems, and given that the non-relativistic approximation is more accurate for the latter, the right comparison might be between lattice and Eq. (3.1). Let us emphasize again that our Cornell model parameters are not obtained from a direct comparison to the static QCD potential, but from fits to experimental data on quarkonium bound states. At the sight of Figs. 7(b) and 7(d) it is clear that a smaller value of α Cornell s would improve the agreement with V MSR static . We turn now our attention to lattice computations of the static potential. On the lattice, the static approximation also breaks down, and this manifests as a linearly divergent term. This divergence is removed by additive renormalization, as studied in Ref. [57]. We will use the results of Ref [19] for a lattice spacing a = 0.04 fm. These cover the range 0.039 fm ≤ r ≤ 0.84 fm, and have an average relative precision of 2.6 %. This precision seems to be roughly proportional to 1/r 2 , but shows a very pronounced peak around r = 0.25 fm. We complement this dataset with results from [20] with lattice spacing a = 0.025 fm, which cover values of the radius as small as 0.024 fm. Since uncertainties in [20] are larger for higher values of r, we only consider data with r ≤ 0.25 fm. Moreover, in the range 0.024 fm ≤ r ≤ 0.25 fm [20] has more density of points than [19]. Our complete dataset is plotted in Fig. 8 as black dots with error bars. The HotQCD collaboration has results for larger lattice spacing, with predictions for the potential covering values for the radius as large as 1.59 fm. However these do not extend as much into the small-r region where perturbation theory dominates, and therefore we do not take them into account in this simple comparison.
Since the potential is arbitrary up to an additive constant (that is, a vertical offset), which in its R-improved version can be parametrized by the subtraction scale R 0 [ see Eq. (5.11) ], and otherwise dependent only on α s , we perform a two-parameter fit of our R-improved static QCD potential with the scale setting shown in Eq. (5.6) to the lattice data. Since a detailed analysis is beyond the scope of this article, we do not include resummation of large ultrasoft logarithms and do not attempt to estimate perturbative uncertainties, simply using µ = R. Furthermore, we do not study the dataset selection dependence and simply include all lattice points in our χ 2 function. Finally, since the correlation matrix is currently unknown, we assume individual lattice measurements are statistically independent. Within this approximation we find : The value of R 0 is in remarkable agreement to the value used to compare with the Cornell model potential. To relate α s with n f = 3 and 5 dynamical flavors, we use the five-loop QCD beta function to perform the running, and the five-loop threshold matching relation with the values m b = 4.2 GeV, m c = 1.3 GeV. We refrain from showing fit uncertainties as they do not reflect the actual accuracy one can reach by this procedure. The resulting QCD potential which uses the best-fit value for α s and the vertical offset is shown as a red line in Fig. 8. This plot also shows as a blue line the R-improved static QCD potential with the same values for the parameters, but using a fully canonical profile µ = 1/r. Whereas our profile functions seems to capture some of the infrared physics, completely agreeing with lattice QCD results up to distances of approximately 1 fm, the canonical profile behaves in an unphysical manner for r 0.2 fm. We believe our findings can have an impact on the precise determination of the strong coupling constant from fits to lattice simulations on the QCD static potential by enlarging the range of data than can be utilized in constructing the χ 2 . Such determinations have been carried out in Refs. [58,59]. To the best of our knowledge, the maximum value of r considered in those fits is r max 0.23 fm.  Fig. 6. We determine α s and R 0 fitting our pNRQCD theoretical expression to the lattice simulations.

Analytic NRQCD Formula for Bottomonium with Massless Quarks
The complete bottomonium spectrum up to n p = 2 for arbitrary m b ≡ m b (m b ) will be constructed using Non-Relativistic Quantum Chromodynamics (NRQCD) up to N 3 LO [60], which will be briefly explained below (further details can be found e.g. in Ref. [26]). In the pole mass scheme [40,60,61], the energy of a non-relativistic QQ bound state, characterized by the (n p , j, , s) quantum numbers and with n massless active flavors, can be written as : 7 Again we neglect the corrections coming from a finite charm quark mass, which start at O(ε 2 ). We also work in the n scheme, see Ref. [26] for more details.
with H n the harmonic number. In Eq. (7.1) ε acts as a bookkeeping parameter that properly implements the so called Υ-expansion [62,63]. The c i,0 coefficients have been computed up to i = 3 [64,65], while the c i,j>0 coefficients can be directly obtained from the latter c i,0 imposing µ independence of the quarkonia mass. We denote the sum in Eq. (7.2) truncated to O(ε n+1 ) as the N n LO result. This formula does not include the resummation of large ultrasoft logarithms. Recently this summation has been carried to N 3 LL precision for P-wave states in Ref. [52]. The u = 1/2 renormalon in the QCD static potential discussed in Sec. 5 is inherited by the perturbative expansion of the quarkonia mass in Eq. (7.1). Exactly as it happened for the potential, such renormalon is canceled by expressing the pole mass in terms of a shortdistance mass. In the past the MS mass has been employed, although the natural scenarios for this scheme are processes where the involved energy scale is much larger than the heavy quark mass. Even though the QCD static potential starts off at O(α s ), the first correction to the quarkonia masses is O(α 2 s ). 8 When switching the pole mass to a short-distance scheme in Eq. (7.1) the first correction from this conversion will be O(α s ), which appears to follow a different pattern. It is at this point when the use of the ε parameter becomes relevant. Since the renormalon cancellation happens already in the static potential, it is crucial to treat the short-distance mass subtractions in Eq. (7.1) exactly as in the static energy, and therefore terms proportional to α n s in the subtraction are regarded as O(ε n ) in the Υ-expansion counting.
Heavy quarkonium mases probe energy scales below m Q , and therefore the relativistic logs showing up in the MS-pole relation series δ SD , namely log(µ/m Q ), may become large if µ is chosen to minimize the non-relativistic logarithm that appears in Eq. (7.2). Therefore it would be better to simultaneously minimize logarithms appearing in δm SD Q and in L n , whose argument is the ratio of a non-relativistic scale and µ. This is in full analogy with the static potential analysis, where r is being replaced by the Bohr radius in bound states masses. Therefore, a low-scale short-distance mass is advisable. Following the results of Sec. 5 and the analysis in Ref. [26] we will employ the MSR mass [33] for our analysis.

Bottomonium masses with a floating bottom mass
The ultimate aim of this article is to calibrate the Cornell model constants (especially the Cornell mass) in terms of the QCD fundamental parameters α s and m b . To achieve this goal one has to scan over these two parameters. Specifically, we need QCD predictions for the eight lowest-lying bottomonium resonances in a reasonable range of values for the bottom mass and the strong coupling constant. For the former we will consider the recent determination of Ref. [26], m b = 4.2 GeV, and vary the MS mass between 4 and 8 GeV in steps of 500 MeV. Since the value of α s that enters the perturbative expansion in Eq. We generate QCD predictions for the bottomonium spectrum varying the two renormalization scales µ and R in a correlated way, as explained in Ref. [26]. This makes sure theoretical correlations are properly propagated but avoids the so called d'Agostini bias [66]. The range in which the two scales are varied comes from analyzing the argument of the logarithm in Eq. (7.2) in the MSR scheme, such that its value varies between 1/2 and 2. Such range will depend on the value of the bottom mass, and therefore we will adapt the results of Ref. [26] for the values covered in our scan. We find the upper limit of both µ and R is independent of m b , and fix it to 4 GeV. The lower bound of these two scales is set equal and found to increase approximately linearly with m b . We find the following approximate expressions : where the subindex 1, 2 denotes the principal quantum number of the heavy quarkonium state. Since both values have a positive slope, the range in which scales are varied decreases as m b increases, and larger scales are probed. This renders smaller perturbative uncertainties for larger values of the bottom mass, as expected. In Fig. 9 the dependence of the mass of the two lowest-lying vector resonances with m b is shown graphically. We remove its main contribution, namely twice the bottom mass in the MSR scheme at the scale R = 1 GeV to better see how the uncertainty shrinks as we increase m b . We also observe that these residual masses slightly decrease as the QCD bottom mass increases.  Once we have generated our ensemble of bottomonium masses for a given bottom mass, we will determine the parameters of the Cornell model that best reproduce the QCD prediction. We generate QCD pseudo-data at LO, NLO, N 2 LO and N 3 LO, which can be viewed as a set of highly correlated experimental measurements. This makes a regular χ 2 fit with a non-diagonal covariance matrix impossible due to the d'Agostini bias. But if we do not include the theoretical covariance matrix there are no other uncertainties left and it is not even possible to write down a χ 2 function. Therefore we will simply make a statistical regression analysis, which will provide us with a dispersion uncertainty. The strategy is then very similar to that of Ref. [26] : the QCD renormalization scales are varied in a correlated way in terms of two dimensionless variables µ n = µ min n + x (4 GeV − µ min n ), R n = µ min n + y (4 GeV − µ min n ), with n = {1, 2}, 0 ≤ {x, y} ≤ 1 and µ i defined in Eq. (8.1). We define our regression χ 2 function as follows : where in practice χ 2 min (x, y) is only known after the minimization is carried out, but renders the χ 2 function dimensionless, and yields the right dimensions to the covariance matrix and parameter uncertainties. With this definition, the reduced χ 2 is exactly one at the minimum. For a given value of m b we scan over all possibles values of {x, y}, and for each pair we determine the parameters of the Cornell model. The average of the best-fit values and regression uncertainties become the central value and fit uncertainty, respectively, while the semi-sum of the largest and smaller values attained for each parameter in the scan is taken as the theoretical uncertainty. The latter dominates over the fit uncertainty. After this procedure is repeated for all values of m b and α s , we obtain the Cornell model parameters as functions of these QCD fundamental quantities.

Results of the Calibration
We start by showing details of the calibration for a given value of the MS bottom mass, using the world average value for α s . We take the value m b = 4.2 GeV, as obtained from NRQCD fits to bottomonium states, which corresponds to m MSR b (R = 1 GeV) = 4.701 GeV at four loops. Given these values of the bottom mass and the strong coupling constant, at each order we get 3719 triads of best-fit values. A useful way of showing these results is through histograms, in which, after appropriate binning, one can see how often a given value of a parameter is generated. For the Cornell mass parameter, we choose our bins equally spaced, with bin-size 10 MeV. A closer look into the histograms reveals that there are large, nearly unpopulated tails extending towards low mass values. This would be of no concern if the renormalization scales were parameters that could be varied in a Gaussian way, since the average and standard deviation automatically dump such effects. This is not our case, and therefore these tails make the uncertainties unnaturally large, and can bias the central values. Therefore we cut off bins whose height is less than 8 % of the highest bin. This translates into discarding between 4 % and 15 % of the best-fit values for each order, leaving always ensembles of more than 3000 elements. The resulting histograms can be viewed in Fig. 10, together with the corresponding value of the MSR mass at the reference value of 1 GeV, which happens to be always within the covered values of the Cornell mass. One can see a very clear peak at the highest order, which gets somehow broader towards lower orders. At N 2 LO we observe two maximums, one of them much narrower and higher than the other. Except at lowest order, the distributions are not symmetric. However, the semi-sum of the maximum and minimum values attained in the scan and the average of all points are very close, being the difference much smaller than perturbative uncertainty. For simplicity we will use the average. Finally, we use the average of the individual (rescaled) fit uncertainties as the global fit uncertainty. The results of this procedure at various orders is shown in Table. 3.
We move to the (arguably) most interesting result of this article, the dependence of the Cornell mass parameter with a short-distance QCD mass. Even though (for convenience) we have generated our QCD predictions in terms of the MS mass, this scheme is far from a kinematic mass, and should rather be thought as a coupling constant. On the other hand the MSR mass for small values of R is a kinetic mass free from renormalon ambiguities, and given the results of Secs. 5 and 6, we will calibrate the Cornell mass against m MSR b (R = 1GeV). As expected, we observe a linear relation between these two parame-   ters, and we find that the slope is very close to (and compatible with) unity : 0.995 ± 0.026 with an intersect compatible with zero within uncertainties : 0.05 ± 0.19. 9 This pattern is found for every value of α s and also at various orders, but for simplicity we show the linear relation only at N 3 LO and for the world average in Fig. 11. We believe this procedure is justified since individual uncertainties are almost 100 % correlated. We have checked that the difference of the regular and weighted average of the individual central values is of the order of half an MeV. In Fig. 12(a) we show the value of the difference between the Cornell and MSR masses at various orders. We find nice convergence and decreasing error bars as the perturbative information is increased, while at each order the result is compatible with zero. Similarly, Fig. 12  . Histograms for these two parameters look less peaky than in Fig. 10, and therefore we trim only those values which appear less frequently than 2 % of the highest peak. Figs. 13(a) and 13(b) show this dependence for the α s Cornell parameter as blue dots with error bands, together with the QCD coupling, depicted as a red solid line, evaluated at the non-relativistic scale µ NR . This scale is chosen such that the argument of the logarithm in Eq. (7.2), once the pole mass is expressed in terms of the MSR scheme, becomes one. Since our fit includes n p = 1, 2 states, we take the average valuen p = 3/2. Also we choose R = µ and determine µ NR by solving numerically the following equation : (d) Figure 13. Dependence of the Cornell parameters α Cornell s (upper two panels) and σ (lower two panels) with the MSR bottom mass (leftmost two panels) and the QCD coupling constant at the Z-pole (rightmost two panels). The two upper plots also show, with a solid red line, the strong coupling constant evaluated at a characteristic non-relativistic scale. (µ NR ) within perturbative uncertainties is observed when scanning either over the bottom mass or the strong coupling constant reference value. Not only the order of magnitude is the same, but also its dependence on m b and α s (m Z ) follows the same pattern (although within uncertainties the dependence on m b could be considered flat). If we take n p = 1(2) the value of α s decreases (increases) roughly by 10 %, leaving our conclusions unchanged. Needles to say, expecting a one-to-one correspondence between these parameters is too naïve, but this simple analysis disfavors quark model analyses in which α Cornell s is assigned a QCD-like running, evaluated at the reduced mass of the QQ pair. However, given that the QCD running is logarithmic and depends only on the ratio of initial and final scales, as long as the bottom pair reduced mass is taken as the boundary condition to obtain α Cornell s for a charmonium analysis (or vice-versa), no serious mistake is committed.
We show the dependence of σ with the bottom mass and the strong coupling constant reference value in Figs. 13(c) and 13(d). The dependence of σ with the former shows exactly what one would expect from a static potential parameter : there is no mass dependence at all and one could simply take the average of all points, obtaining σ Cal. = 0.176 ± 0.088 GeV 2 . This value compares well with the one obtained from the fit to the bottomonium experimental data in Eq. (3.1) σ fit = 0.207 ± 0.011 GeV 2 . It appears there is some dependence on the value of α (n f =5) s (m Z ), although drawing strong conclusions is not possible given the size of the uncertainties, which grow for larger values of α s . In any case one expects some dependence of σ with α s , since as argued in Refs. [9,10], the linear rising term in the static potential is of perturbative nature.

Conclusions
In this article we have confronted a simple version of the Cornell model with both experimental data and QCD. This model contains three terms, with one parameter associated to each one of them : a constituent quark mass, a linear raising potential, and a Coulomb-like potential. We have solved numerically the Schödinger equation for the Cornell model using the Numerov algorithm, performing several checks to make sure the accuracy of the approximation is much smaller than any other uncertainty involved. The Cornell model includes the static approximation (solved exactly) plus the leading non-relativistic correction, in the form of 1/m 2 Q suppressed (angular-momentum-dependent) potentials, whereas in QCD we include as many non-relativistic corrections as necessary to achieve up to N 3 LO accuracy.
As a warm-up exercise we have determined the Cornell potential parameters for the bottomonium and charmonium systems from fits to the 8 states with lower mass. We have confirmed that this simple version of the Cornell model fails to predict states with higher mass, possibly because it does not include string-breaking effects, but gives reasonable post-dictions of the masses of the states that enter the fit. NRQCD can predict the first 8 bottomonium states and the first 2 charmonium states within perturbation theory, and therefore it is possible to calibrate the Cornell model using bottomonium QCD predictions. Of course this only makes sense if the Cornell model is solely of perturbative nature.
The QCD static potential suffers form a u = 1/2 renormalon, identical to that of the quark pole mass up to a factor of 2 and a sign. Therefore one can cancel the renormalon in the static energy (the sum of the static potential and twice the pole quark mass) by expressing the quark mass in a short-distance scheme. Since we are dealing with a thresholdlike problem (bound states of a quark-antiquark pair), a low-scale short-distance mass should be used. To keep small logs in the static potential and in the subtraction series relating the pole and short-distance masses, it is compulsory to use a scheme with a tunable subtraction scale. The MSR mass is a scheme that satisfies these two criteria, and has already been used in the context of quarkonium. Therefore we also employ this scheme to make predictions for the static energy. Furthermore, we use R-evolution to sum-up large renormalon-type logarithms in the static potential, whose argument is the ratio of the subtraction and renormalization scales. In this way, we define an R-improved MSR static potential, which shows nice order-by-order convergence, and has the same linear rising behavior as the Cornell model. In fact, if the renormalon subtraction scale is chosen close to 1 GeV, the Cornell and R-improved static potential at N 3 LO nicely agree for moderate and large values of r. Since the linear raising potential is of perturbative nature and agrees with that of the Cornell model, we conclude that all ingredients in this model for bottomonium are perturbative and a calibration against QCD makes sense.
We have compared our R-improved MSR static potential with lattice simulations of the same quantity, performing a fit for the strong coupling constant and the renormalon subtraction reference scale. After the fit, we find very nice agreement with lattice QCD results for the entire set of r values covered in the simulation, which includes distances as large as 0.84 fm. To make this agreement possible, it is essential to use a "profile function" for the renormalization scale µ, which freezes to 1 GeV for values of the radius larger than 1 fm.
To calibrate the Cornell model we generate templates for the 8 lightest bottomonium bound states with NRQCD. These predictions depend on two scales : the renormalization scale µ and the renormalon subtraction scale R. These take different values for the various bound states, but are varied in a correlated way. For a given value of the bottom quark mass and strong coupling constant we generate templates at LO, NLO, N 2 LO and N 3 LO, which contain the QCD prediction for the masses of the bound states in a two-dimensional grid of µ and R. To avoid the d'Agostini bias, we adjust the Cornell model parameters to every entry on each template, effectively scanning over the two renormalization scales. Since there are no uncertainties in the template that can be used to construct a χ 2 function, we use a regression algorithm to assign "fit-incertitudes" to the adjusted parameters, normalizing the χ 2 such that equals the number of degrees of freedom at its minimum. On top of these, there are theoretical uncertainties from the scale scan. To figure these out, we first trim away values of the Cornell model parameters that are found in the regression less often, discarding those that are less than 8 % (2 %) less frequent than the most likely occurrence for m Cornell b (σ, α Cornell s ). After applying this procedure, we take the average of the remaining values as the central value for the parameter, and half the sum of the maximum and minimum values as the perturbative uncertainty. We also take the average of the regression errors as our final fit uncertainty.
We find an almost perfect linear relation between the Cornell and MSR masses with slope compatible with 1 within 0.19 standard deviations, and if one chooses the MSR reference scale R = 1 GeV the intersect of this relation is compatible with zero within 0.26 standard deviations. This pattern is replicated in a wide range of α s values and for all perturbative orders considered, and nicely complies with the agreement found between the Cornell and static potentials for this particular choice of the reference scale. Our calibration exercise also reveals that the confining parameter σ does not depend on the value of the QCD quark bottom mass, as expected, but the precision of the analysis cannot discard some dependence of this parameter with α s (m Z ). On the other hand, the Coulomb-like Cornell parameter α Cornell s is found to agree within uncertainties with α (n f =4) s (µ NR ), that is with the QCD strong coupling constant evaluated at a typical non-relativistic scale.
Our analysis could be extended in several directions : On the QCD side one could consider more advanced NRQCD predictions, for instance using pNRQCD resummation as done in Ref. [52], or including non-perturbative effects as in Ref. [67]. On the Cornell side, one could consider more sophisticated models, for example incorporating string breaking effects or coupled channels. The calibration itself could be carried out for mixed bc states in which both masses are varied such that their ratio remains constant. As for the comparison with lattice QCD results on the static potential, the next step is studying perturbative uncertainties, order and dataset dependence, incorporating results for other lattice spacings, and taking into account the lattice correlation matrices once they are known.
We close this article raising a concern. We have shown in Sec. 2 that the confining part of the Cornell potential cannot be treated by any means as a perturbation of the Coulomb potential. We have also seen in Sec. 5 that the σ term can be entirely described in perturbative QCD once the renormalon has been canceled. Finally we have presented an analytic formula for the QQ mass in Sec. 7, which is based on solving the Schrödinger equation perturbatively around the lowest order result, that is, the Coulomb potential, including both radiative and relativistic corrections. Therefore one could call into question this perturbative treatment of the static QCD potential since it is known to fail for the not-so-different Cornell potential. On the other hand, pNRQCD power counts perturbative and non-relativistic corrections on an equal footing, as can be seen in Eq. (7.1), where there is no distinction between the former and the latter. One way of shedding light on this apparent puzzle is though a numerical, exact, resolution of the Schrödinger equation for the QCD static potential. A step in this direction has been taken in Refs. [50][51][52].

A Numerical resolution of the Cornell Potential
In this work, the Numerov algorithm [21,22] (or Cowell's method) is used to solve numerically the Cornell potential. Such algorithm can be used to solve ordinary differential equations of second order in which the first derivative does not appear, and therefore is particularly well suited to solve the Schrödinger equation [68]. We follow the specific implementation of Ref. [69], to which the reader is referred to for a complete and pedagogical explanation. We outline the main points of the algorithm in what follows. The method consist in numerically solving the ordinary differential equation in Eq. (2.2) in the range r ∈ [r min , r max ], discretized in N + 1 nodes with step-size h = (r max − r min )/N . This implies r n = r min + n h with 0 ≤ n ≤ N and r n+1 = r n + h. Note that one cannot include the point r = 0 as starting point for a numerical solution since the potential is singular at the origin. Nevertheless we know that for very small r the solution behaves as in Eq. (2.5). The reduced wave function u(r) and the kernel k(r) can be discretized as well, and we use the notation u n = u(r n ), k n = k(r n ). One can Taylor expand u Taking the sum and the difference of the plus and minus equations we isolate even and odd terms, respectively : One can solve for u n in the above equations to find : The point r c will be chosen as the distance where k(r c ) = 0, that is which marks the distance where the classically forbidden region begins and the asymptotic exponential behavior is expected to start dominating. It depends on and n. The condition in Eq. (A.9) as well as r c depend on the energy. For a given value of we use the classic bisection method to find the root of Eq. (A.9), requiring an accuracy of 1 eV, but limiting the number of steps to 1 000. Once the ground state is found, the same method can be used to find the energy of excited states. Iterating the process for various values of the orbital angular momentum one figures out the complete (discrete) spectrum of the static Cornell potential.
Regarding the calculation of the leading relativistic corrections of the Cornell model, the radial integration for the V OGE SS operator can be performed analytically, and we find : Since u n (r → 0) = A n, r +1 this matrix element is non-zero only if = 0, and the limit is simply |A n,0 | 2 . In practice we compute A n,0 from the value of the reduced wave function at the fist node : A n,0 = u n0 (r min )/r min . For rest of operators in Eq. (2.11) we perform a numerical integration between 0 and r max . For the integration we use a Fortran 90 implementation of the QUADPACK package [70]. We reconstruct the function u n (r) with an interpolation using the values computed in the nodes by the Numerov method between r min and r max , while for r < r min we assume it follows the same patter as the first two nodes, u n (r → 0) = u n (r min ) r +1 /r +1 min . For the interpolation we use a Fortran 2008 implementation of the algorithm described in [71].