Simple high-accuracy method for solving bound-state equations with the Cornell potential in momentum space

The well-known Cornell quark-antiquark potential in momentum space contains singularities both in its one-gluon-exchange (OGE) and linear confining parts, which prevents a direct use of the convenient Nystr\"om method to solve the corresponding bound-state integral equation for the meson masses. While it has been known for a long time how the Coulomb-type singularity in the OGE potential can be treated with a subtraction technique, only very complicated methods have been developed to deal with the stronger singularity in the linear potential. In this work, we present a simple subtraction method to remove this singularity from the kernel, such that the Nystr\"om method becomes applicable. Derivatives of the wave function, that appear as a result of the subtraction, are represented by means of interpolating functions, for which we found Lagrange polynomials to be very efficient. Test calculations show excellent agreement with exactly known energy eigenvalues. By increasing the number of integration points and the order of the Lagrange interpolation polynomials, extremely high accuracy can be achieved. This method can also be extended to relativistic Bethe-Salpeter type equations with singular kernels.