You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
abstract: 'A semi-spectral Chebyshev method for solving numerically singular integral equations is presented and applied in the quarkonium bound-state problem in momentum space. The integrals containing both, logarithmic and Cauchy singular kernels, can be evaluated without subtractions by dedicated automatic quadratures. By introducing a Chebyshev mesh and using the Nystrom algorithm the singular integral equation is converted into an algebraic eigenvalue problem that can be solved by standard methods. The proposed scheme is very simple to use, is easy in programming and highly accurate.'
address: 'Institute for Nuclear Studies, Warsaw, Poland'
author:
- 'A. Deloff'
title: 'Quarkonium bound-state problem in momentum space revisited'
---

semi-spectral method,Schrödinger equation, quarkonium 12.39.-x,03.65.Ge ,02.30.Rz

Introduction
============

In a recent work [@deloff] we have advocated the Chebyshev semi-spectral method demonstrating its efficiency in solving some typical differential and integral equations emerging in quantum mechanics. The present paper is in the same vein but here we wish to focus our attention solely on the heavy quarkonium momentum space bound-state problem. Admittedly, the problem is not new but our incentive here is to examine the effectiveness of the semi-spectral approach in solving strongly singular integral equations. Since the latter topic was beyond the scope of [@deloff], this work may be regarded as an immediate continuation of the previous paper.

We would like to believe that the presented method will be useful also outside quantum mechanics, especially that strongly singular integral equations are encountered in many areas of science and engineering. The well known physical applications comprise the quantum mechanical scattering problem, the Omnes formulation [@omnes] of the final-state-interaction, radiative transfer, neutron transport [@chandra] etc. The list of engineering applications is by no means restricted to the widely known aerofoil problem [@wing] and, indeed, many important problems of engineering mechanics like elasticity, plasticity, fracture mechanics, etc. may be also efficiently expressed in terms of singular and hypersingular integral equations. Because it is not always possible to find explicit solutions to the problems posed, much attention has been devoted to approximate methods. It is interesting to note that even when an analytic solution is known, quite often the latter takes the form of a singular integral whose numerical evaluation might be more complicated than a numerical solution of the integral equation.

A hypersingular integral equation arises in quantum mechanics already at a quite elementary level when the linear potential bound-state problem, easily tackled in configuration space, is approached in momentum space. This problem is far from academic since the linear potential plays an important role not only in atomic physics where it is associated with the hydrogen radial Stark effect but also in particle physics serving as a simple confinement model of QCD. Although, in principle, QCD alone should describe the spectroscopy of heavy quarkonia but the implementation of such program is very difficult and instead various phenomenological models incorporating some QCD properties have been developed (for a recent review of quarkonium physics and references to the literature cf. [@cern]). The QCD motivated quark potential models have played a prominent role in understanding quarkonium spectroscopy and are capable of reproducing with surprising accuracy a sizable part of the meson and baryon properties. The non-relativistic potential approach may be justified by the fact that the bottom quark and, perhaps to a lesser extent, also the charmed quark have masses that are large in comparison with $\Lambda\;$– the typical QCD hadronic mass scale. The quark–antiquark potential has been tailored to mock up the properties expected from QCD and the different potential shapes set up in the early days after years of research have evolved to a common form that one might expect from the asymptotic limits of QCD. The prototype for these potentials is still the popular Cornell potential [@cornell] including the one-gluon-exchange Coulomb potential supplemented by a linear potential simulating confinement, as expected from QCD. Therefore, this potential will be also considered in this paper.

Obviously, the non-relativistic potential model can not be pushed beyond certain limits and for systems containing one light quark a complete disregard of relativistic effects might be a serious omission. In addition to that, it was somewhat embarrassing when people realized [@non-rel] that within the non-relativistic formalism the mesons containing a light quark might be more massive than a meson composed with heavier quarks. These difficulties could be aleviated at the expense of a semirelativistic treatment where the relativistic expression for the energy is used. A popular relativistic extension of the Schrödinger equation is the spinless Salpeter equation $$\label{I1}
\left[\sqrt{\bm{p}^2+m_1^2}+\sqrt{\bm{p}^2+m_2^2} + V(r)\right]\,\Psi(\bm{r}) = E\,\Psi(\bm{r})$$ where $m_1,\, m_2$ are the quark masses, $\bm{p}$ is the c.m. momentum, $V(r)$ denotes the quark-antiquark potential and $E$ is the eigenenergy. Since in such case the Laplacian operator appears under a square root, the coordinate space is rather unwieldy for solving the bound state problem and the momentum space seems to be the most natural alternative. Indeed, in momentum space the energy operator is diagonal and the difference in computational effort between non-relativistic and semi-relativistic treatment is minor. Although the momentum space approach solves some problems automatically but at the same time it does create another difficulty in that the quark-antiquark potential gives rise to a singular kernel in the appropriate integral equation. Whilst the Coulomb potential yields a kernel with a logarithmic singularity that can be removed by subtraction [@lande], the kernel associated with a linear potential exhibits a double-pole singularity for which the subtraction scheme is insufficient. To clarify this important point let us consider just the linear potential for simplicity restricting our attention to a zero orbital momentum state. The potential term that enters the appropriate wave equation involves the integral with a double pole singularity $$\int_0^\infty\dfrac{k^2\,\phi(k)\,{\textstyle{\rm d}}k}{(k^2-p^2)^2}=\displaystyle
\int_0^\infty\left\{k^2\,\dfrac{\phi(k)-\phi(p)}{k^2-p^2}-{\textstyle\frac{1}{2}}p\, \phi^{\,\prime}(p)\right\}
\dfrac{{\textstyle{\rm d}}k}{k^2-p^2},
\label{I0}$$ where $\phi(k)$ is the wave function, $\phi^{\,\prime}(p)$ denotes the derivative and $p$ is a real parameter. It may be easily verified that the two extra terms occurring on the right hand side of (\[I0\]) can be supplemented with impunity because the integrals multiplying, respectively, $\phi(p)$ and $\phi^{\,\prime}(p)$ are both bound to vanish. The integral on the right hand side is non-singular and in the limit $k\to p$ the integrand goes to a finite limit ${\textstyle\frac{1}{2}}\phi^{\,\prime}(p)/p+\frac{1}{8}\phi^{\,\prime\prime}(p)$. This demonstrates explicitly that by using a subtraction technique it is perfectly possible to remove the singularity converting the integral to a form amenable for computation. Nevertheless, the subtraction scheme (\[I0\]) would be insufficient for solving an integral equation as it introduces unknown first $\phi^{\,\prime}(p)$ and second derivative $\phi^{\,\prime\prime}(p)$ at the top of the unknown function. This also explains why the Nystrom method, which has been rather efficient in solving the Coulomb bound state problem in momentum space [@lande], does not work for the linear potential. ( The calculation using Nystrom method presented in [@nys] is incorrect because the infinite diagonal term in the potential matrix has been simply omitted whereas the proposed correction, given in their eq. (34), is proportional to a logarithmically diverging integral.)

In the early attempts to overcome this difficulty the singularity was removed by hand, by introducing an arbitrary cut-off [@chiu][@eyre] in the potential. The resulting non-singular integral equation involving the modified potential could be then solved by standard methods. The unwelcome arifacts of the cutoff might be eventually disposed of by perturbative methods [@eyre]. However, a more promising approach is to seek the wave function in the form of an expansion in terms of a complete set of orthogonal basis functions. The most common choice here has been the oscillator or Sturmian basis both of which have analytic Fourier-Bessel transforms making them well suited in calculations where it is advantageous to work in configuration and momentum space simultaneously. In a variational Ritz-type approach the upper bounds of the true eigenvalues could be computed by diagonalizing the corresponding Hamiltonian matrix (cf. [@var1], [@var2],[@var3]). The expectation values of the energy can be evaluated in momentum space and the potential expectation values in configuration space. The expansion method could be used in a similar fashion to solve the momentum space integral equation by means of the Galerkin method [@spline], [@maung]. With a judicious choice of the basis functions, the singular integrals can be calculated analytically, or numerically. Note, that in this case the integrand is a known function and, therefore, the subtraction technique, like the one outlined in (\[I0\]), is fully applicable. There are also non-variational approaches based on eigenfuction expansion such the collocation method [@spline],[@jean], or the Multhopp [@wing][@BB] technique. Keeping $N$ terms of the truncated expansion, the $N$ expansion coefficients can be determined from the requirement that the integral equation be exactly satisfied at $N$ distinct values of the momentum variable. The semi-spectral Chebyshev method developed in this paper also belongs to the last group. However, the Chebyshev series, after reshuffling takes the form of an interpolative formula. In consequence, the expansion coefficients and the function values taken at the mesh-points are connected by a linear relation (cf. [@deloff]). Thus, put in a nut-shell, the underlying idea is to solve the integral equation [*exactly*]{} on the Chebyshev mesh and, subsequently, interpolate by means of a high degree polynomial. The plan of the presentation is as follows. In the next section we set the necessary background deriving the hypersigular integral equation associated with the Coulomb-plus-linear potential in momentum space. Upon introducing the Chebyshev mesh and using the interpolative formula for the wave function, the integral equation is converted into an algebraic eigenvalue problem. This is the ultimate form because the eigenvalue problem can be solved with the aid of standard library procedures. Section 3 is devoted to a numerical test where we compare the momentum space calculations with the results obtained by solving the Schrödinger equation in configuration space. Finally in the last section we present our conclusions.

Solution of the singular integral equation
==========================================

The Coulomb-plus-linear potential considered in this paper is $V(r)= V^{(C)}(r)+V^{(L)}(r)$ with $$V^{(C)}(r)= -\alpha/r; \quad V^{(L)}(r)=r/a^2
 \label{prel1}$$ where the ”coupling” $\alpha$ is dimensionless and the parameter $a$ has a dimension of length ($\hbar=c=1$ units are adopted hereafter). Both parameters are assumed to be provided. In momentum space the wave function $\phi_\ell(k)$ with orbital momentum $\ell$ obeys the partial wave Schrödinger equation $$(E- k^2/2\mu)\;\phi_\ell(k)=
\int_0^\infty\,V_\ell(k,k')\,\phi_\ell(k')\,k^{\prime\, 2}\,{\textstyle{\rm d}}k'
  \label{prel2}$$ where $\mu$ is the quark-antiquark reduced mass, $E$ is the binding energy and $V_\ell(k,k')$ denotes the $\ell$-th partial wave projection of the local potential $V(r)$ $$V_\ell(k',k)= \dfrac{2}{\pi}
        \int_0^\infty j_\ell(k'r)\,V(r)\,j_\ell(kr)\,r^2\,{\textstyle{\rm d}}r,
          \label{prel3}$$ where $j_\ell(x)$ is the spherical Bessel function [@abramowitz]. Strictly speaking, upon inserting (\[prel1\]) in (\[prel3\]), we obtain a divergent integral but a customary regularizing procedure to overcome this difficulty is first to multiply $V(r)$ by a screening factor ${\textstyle{\rm e}}^{-\eta r}$ enforcing convergence and then set $\eta\to 0$ in the result. Applying this procedure, the Fourier transform (\[prel3\]) of a power–law potential $v(r)=r^{2n-1},\;n=0,1,2,...$ can be effected in an analytic form [@BB] $$\lim_{\eta\to 0}\,\dfrac{2}{\pi}
     \int_0^\infty j_\ell(k'r)\,{\textstyle{\rm e}}^{-\eta r}
        \,r^{2n+1}\,j_\ell(kr)\,{\textstyle{\rm d}}r=
  \dfrac{(2n)!}{2^n\,n!\,\pi (k\,k')^{n+1}}\; Q_\ell^n (z)
          \label{prel4}$$ where $z=(k^2+k^{\prime\,2})/2kk'$ and the $Q_\ell^n(z)$ denotes $n$–th derivative of the Legendre function of the second kind with respect to the argument $z$ (formula (5) in [@BB] contains a misprint). Setting $n=0$ and $n=1$ in (\[prel4\]) we obtain, respectively, the kernels for the Coulomb (C) and the linear potential (L) $$V^{(C)}_\ell(k,k')=-\alpha \, Q_\ell(z) /(\pi kk');
 \quad
 V^{(L)}_\ell(k,k')= Q_\ell^\prime(z)/[\pi (akk')^2].
\label{prel5}$$ The Coulomb part of the kernel exhibits a logarithmic singularity for $k'=k$ contained in the Legendre function. Indeed, the latter can be written as $$Q_\ell(z)=P_\ell(z)\,Q_0(z) - w_{\ell-1}(z)
\label{prel6}$$ where $$Q_0(z)={\textstyle\frac{1}{2}}\log|(1+z)/(1-z)|=\log|(k+k')/(k-k')|
\label{prel6a}$$ with $ P_\ell(z)$ being a Legendre polynomial. It is understood that the last term in (\[prel6\]) should be absent for $\ell=0$ whereas for $\ell>0$ it assumes the form of a polynomial in $z$ (cf. [@abramowitz]) given by the expression $$w_{\ell-1}(z)=
 \displaystyle\sum_{n=1}^\ell\frac{1}{n}P_{n-1}(z)\,P_{\ell-n}(z).
\label{prel7}$$ The kernel associated with the linear potential given in (\[prel5\]), in addition to the logarithmic singularity, exhibits also a second order pole, as may be seen by performing explicitly the differentiation in (\[prel6\]) $$Q_\ell^\prime(z)=P_\ell^\prime(z)\,Q_0(z)
 + P_\ell(z)\,Q_0^\prime(z)
 - w_{\ell-1}^\prime(z)
\label{prel8}$$ with $$Q_0^\prime(z)=\dfrac{1}{1-z^2}
=-\left(\dfrac{2kk'}{k'+k}\right)^2\dfrac{1}{(k'-k)^2}.
\label{prel8a}$$ The second order pole given by (\[prel8a\]) can be eliminated from the integral equation (\[prel2\]) and to this end integration by parts is applied to this term. Quite generally, this procedure gives $$\int_0^\infty \dfrac{f(k,k')\,\phi_\ell(k')\,{\textstyle{\rm d}}k'}{(k'-k)^2}=
\int_0^\infty \dfrac{{\textstyle{\rm d}}k'}{k'-k}\,  \dfrac{\partial }{\partial k'}\, \left[f(k,k')\,\phi_\ell(k')\right]
\label{prel9}$$ where the unspecified function $f(k,k')$ needs to be integrable. The above formula holds because the wave function $\phi_\ell(k')$ vanishes when $k'$ tends to either of the integration end points. The resulting Cauchy principal value integral in (\[prel9\]) can be computed by using the dedicated Chebyshev quadrature given in [@deloff]. Nevertheless, the lowering of the order of the pole outlined above has its price and in the integral on the right hand side of (\[prel9\]) the derivative of the unknown wave function will appear. As we shall see in a moment, the semi-spectral Chebyshev method is well suited to handle such situation.

It will be convenient for us using $1/a$ as the unit of energy, passing to dimensionless quantities: $ \epsilon\equiv Ea,\;x\equiv ka,\;x' \equiv k'a$. The resulting integral equation $$\label{prel10}
\begin{split}
\left(\epsilon -\dfrac{x^2}{2\mu a} \right)\phi_\ell(x)=
\dfrac{1}{\pi x^2}\int_0^\infty\left\{P_\ell^{\,\prime} (z) \, \log\left|\dfrac{x'+x}{x'-x}\right|
-w_{\ell-1}^{\,\prime}(z)\right\}\phi_\ell(x')\,{\textstyle{\rm d}}x'
\hspace*{.2\textwidth}
\\ \mbox{}-
\dfrac{4}{\pi}\int_0^\infty
\dfrac{{\textstyle{\rm d}}x'}{x'-x}\;\left\{\chi_\ell(x')+\phi_\ell(x')\,\dfrac{\partial}{\partial x'} \right\}
 \dfrac{x^{\prime\, 2} P_\ell(z)}{(x'+x)^2}-
\hspace*{.2\textwidth}
\\ \mbox{}-
 \dfrac{\alpha}{\pi x}\int_0^\infty\left\{P_\ell(z)\, \log\left|\dfrac{x'+x}{x'-x}\right|
-w_{\ell-1}(z)\right\}\phi_\ell(x')x'\,{\textstyle{\rm d}}x' 
\hspace*{.2\textwidth}
\end{split}$$ involves two dimensionless parameters: $\alpha$ and $2\mu a$. Prime on a function of $z$ denotes in (\[prel10\]) the derivative with respect to the argument. The derivative of the wave function appearing in the integrand of the second integral in (\[prel10\]) has been regarded as an additional function $\chi_\ell(x)$ to be determined. In order to complete our scheme the integral equation (\[prel10\]) needs to be supplemented with a complementary equation $${\textstyle{\rm d}}\phi_\ell(x)/{\textstyle{\rm d}}x=\chi_\ell(x)
\label{prel11}$$ and we end up with two equations for two unknown functions: $\phi_\ell(x)$ and $\chi_\ell(x)$.

The system (\[prel10\])–(\[prel11\]) is amenable for computation and the integral equation will be turned into a finite matrix equation. As a preliminary step, the semi-infinite domain of the independent variable $x$ will be mapped onto a finite interval $(-1,1)$. Among endless possibilities perhaps the simplest is the rational mapping $$\label{prel12}
x=\sigma (1+t)/(1-t),$$ where $t\in (-1,1)$ and $\sigma$ is a numeric parameter at our disposal providing additional control of the rate of convergence. We tried some other mappings, specifically trigonometric ($x=\sigma \tan[(\pi/4)(1+t)]$), or logarithmic ($x=\sigma \log[(3+t)/(1-t)]$) but they did not bring noticible improvement in the problem under consideration. The semi-spectral Chebyshev method uses Chebyshev polynomials as the basis functions. The Chebyshev polynomial of the first kind $T_N(t)$ of the order $N$ is defined by the formula $$T_N(t)=\cos[N\, \arccos(t)]
\label{prel15}$$ and has $N$ zeros in the interval $(-1,1)$, located at the points $$\label{prel13}
t_i = \cos[\pi(i-{\textstyle\frac{1}{2}})/N];   \quad    i=1,2,...,N.$$ In the following the variable $t$ will be discretized by using the classical Chebyshev mesh (\[prel13\]) in which case $N$ becomes the order of approximation to be selected by the user. The semi-spectral Chebyshev method interpolates the unknown function $f(t)$ on the Chebyshev mesh (\[prel13\]) $$f(t)=\sum_{i=1}^N f(t_i)\,G_i (t),  
\label{prel14}$$ where $G_i(t)$ denotes the cardinal function with the property $G_i(t_j)=\delta_{ij}$. These functions can be constructed as superpositions of Chebyshev polynomials $$G_j(t)=\dfrac{2}{N}\sum_{i=1}^N {\rm '\mbox{} } \ T_{i-1}(t_j)\,T_{i-1}(t),
\label{prel16}$$ where the primed sigma denotes a summation in which the first term should be halved. By taking advantage of the interpolative formula (\[prel14\]), the differentiation or integration of a function reduces to differentiation or integration of Chebyshev polynomials which in most cases is elementary and can be performed in an analytic form. In consequence, the array containing the values of the derivative computed at the grid-points will be connected to similar array representing the function by a linear transformation $$\left\{\dfrac{{\textstyle{\rm d}}f(t)}{{\textstyle{\rm d}}t}\right\}_{t=t_i}=\sum_{j=1}^N D_{ij} \, f(t_j),\quad i=1,2,...,N
\label{prel17}$$ where $D_{ij}$ is easily computed numerical matrix (cf. [@deloff]). There are also various integration rules available. Assuming that the function $f(t)$ is non-singular in the integration domain, we have $$\int_{-1}^1 f(t) \,{\textstyle{\rm d}}t =\displaystyle \sum_{i=1}^N w_i\, f(t_i),
\label{prel18}$$ which is Gauss-Chebyshev integration in which the weighting function is equal to unity. The weights $w_i$ are all positive and their sum equals to 2. Similar rules can be derived for singular integrals. The Cauchy principal value integration can be performed using the automated quadrature rule $$\int_{-1}^1\dfrac{ f(t) \,{\textstyle{\rm d}}t}{t-\tau} = \displaystyle\sum_{i=1}^N \omega_i(\tau)\, f(t_i),
\label{prel19}$$ where it is assumed that $\tau\in<-1,1>$. When $\tau$ coincides with either of the integration end-points the integral is undefined. The dedicated weighting functions $\omega_i(\tau)$ can be calculated analytically and exhibit logarithmic end-point singularity for $\tau=\pm 1$. Similar rule can be obtained for a weakly singular integral $$\int_{-1}^1 f(t) \,\log|t-\tau|\,{\textstyle{\rm d}}t = \displaystyle\sum_{i=1}^N \Omega_i(\tau)\, f(t_i),
\label{prel20}$$ where it is assumed that $\tau\in(-1,1)$. In contrast with the previous case, $\log|t-\tau|$ singularity is integrable and the dedicated weighting functions $\Omega_i(\tau)$ do exist even when $\tau$ coincides with either of the integration end-points. For explicit analytic expressions for all of the weighting functions introduced above the reader is referred to [@deloff].

To arrive at the ultimate finite matrix eigenvalue problem, as the first step, we map both, the external ($x$), and the internal ($x'$) variable onto the $(-1,1)$ interval with the aid of (\[prel12\]). Subsequently, the problem is discretized by putting the external variable on the Chebyshev mesh (\[prel13\]), at the same time replacing all integrations in (\[prel10\]) by summations, following the appropriate Chebyshev rules listed above. In practice this procedure leads to a chain of substitutions to be made in the integrals occurring in (\[prel10\]), [*viz.*]{} $$x\to x_i =\sigma(1+t_i)/(1-t_i);\quad     \phi_\ell(x)\to\phi_\ell(x_i)\equiv X_i;$$ and $$x'\to x_j =\sigma(1+t_j)/(1-t_j);\quad     \phi_\ell(x')\to\phi_\ell(x_j)\equiv X_j;$$ where $X_i$ are the unknown mesh values of the wave function to be determined. The derivative $\chi_\ell(x_j)$ is eliminated in favor of $X_j$ with the aid of the $D_{ij}$ matrix, accounting for the change of variables $$\chi_\ell(x_j) = \dfrac{(1-t_j)^2}{2\sigma}\sum_{k=1}^N D_{jk}\,X_k.$$ Further substitutions associated with integration, respectively, are $${\textstyle{\rm d}}x'                           \to  2\sigma w_j/(1-t_j)^2,$$ for non-singular integrals $$\dfrac{{\textstyle{\rm d}}x' }{x' - x}                \to  \omega_j   (t_i)\,   \dfrac{1 - t_i }{1 - t_j },$$ for principal value integral, and $$\log \left| \dfrac{x'+x}{x'-x} \right|\, {\textstyle{\rm d}}x'     \to  2 \sigma \,\dfrac{w_j\,\log|1-t_i\,t_j| - \Omega_j (t_i)}{(1-t_j)^2}$$ for integrals involving logarithmic singularity. Finally, all integrations will be effected by carrying out a summation over $j$. It is worth noting that the diagonal terms $i=j$ are always finite and all singularities are under control.

When the indicated above manipulations have been accomplished, we end up with a homogeneous system of $N$ algebraic equations in which the $N$ unknowns are the mesh-point values of the wave function ($X_j$) and the Schrödinger equation takes the desired finite matrix form $$\sum_{j=1}^N \left(
  V_{ij}+\dfrac{x_i^2}{2\mu a}\, \delta_{ij} -\epsilon \right)\,X_j = 0.
\label{prel21}$$ The non-symmetric matrix $V_{ij}$ represents here the potential and results from evaluating the integrals occurring on the right hand side of (\[prel10\]) (the explicit form of $V_{ij}$ is rather lengthy and will not be quoted here). When the kinetic energy term is lumped together with $V_{ij}$ into a single matrix, eq. (\[prel21\]) presents a standard algebraic eigenvalue problem. If need arises, the non-relativistic Schrödinger equation (\[prel21\]) can be easily converted to the relativistic form (\[I1\]) in the center-of-mass frame by changing just the kinetic energy term $$x_i^2/(2\mu a)\to \sqrt{x_i^2+(a\,m_1)^2}+\sqrt{x_i^2+(a\,m_2)^2}-a\,(m_1+m_2).$$ Our calculational scheme is now complete and for assigned values of $\ell$ and two dimensionless parameters $s\equiv 1/2\mu a$ and $\alpha$ specifying the strength of the two potentials in (\[prel1\]), we are in the position to determine numerically the value of the binding energy $\epsilon(\ell,s,\alpha)$. In the particular case $\ell=0$ and $\alpha=0$ the exact result is known and the binding energy is $\epsilon(0,s,0)=-s^{2/3}\, z_\nu$ where $z_\nu$ with $\nu=1,2,3,...$ denotes a zero of the Airy function $Ai(z)$ (cf. [@abramowitz]).

Numerical test
==============

We start the numerical test with the Coulomb bound state problem leaving out the first two integrals on the right hand side of (\[prel10\]). The hydrogen-like bound state problem in momentum space has already been considered in [@deloff] but to determine the bound states we solved the secular equation. It is therefore of interest to repeat the Coulomb bound-state calculation in which the energy spectrum is obtained by solving the algebraic eigenvalue problem (\[prel21\]). The latter procedure is much simpler as there is no need to solve a transcendental equation. In all our computations we were using the linear algebra package LAPACK [@lapack] as our eigenvalue solver. The results for the Coulomb potential are displayed in Table \[table1\]. Since in this case the exact eigenenergies are known analytically we present the absolute value of the relative error on each level as a function of the mesh size $N$. The nodal quantum number $n$ enumerates the the different bound states for a fixed $\ell$ with $n=0$ corresponding to the ground state. We wish to recall that with non-symmetric matrices the accuracy of the standard library procedures is believed to be not as good as in the case of symmetric matrices. Nevertheless, as seen from Table \[table1\], the convergence rate is exponential and $N=80$ is sufficient for securing machine accuracy. There are not very many methods available that would be capable of achieving such a high precision. For comparison, in the last raw (entries in parenthesis) we give the relative error corresponding to the traditional method using the subtraction scheme [@lande] in which case the resulting eigenvalue problem is symmetric. The advantage of the semi-spectral method is manifest.

As our next test we take on the linear potential setting $\alpha=0$ in (\[prel21\]) and putting for simplicity $s=1$ in our computations. The resulting binding energies $\epsilon$ for different $\ell$ values are displayed in Table \[table2\] using the same conventions as in Table \[table1\]. For $\ell=0$, as the exact values we take the zeros of the Airy function tabulated in [@abramowitz]. For $\ell>0$ the values marked as exact have been computed by solving the appropriate Schrödinger equation in configuration space. For this purpose we used the ingenious algorithm developed in [@wien]. The code from [@wien] has been revamped for obsolescent features and the original Runge-Kutta driver advancing the solution from $x$ to $x+h$ has been replaced by a more accurate driver based on Chebyshev approximation as described in [@deloff]. After the above changes, the typical relative error in all considered here cases was estimated to be of the order of $10^{-11}$. As a cross-check, we succeeded in reproducing the exact results for $\ell=0$ up to eleven significant digits. To obtain the entries in table \[table2\] for each $\ell$ value we needed to solve the algebraic eigenvalue problem (\[prel21\]) and in nearly all considered here cases we managed to get seven significant figures which is more than adequate in all practical applications. Our results have been obtained keeping quite moderate approximation order $N=100$. Only the $\ell=0$ case which was more stubborn forced us to go to larger $N$. It is apparent from table \[table2\] that the solutions are very stable with respect to increasing $N$ albeit the rate of convergence is no longer exponential. In fact, it is quite slow when compared with the Coulomb case. Making such comparison, however, it has to be kept in mind that in the linear potential case we need to determine two unknown functions (wave function and its derivative) rather than one and therefore $N$ should have been doubled if we wanted to keep the same number of points per function. Other than that, there is probably a good deal of cancellation across the pole and this might be responsible for some loss of accuracy.

Finally, we are going to consider the case where both, the Coulomb and the linear potential are present. The quark-antiquark potential has been adopted from a realistic study [@quigg] of charmonium $(c\bar{c})$ and bottomium $(b\bar{b})$ $V(r)=-\alpha/r+\beta\,r$ where we stick to the parameter values provided in [@quigg], namely $$\alpha=0.50667,\quad \beta=0.1694\;GeV^2,\quad m_c=1.37\;GeV,\quad m_b=4.79\;GeV. 
\label{q}$$ The results of our computations are presented in Table \[table3\]. The quarkonium masses $M$ displayed there have been obtained from the expression $M=2m_q+E$ where $m_q$ is the quark(antiquark) mass. To determine the binding energy $E$ the appropriate non-relativistic Schrödinger equation was solved in both, the momentum and the configuration space. As seen from Table \[table3\] there is excellent agreement between these two approaches.

Summary and Conclusion
======================

The aim of this paper was to demonstrate the strength of the semi-spectral Chebyshev method in solving integral equations whose kernels exhibit singularities of the Cauchy or the logarithmic type. Such equations may be encountered in quantum mechanics as has been exemplified by considering the Coulomb-plus-linear potential bound state problem in momentum space. The latter problem is considered in this work for illustrative purposes and therefore we have gone in some details. The semi-spectral Chebyshev method has many advantageous features. First, it is very easy to use since it is based on a polynomial interpolation where both, the mesh and the polynomials, can be readily obtained in an analytic form. Second, the programming is exceedingly simple because differentiation or integration of polynomials can be performed analytically and on a mesh these operations take the form of matrix multiplications. The presented method is well suited to handle singular integral equations (with Cauchy or logarithmic singularities) because automatic quadratures are provided for evaluating singular integrals. This allows for a quick and seamless discretization and since the integrals involving singular kernels have finite diagonal elements the Nystrom method is still applicable. Ultimately, the integral equation is converted into an algebraic eigenvalue problem which can be solved directly by standard library procedures. There is no need to solve a complicated transcendental equation. Third, the method is highly accurate. This is because the approximation is global basing on a polynomial of a high degree. The eigenvectors contain the wave function values on the mesh and can be used to calculate various expectation values. If this is not enough, once the integral equation has been solved exactly on the mesh, the solution at an arbitrary point may be immediately obtained by interpolation. In conclusion, with the aid of the semi-spectral Chebyshev method the solution of a singular integral equation becomes no more difficult than the solution of a Fredholm equation.

[99]{} A. Deloff, Ann. Phys. (N. Y.), [*in press*]{} R. Omnes, Nuovo Cimento 8 (1958), 316. S.S. Chandrasekarn, [*Radiative Transfer*]{}, Dover, New York, 1960. A. Robinson, [*Wing Theory*]{}, Cambridge U.P. Cambridge, 1958; K. Karamachetti, [*Principles of Ideal Fluid Dynamics*]{}, New York, 1966. N. Brambilla et al., [*Heavy quarkonium physics*]{}, CERN Yellow Report, CERN-2005-005, available from the archive as hep-ph/0412158 E. Eichten and K. Gottfried, Phys. Lett. B 66 (1977) 286; E. Eichten et al., Phys. Rev. D 17 (1978) 3090; E. Eichten et al., Phys. Rev. D 21 (1980) 203 D. Flamm, F. Schöberl and H. Uematsu, Phys. Rev. D 36 (1987); S. B. Elegba and M. A. Rashid, Phys. Rev. D38 (1988) 2911. D. P. Heddle, Yong Rae Kwon and F. Tabakin, Comp. Phys. Commun. 38 (1985) 71; Yong Rae Kwon and F. Tabakin, Phys. Rev. C18 (1978) 932 A. Tang and J. W. Norbury, Phys. Rev. E63 (2001) 06670 T. W. Chiu, J. Phys. A 19 (1986) 2537 D. Eyre and P. Vary, Phys. Rev. D34 (1986) 3467 W. Lucha, H. Rupprecht, and F. F. Schöoberl, Phys. Rev. D45 (1992) 1233; W. Lucha and F. F. Schöberl, Phys. Rev. A56, (1997) 139 L. P. Fulcher, Z. Chen and K. C. Yeong, Phys. Rev. D47, (1993) 4122; L. P. Fulcher, ibid. D50,(1994) 447 L. J. Nickisch, L. Durand, and B. Durand, Phys. Rev. D30 (1984) 660; ibid D30 (1984) 1995(E) L. Durand and A. Gara, J. Math. Phys. 31 (1990) 2237 J. R. Spence and P. Vary, Phys. Rev. D35 (1987) 2191 K. H.  Maung, D. E.  Kahana, J. W.  Norbury, Phys. Rev. D47 (1993) 1183. H. C. Jean, D. Robson and A. G. Wiliams, Phys. Rev D50 (1994) 5873 S. Boukraa and J. L. Basdevant, J. Math. Phys. 30 (1985) 1060 M. Abramowitz and I. Stegun, (Eds.), [*Handbook of Mathematical Functions*]{} (Dover, New York 1972) E. Anderson, Z. Bai, C. Bischof, C. S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney and D. Sorensen, [*LAPACK Users’ Guide*]{}, Third Edition, Society for Industrial and Applied Mathematics, (1999), Philadelphia, PA P. Falkensteiner, H. Grosse, F. Schöberl and P. Hertel, Comp. Phys. Commun. 34 (1985) 287 C,  Quigg and J. L.  Rosner, Phys. Rep. 56 (1979) 167

  -------------- ---------------------- ---------------------- ---------------------- ---------------------- ----------------------
                                                               $\ell=0$                                      
           $N\;$ $\; n=0$               $\; n=1$               $\; n=2$               $\; n=3$               $\; n=4$
            $40$ $4\times 10^{-12}$     $3\times 10^{-10}$     $4\times 10^{ -9}$     $3\times 10^{ -8}$     $2\times 10^{ -7}$
            $60$ $2\times 10^{-13}$     $1\times 10^{-11}$     $2\times 10^{-10}$     $1\times 10^{ -9}$     $6\times 10^{ -9}$
            $80$ $2\times 10^{-14}$     $1\times 10^{-12}$     $2\times 10^{-11}$     $1\times 10^{-10}$     $6\times 10^{-10}$
            $80$ $(4\times 10^{ -5})$   $(9\times 10^{ -5})$   $(2\times 10^{ -4})$   $(4\times 10^{ -4})$   $(6\times 10^{ -4})$
                                                               $\ell=1$                                      
           $N\;$ $\; n=0$               $\; n=1$               $\; n=2$               $\; n=3$               $\; n=4$
            $40$ $8\times 10^{-13}$     $2\times 10^{-13}$     $3\times 10^{ -9}$     $6\times 10^{ -8}$     $5\times 10^{ -7}$
            $60$ $2\times 10^{-14}$     $4\times 10^{-13}$     $3\times 10^{-12}$     $1\times 10^{-11}$     $2\times 10^{-10}$
            $80$ $2\times 10^{-15}$     $4\times 10^{-14}$     $3\times 10^{-13}$     $1\times 10^{-12}$     $2\times 10^{-12}$
            $80$ $(7\times 10^{ -6})$   $(5\times 10^{ -5})$   $(2\times 10^{ -4})$   $(5\times 10^{ -4})$   $(1\times 10^{ -3})$
                                                               $\ell=2$                                      
           $N\;$ $\; n=0$               $\; n=1$               $\; n=2$               $\; n=3$               $\; n=4$
            $40$ $3\times 10^{-12}$     $2\times 10^{-10}$     $2\times 10^{ -7}$     $5\times 10^{ -6}$     $8\times 10^{ -5}$
            $60$ $3\times 10^{-15}$     $6\times 10^{-14}$     $1\times 10^{-12}$     $7\times 10^{-10}$     $2\times 10^{ -8}$
            $80$ $2\times 10^{-16}$     $2\times 10^{-15}$     $4\times 10^{-14}$     $4\times 10^{-13}$     $6\times 10^{-12}$
            $80$ $(1\times 10^{-5})$    $(6\times 10^{-5})$    $(3\times 10^{-4})$    $(7\times 10^{-4})$    $(2\times 10^{-3})$
                                                               $\ell=3$                                      
           $N\;$ $\; n=0$               $\; n=1$               $\; n=2$               $\; n=3$               $\; n=4$
    \[1mm\] $40$ $2\times 10^{ -9}$     $2\times 10^{ -7}$     $9\times 10^{ -6}$     $3\times 10^{ -4}$     $3\times 10^{ -3}$
            $60$ $1\times 10^{-12}$     $3\times 10^{-12}$     $8\times 10^{-11}$     $2\times 10^{ -8}$     $6\times 10^{ -7}$
            $80$ $1\times 10^{-13}$     $2\times 10^{-12}$     $6\times 10^{-12}$     $1\times 10^{-10}$     $5\times 10^{-10}$
            $80$ $(8\times 10^{-5})$    $(3\times 10^{-5})$    $(3\times 10^{-4})$    $(1\times 10^{-3})$    $(1\times 10^{-3})$
  -------------- ---------------------- ---------------------- ---------------------- ---------------------- ----------------------

  : Relative errors on the computed Coulomb binding energies. The corresponding errors appropriate to traditional method based on subtraction are given in parenthesis.

\[table1\]

  ------------- ------------- ------------- ------------- ------------- -------------
                                            $\ell=0$                    
          $N\;$ $\quad n=0$   $\quad n=1$   $\quad n=2$   $\quad n=3$   $\quad n=4$
             50 2.338034      4.087928      5.520416      6.786654      7.943940
            100 2.338099      4.087947      5.520543      6.786702      7.944111
            150 2.338105      4.087949      5.520555      6.786706      7.944127
            200 2.338106      4.087949      5.520558      6.786707      7.944131
            250 2.338107      4.087949      5.520559      6.786708      7.944132
            300 2.338107      4.087949      5.520559      6.786708      7.944133
    [*exact*]{} 2.338107      4.087949      5.520560      6.786708      7.944134
                                            $\ell=1$                    
          $N\;$ $\quad n=0$   $\quad n=1$   $\quad n=2$   $\quad n=3$   $\quad n=4$
             50 3.361254      4.884452      6.207617      7.405649      8.515212
            100 3.361255      4.884452      6.207623      7.405665      8.515234
    [*exact*]{} 3.361254      4.884452      6.207623      7.405665      8.515234
                                            $\ell=2 $                   
          $N\;$ $\quad n=0$   $\quad n=1$   $\quad n=2$   $\quad n=3$   $\quad n=4$
             50 4.248183      5.629693      6.868774      8.009828      9.075383
            100 4.248182      5.629708      6.868883      8.009703      9.077003
    [*exact*]{} 4.248182      5.629708      6.868883      8.009703      9.077003
                                            $\ell=3$                    
          $N\;$ $\quad n=0$   $\quad n=1$   $\quad n=2$   $\quad n=3$   $\quad n=4$
             50 5.050918      6.331874      7.504206      8.593338      9.632163
             80 5.050926      6.332115      7.504646      8.597127      9.627263
    [*exact*]{} 5.050926      6.332115      7.504646      8.597117      9.627267
  ------------- ------------- ------------- ------------- ------------- -------------

  : Binding energy $\epsilon $ for a linear potential.

\[table2\]

                          $m_c$ = 1.37  
  ------------ --------- -------------- ---------
                   $n=0$     $n=1$        $n=2\;$
  $\ell = 0$      3.0869     3.6748        4.1094
                  3.0869     3.6748        4.1093
  $\ell = 1$      3.4988     3.9544        4.3388
                  3.4987     3.9543        4.3388
  $\ell = 2$      3.7868     4.1868        4.5407
                  3.7868     4.1868        4.5407
                          $m_b$ = 4.79  
                   $n=0$     $n=1$        $n=2\;$
  $\ell = 0$      9.4550    10.0105       10.3423
                  9.4547    10.0104       10.3422
  $\ell = 1$      9.9171    10.2582       10.5318
                  9.9170    10.2581       10.5318
  $\ell = 2$     10.1555    10.4385       10.6838
                 10.1554    10.4385       10.6410

  :  Charmonium $(c\bar{c})$ and bottomium $(b\bar{b})$ masses (all entries in GeV) computed from the Coulomb-plus-linear potential [@quigg] $V(r)=-\alpha/r+\beta r$ with the parameters given in (\[q\]). The upper (lower) values result from a calculation conducted in momentum (configuration) space using non-relativistic Schrödinger equation. In all momentum space computations the mesh size was $N=80$.

\[table3\]