{'fraction_non_alphanumeric': 0.03655709445402128, 'fraction_numerical': 0.010461867681150806, 'mean_word_length': 4.173431734317343, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Nystrom plus Correction Method for Solving Bound State Equations in Momentum Space

A new method is presented for solving the momentum-space Schrodinger equation with a linear potential. The Lande-subtracted momentum space integral equation can be transformed into a matrix equation by the Nystrom method. The method produces only approximate eigenvalues in the cases of singular potentials such as the linear potential. The eigenvalues generated by the Nystrom method can be improved by calculating the numerical errors and adding the appropriate corrections. The end results are more accurate eigenvalues than those generated by the basis function method. The method is also shown to work for a relativistic equation such as the Thompson equation.


I. INTRODUCTION
The momentum space Schrödinger equation has a singular kernel for both the Coulomb and linear potentials. The Coulomb singularity is removed with the Lande subtraction method [1,2]. Previous work [3][4][5] showed how to remove the singularity from the linear potential using a subtraction method with basis functions. A problem with this method is that one must guess a suitable set of basis functions in advance. In this paper, we show that Nystrom method [7] can solve the same problem more simply amd accurately. We begin with a review of the basis function method. Then we introduce the Nystrom method and apply it to the s-state momentum space Schrödinger equation with a linear potential. We use our new numerical results to show that the Nystrom plus correction method is more accurate than the basis function method. At the end, we generalize the Nystrom method to higher angular momentum quantum numbers (l > 0).

II. THE BASIS FUNCTION METHOD
We begin this paper with a discussion of the basis function method to give the proper theoretical motivation. We shall use the simplest momentum space Schrödinger equation to illustrate the principles of the numerical methods, which is the s-state equation.
The momentum space Schrödinger equation is related to an integral equation of the Fredholm type K(p, p ′ )φ(p ′ )dp ′ = λφ(p). (1) Suppose that the wavefunction φ can be expanded in a set of basis functions {g i }, such that where C i are constant coefficients. Substitute Eq. [2] into Eq. [1] to obtain Now multiply both sides of Eq. [3] with g j (p) and integrate over p to symmetrize the equation over i and j, N i=1 C i K(p, p ′ ) g i (p ′ )g j (p) dp ′ dp and the result is a matrix equation, where C i is the eigenvector and λ the eigenvalue. The indices i and j correspond to the quadrature points p and p ′ . N represents the number of mesh points. In the case of the momentum space Schrödinger equation with a Coulomb or linear potential, the kernel A is singular. A simple example is the momentum space Schrödinger equation with a linear potential in the S-state [3,5], where y = (p 2 + p ′2 )/2p ′ p, and Lande subtraction [3,5,6] involves subtracting a zero term from Eq. [6] such that Using Eqs. [7,8], the integral in Eq. [10] for p > 0 in the limit of y → 1 can be shown to The order of the limits in Eq. [11] is important. The reverse order will lead to the nonsensical result Q ′ 0 (y) dp ′ = 0. Next, in the limit of p, p ′ → 0, (p + p ′ ) 2 = (p − p) 2 . By substituting this equality into Eqs. [7,8], it can be shown again that the integral in Eq. [10] where Q ′ 0 (y) = 1/(1 − y 2 ). As mentioned before, φ is expanded in basis functions, followed by integrating Eq. [12] over p to generate a matrix equation. The basis functions used in previous publications [3,5] are and respectively, where M is the maximum number of basis functions used. M has a maximum because the code crashes when too many basis functions are used. The basis functions g A i (p) and g B i (p) have M = 18. The singularity of the kernel is a major challenge in solving the integral equation with a linear potential. It was shown [3] that a simple pole remains even after subtraction. The role of the basis functions is to make possible the evaluation of the Cauchy principal value of the subtracted integral using the Sloan method [8]. To illustrate the Sloan method, we suppose that f (x) has a simple pole such that where g(x) is regular. The Cauchy principal value of the subtracted integral of f (x) can be evaluated if the range of integration is symmetric. For example, the numerical integration yields the Cauchy principal value because the point x = 0 is skipped when quadrature points are generated in the symmetric interval (−1, 1). The subtraction term has zero contribution The purpose of this term is to justify the existence of the Cauchy principle value and to reduce numerical errors. In order to apply the Sloan method to Eq. [12], the integration variable is transformed from p to x such that x is centered at zero and its range is symmetric.
In the case of the Coulomb potential, the kernel has a logarithmic singularity, which is completely removed by Lande subtraction [1,3] because no simple pole remains after the subtraction.
The key to the success of the basis function method is the availability of a suitable set of basis functions for a particular problem. Unfortunately there is no a priori reason why the same set of basis functions will work in every situation. For this reason, it may be advantageous to have a method (such as the Nystrom method) that does not depend on a choice of the basis functions.

III. THE NYSTROM METHOD
In general, an integral equation of the Fredholm type can be rewritten as a matrix equation as where K ij is the kernel and i and j are now indices corresponding to p and p ′ . Instead of integrating over p from 0 to ∞, we integrate over x from −1 to 1. Transform x i to p i by the transformation The mesh points x i and the weights wt i are generated by the gaussian quadrature rule using the routine gauleg from Numerical Recipes [7]. In order to integrate along Changing the dummy variable inside the integral and substituting the differentiation of Eq.
[24] can now be written as a matrix equation, The left hand side of Eq.
[25] is the kernel times the eigenvector and the right hand side is the eigenvalue times the eigenvector. The sum over k is independent of the eigenvector, which is just a scalar. The terms on the left that have only one index i belong to the diagonal elements K ii . The terms with mixed indices i and j make up the off-diagonal elements, K ij .
More explicitly, the matrix elements of the kernel are So far the kernel K is asymmetric under the interchange of i and j. We can improve the stability and the efficiency of the numerical solutions by symmetrizing Eq.
[25]. We do so by multiplying the equation with p 2 i (1 + p 2 i ). It will change the original matrix equation to an equivalent matrix equation where C is a diagonal matrix and K ′ = C · K. If C is positive definite, meaning then C can be Cholesky-decomposed as where L is a unique lower triangular matrix. The reason for performing Cholesky decomposition is that the new matrix is real symmetric and yields the same eigenvalues as Eq.
[29]. In the case of C is guaranteed to be positive definite. After symmetrization, Eq.
[26] does not change The eigenvalues of K ′′ can be calculated by using standard packages such as EISPAK.
In this paper, we use the tred2 and tqli routines in Numerical Recipes [7].

IV. THE CORRECTION METHOD
Maung, Norbury and Kahana [3,5] have shown that the subtraction method does not completely remove the singularity at y = 1. There is a residual simple pole term remaining after the subtraction. The basis function method evaluates the Cauchy principal value by the Sloan method as described in Section 2. The Sloan method eliminates the simple pole term by integrating symmetrically around the singularity. Symmetrical integration involves splitting the range of integration is into two intervals, The singularity at p = p ′ is contained in the first term on the right hand side of Eq. [35] which is assigned a symmetric transformation rule (dp ′ /dx) 1 . The second term generally has a different transformation rule (dp ′ /dx) 2 because it is mapping between two different sets, Notice that the division of the range of integration moves with p. If two transformation rules are used with a moving division, each row (column) of the kernel has a different way But the eigenvector φ(p) must be mapped to φ i in a unique way.
This mismatch between the mappings of the kernel and the eigenvector does not affect the basis function method (see Eq. [2.24] of reference [3]) because the eigenvector C i is an N-tuple of the coefficients of the basis function expansion of the wavefunction φ(p) and is independent of the transformation rules. In the case of the Nystrom method, the problem is real, at least for the range of integration that we are interested in. Therefore we cannot evaluate the Cauchy principal value by symmetric integration in the Nystrom method. In other words, a new method is needed to treat the errors arising from the simple pole term.
So far the error term Eq.
[34] is not included in the Nystrom kernel in our derivation and is contributing to the errors of the eigenvalues. Since the error term Eq.
[34] involves dφ/dp, we associate it with the error of the wavefunction ∆φ = ∆p dφ dp where the mesh size ∆p has an N −1 dependence. This fact leads to an estimate of the N dependence of the error of the eigenvalue, ∆E. Let the approximate eigenvector to be φ ′ and the approximate eigenvalue E ′ . It is reasonable to say that an approximate kernel K acting on an approximate eigenvector φ ′ yields an approximate eigenvalue E ′ as in It is easy to see that It is safe to assume that (K − E) dφ/dp << 1. φ −1 can be interpreted as the normalization.
The product of all of the pseudo-constants is labelled as the coefficient ǫ. The approximate eigenvalue E ′ produced in the background of Eq.
[34] is related to the true eigenvalue E by where n is the principal quantum number, ǫ a constant and f n (N) is a function approximately equals to N −1 . In general, f n (N) varies slightly depending on the type of integral equation and the potential. As a first order approximation, assume that The exponent of Eq.
[43] is a first order Taylor series expansion of some negative unity function around n = 1. The constant α is always taken to be small. More particularly, choose an α such that the variance of E n and ǫ and χ 2 are minimized in the linear fit.
Finally the refinement of an eigenvalue involves generating a set of E ′ n for various N by the Nystrom method and then extrapolating E n by a χ 2 linear fit in the graph of E n versus f n (N). In the case of Eq. [12], α = 0.004 is an optimal choice. The numerical results are explained in Section 6.
The order of the Nystrom algorithm is derived from those of tred2 and tqli, which is

V. EXACT S-STATE SOLUTION
The eigenvalue of Eq. [12] can be solved exactly in configuation space. We shall use the analytic results to check our numerical results. The non-relativistic Schrödinger equation can be written as Let S ≡ r R, then Eq.
[44] can be simplified as Define a new variable such that Eq.
[44] can be transformed as which is the Airy equation. The solution which satisfies the boundary condition S → 0 as x → ∞ is the Airy function Ai(x). It is easy to show that the eigen-energy formula is where x n is the n-th zero of the Airy function counting from x = 0 along −x. In reference [5], the values λ L = 5 and µ = 0.75 are used. In this case, the eigen-energy formula is

VI. NUMERICAL RESULTS FOR THE S-STATE
The accuracy of the Nystrom plus correction method is sensitive to the range of N. In where µ is the reduced mass and m is the mass of each of the two equal mass elementary particles. The numerical results obtained using the Thompson equation is compared against those using the non-relativistic Schrödinger equation in Table [II] calculating to 2 decimal places. Our new results are exactly the same as the previous results obtained in reference [3] that uses basis functions g A i (p) from Eq. [13].

VII. L = 0 KERNELS
The l = 0 kernels for the linear and Coulomb potentials contain the Legendre function of the second kind Q l (y) and its derivative respectively. There are several mathematical issues that need to be addressed before constructing the l = 0 kernels. First of all, the definition of y ≡ is easily seen to yield y > 1 for p, p ′ > 0. In reference [3], Maung et al. use the Legendre identity Q l (y) = P l (y)Q 0 (y) − w l−1 (y), which is valid for −1 ≤ y ≤ 1 [9,10] but can be extended to y > 1 by analytic continuation [11]. Q ′ l (y) is easily obtained by straightforward differentiation. The derivative of Legendre polynomial can be calculated from one of the recurrence formulas, dP l (y) dy = y dP l−1 (y) dy which can be computed numerically by a recursive call. The Legendre function can be generated by modifying the routine plgndr in Numerical Recipes [7] to allow y > 1. The accuracy of Eq.
[52] and its derivatives are generally sufficient. Slightly more accurate results can be obtained by the explicit evaluation of the Neumann integral, with derivative The first few Q l (y) are Q ′ l (y) can be obtained by the direct differentiation of Q l (y), such that As y → ∞, it is easily seen that Q 0 (y) = Q ′ 0 (y) → 0. This limit is true for all Q l (y) and Q ′ l (y) from applying the L'Hopital rule. Unfortunately straightforward numerical calculation of Q l (y) and Q ′ l (y) by using Eqs. [56-67] leads to serious numerical errors as y → ∞. At the same time, it is observed that the numerical integration of Eq. [54,55] are reasonably accurate in the same regime. Therefore the two representations are combined to minimize numerical error by using the Neumann integrals for y > y 0 and the explicit formulas for y ≤ y 0 . Our codes use y 0 = 50.
The subtracted momentum space NRSE with a linear potential is given in reference [3], which can be simplified as The matrix elements of a symmetric kernel for arbitrary l are Despite our method to control numerical noise, numerical errors still manifest themselves in the form of spurious large negative eigenvalues for l ≥ 8. Fortunately the rest of the positive eigenvalues are accurate. Some sample eigenvalues for 0 ≤ l ≤ 5 are shown in Table [III] which also compares the eigenvalues generated by both the p-space and r-space codes. The r-space eigenvalues are calculated by solving NRSE using the relaxation method [13].
The Lande-subtraced momentum space NRSE equation with a Coulomb potential is also given in reference [3] and is simplified as The kernel of a Coulomb potential can be symmetrized the same way as that of a linear potential. The matrix elements are The correction method which we have developed for the linear potential cannot be used in the Coulomb case. The only available technique of refining the eigenvalues of a Coulomb potential is by the way of increasing the number of mesh steps N. Some sample eigenvalues are shown in Table [IV]. Since both the linear and Coulomb potentials can be symmetrized using the same formalism, we can easily splice the two kernels together to calculate the eigenvalues of the Cornell (linear plus Coulomb) potential It is not surprising that the correction method derived for the linear potential may also work for the Cornell potential because we expect that the error of the Cornell potential is dominated by the error of the linear potential term. But it is a surprise that the correction method works more accurately with the Cornell potential than the linear potential as evidenced by vanishingly small variance and χ 2 .

VIII. CONCLUSION
The basis function method requires a priori knowledge of the eigenfunctions in order to pick out an appropriate set of basis functions. The advantage of the Nystrom method is that no such prior knowledge of the eigenfunctions is needed. The kernel constructed by the Nystrom method is also much simpler than that by the basis function method. The eigenfunctions can be generated by the same Nystrom routines that compute the eigenvalues.
The Nystrom plus correction is more accurate than the basis function method in the cases studied in this paper. In other words, the new method has all of the advantages-elegance, accuracy and versatility. In addition, the kernel of the relativisitc and non-relativistic equation of motion with the Coulomb and linear potential can be symmetrized in exactly the same manner. It allows the calculation of the eigenvalues of a Cornell potential readily. Since the Nystrom method can be generalized for higher l's, we can use it to calculate the Regge trajectories. Since the Thompson equation which we have solved is a 2-body equation, we can use it to analyze the experimental meson Regge trajectories [14]. This will be pursued in later work.