Spectral-shift and scattering-equivalent Hamiltonians on a coarse momentum grid

The solution of the scattering problem based on the Lippmann-Schwinger equation requires in many cases a discretization of the spectrum in the continuum which does not respect the unitary equivalence of the S-matrix on the finite grid. We present a new prescription for the calculation of phase shifts based on the shift that is produced in the spectrum of a Chebyshev-angle variable. This is analogous to the energy shift that is produced in the energy levels of a scattering process in a box, when an interaction is introduced. Our formulation holds for any momentum grid and preserves the unitary equivalence of the scattering problem on the finite momentum grid. We illustrate this procedure numerically considering the non-relativistic NN case for $^1S_0$ and $^3S_1$ channels. Our spectral shift formula provides much more accurate results than the previous ones and turns out to be at least as competitive as the standard procedures for calculating phase shifts.


Introduction
Scattering experiments provide usually the most direct phenomenological approach to constrain the corresponding dynamics in hadronic systems. The standard procedure is to determine the scattering observables, which in the case of rotationally invariant interactions reduces to the determination of phase-shifts via a partial wave analysis. The usefulness of the Hamiltonian method becomes more evident when dealing with the fewbody problem, where one expects to determine bound states and resonances of multihadron systems in terms of corresponding potentials. In practice the calculation of phase shifts requires in general solving an scattering integral equation, such as the Lippmann-Schwinger (LS) equation [1] in the non relativistic case where for any given potential V the corresponding S-matrix and hence the phase shifts are determined. While this is a valid perspective for theories where the interaction is known ab initio, ambiguities arise when one tries to infer the potential from scattering information, as it is usually the case in nuclear and hadronic physics. In fact, under an unitary transformation of the Hamiltonian, the corresponding phase-shifts remain invariant. This generates a whole equivalence class of Hamiltonians which are actually compatible with the known scattering information. The notion of equivalent Hamiltonians was introduced by Ekstein in 1960 [2] (see also [3] and [4] for an early review and [5,6] for implications at the relativistic level). This is similar to the equivalence under change of variables in quantum field theory [7,8].
Only in few cases, however, can one provide an analytical solution of the scattering problem, so that one employs often a numerical method which implies a discretization procedure [9,10].
Email addresses: mgomezrocha@ugr.es (María Gómez-Rocha ), earriola@ugr.es (Enrique Ruiz Arriola) In the present work we will be concerned with the discretization in momentum space, since it has been widely used in the past and is the only practicable method for non-local interactions which have been and are commonplace in low and intermediate energy hadronic physics. This requires the introduction of a momentum or energy grid, which has a similar effect to introducing radial boundary conditions [11,12], but is in fact a more general scheme [13] (see [14][15][16][17] for a modern perspective). The interest for this kind of methods for the LS equation started with the work of Haftel and Tabakin [18] where actually the interest was more in providing a method to solve the Bethe-Goldstone equation, which is genuinely non-local, even if the original interaction is local. The computation of elastic-scattering phase shifts via analytic continuation of Fredholm determinants, which are isospectral, constructed using an square integrable basis was introduced in Ref. [19] (for a review see e.g. [20] and references therein). A relevant question is the choice of the particular grid (see e.g. [21] for consideration of adaptive mesh). We will take Gauss-Chebyshev grid whose energy eigenvalue problem corresponds to diagonalizing in a Laguerre basis [22]. We refer to [23] for a comprehensive and self-contained exposition on Chebyshev methods within the present context. To avoid confusion, ours corresponds to a radial one-dimensional momentum grid and not to a three-dimensional momentum grid proposed by the so-called KKR method [24,25] in solid state physics [26] and also called the Lüscher formula in the relativistic case because of more recent popularity within lattice-QCD calculations [27,28].
While the physics of a finite momentum grid is that of the bound states, the continuum limit provides a clear distinction between bound and scattering states. Actually, in a finite momentum grid important properties such as the intertwining properties of the Moller wave operators do not hold [29]. Moreover, the momentum grid solution of the LS equation is not invari-ant under unitary transformations on the finite grid (see for instance [15,17]). The question, of course, is that since one expects that with a sufficiently fine grid the continuum limit will be recovered the number of grid points may be unnecessarily large.
In this letter we provide a method which is in fact rather accurate for a coarse grid, preserves phase-equivalence and in fact is more accurate than the standard numerical solution of the scattering problem based on the (phase-inequivalent) LS equation.
While the problem we address is fairly general, for illustration purposes we consider the toy model separable Gaussian potential discussed previously [30,31] which provides a reasonable description of the NN system in the 1 S 0 and 3 S 1 partial-wave channels at low-momenta, and supports none or one (deuteron) bound state, respectively. The extension to coupled channels (including tensor forces) or resonant and relativistic systems such as ππ or πN scattering requires some modifications and will be discussed elsewhere (See [32] for preliminary results).

Scattering on a finite momentum grid
As already mentioned, scattering occurs in the continuum but numerical approximation schemes discretize it by a finite momentum grid. We review here some well known aspects of both formulations in order to fix our notation.

Continuum formulation
Quite generally we will consider non-relativistic scattering of two particles with masses m 1 and m 2 where H = H 0 + V , H 0 = p 2 /2µ and µ = m 1 m 2 /(m 1 + m 2 ) [33]. Along the paper we will work in unitsh = c = 2µ = 1 so that the free energy is given by E = p 2 . We will assume that the potential is rotationally invariant so that the total Hilbert space can be decomposed as H = ⊕ ∞ l=0 H l and work on the partial wave basis |p, l, m l 1 which is assumed to fulfill the completeness relation in the Hilbert subspace H l Thus, the action of the Hamiltonian on a given state in momentum space is given by where V l (p , p) are the matrix elements and the corresponding Schrödinger equation in momentum space reads, E = k 2 ≥ 0 and a bound-state with (negative) eigenvalue P 2 α = −B α corresponds to a pole in the scattering amplitude at imaginary momentum P α = iγ.
In the continuum the S-matrix is defined as a boundary value problem for E ≥ 0 where we have introduced the T -matrix which satisfies the scattering equation in operator form, where in the second equality we write the exact result. Other (complex) energy values are defined by analytical continuation. This operator satisfies the reflection property (5) and hence the unitarity condition, S(E + iε)S(E + iε) † = 1, follows also from V = V † in Eq. (4). From its definition [S, H 0 ] = 0, and the phaseshift is defined in terms of the eigenvalues of the S-matrix, and introducing the reaction matrix via so that

Discrete formulation
The previous equations can be solved numerically by restricting the Hilbert space H l to a finite N-dimensional space H l,N . On a N-dimensional momentum grid, p 1 < · · · < p N , by implementing a high-momentum ultraviolet (UV) cutoff, p max = Λ, and an infrared (IR) momentum cutoff p min = ∆p [34]. The integration rule becomes and the completeness relation in discretized momentum space reads For instance, the eigenvalue problem on the grid may be formulated as where the matrix representation of the Hamiltonian reads where H nm = H(p m , p m ) and V nm = V (p m , p m ) have been defined.
In practice we will use the Gauss-Chebyshev quadrature method taking the corresponding grid points [35,36], which after re-scaling to the interval [0, Λ] read, where we have introduced the Chebyshev angle, φ n , which will play a crucial role in our considerations bellow. From these definitions we have As it is well known, this grid choice guarantees an exact result for polynomials in p to order M ≤ N, i.e.
Taking matrix elements on the momentum grid of the LS equation in operator form we get where p 2 is the scattering energy. The on-shell limit is obtained by taking p = p l on the grid. As usual we switch to the reaction matrix which on the grid yields the equation for the half-on-shell amplitude where the excluded sum embodies the principal value prescription of the continuum version in the limit ε → 0. This equation can be solved by inversion for any grid point p n and thus we may obtain the phase-shifts where the supper-script LS denotes that these phase-shifts are obtained from the solution of the LS equation on the grid. Of course, the limit N → ∞ should be understood in the end. As mentioned, one drawback of the LS formulation is the fact that if we undertake a unitary transformation on the grid, U nm of the Hamiltonian the, the T nm (p) in Eq. (20) still transforms as its continuum counterpart (for any ε) but the phase-shifts given by Eq. (22) do not remain invariant due to the principal value prescription 2 . Hence the LS phase-shifts are not isospectral. 2 This has been explicitly shown within the context of the similarity renormalization group [15,17].

Spectral shifts: the transition from the discrete to the continuum
In the previous section we have described how the scattering equations can be solved by discretizing the spectrum, thus violating the isospectrality of the phase-shifts. In this section we discuss three alternative spectral shifts: the momentum shift, the energy shift and the Chebyshev-angle shift which enjoy the isospectrality of the phases.

Momentum-shift prescription: scattering in a spherical box
The momentum-shift prescription proposed by Fukuda and Newton long ago [11,12] is the simplest and assumes that the scattering process takes place in configuration space in a large spherical box of radius R. For a potential with a finite range a the reduced wave function u l (r) has the following asymptotic behavior for R ≥ r a and must vanish for r = R, so that u l (R) = 0, which implies that In a finite momentum grid of N points the equation yields N eigenfunctions, and thus Eq. (24) holds for every p n . In absence of interactions, δ l = 0, and one has with n = 1, . . . , N, and ∆p n ≡ p n+1 − p n = π R is the separation of the grid points. Representing now the distorted momentum by P n , and using Eq. (25) to replace nπ + lπ/2, one can write which, using that ∆p n = π/R becomes This is the momentum-shift formula which strictly speaking holds for an uniformly distributed momentum grid. This is equivalent to a trapezoidal rule quadrature which is generally a poor integration method. In Section 4 we will consider Eq. (27) for the Chebyshev grid quadrature method, namely 3.2. Scattering in an energy-equidistant discretized spectrum An alternative approach was proposed simultaneously de De-Witt [13]. Let us present DeWitt's argument in a slightly different way so that our points can be easily formulated. For the sake of clarity and the benefit of the reader we try to be pedagogical here since we found some parts hard to follow. Let us consider the eigenvalue problems with H = H 0 + V . Our notation is such that E n → E (0) n when V → 0, and as according to the Landau-von Neumann theorem, there is no crossing between non-degenerate levels [37].
The cumulative number associated to the discretized Hamiltonian H reads where we have introduced the trace TrA = ∑ n ψ n |A|ψ n . For a discrete spectrum this function represents a staircase, with unit jumps at any eigenenergy, E n . In order to proceed further it is important to separate the states into positive and negative energy states. In the continuum limit the negative energy states will become bound states whereas the positive energy states will become scattering states. Starting from the free Hamiltonian where all energies are positive, E n > 0, with a gradually increasing attractive interaction some energy levels may drive into the negative energy spectrum. From a variational point of view the discretized energy provides an upper bound of the true spectrum, since the discretization procedure may be viewed as a restriction on the physical Hilbert space, E N n ≥ E ∞ n and thus the net effect of the finite grid is repulsive. In what follows we will assume that the grid is fine enough so that no positive energy level will cross zero.
The step function in Eq. (31) can be regularized as proposed in Ref. [38,39], namely introducing a small imaginary energy iε → i0 + , as follows where for a general complex number z = ρe iθ with −π ≤ θ ≤ π we define the logarithm where the branch cut runs along the negative real axis. Deriving with respect to x we also get where it is important in what follows to keep a finite ε and to take the continuum limit N → ∞ with ∆e n ≡ E n → 0. Then we have, the regularized cumulative number where the identity log DetA = Tr log A has been used. Using a similar equation for H 0 and its corresponding cumulative number N 0 we get for the difference, where G 0 (E − iε) = (E − iε − H 0 ) −1 is the resolvent of the free Hamiltonian and the cyclic property of the trace can be used to allocate V to the left or to the right of G 0 . The discontinuity of this function is defined as Direct application of the eigenvalues yields, In order to carry out the sum we use the trigonometric identity, so that we get In this formula the continuum limit ∆e → 0 has to be taken before the limit ε → 0. In the limit ∆e/ε 1 we change the summation into an integral, ∑ n → dn. One crucial aspect in DeWitt's formulation is the explicit use of a uniform energy spectrum, so that one takes E (0) n = n∆e and thus E n = n∆e + ∆E n where ∆E n is the energy shift. Clearly, in the continuum limit ∆e → 0 the energy shift vanishes ∆E n , but the ratio ∆E n /∆e remains finite. Defining the change of variables t = n∆e/ε, in the limit ∆E n /ε → 0, one gets, after shifting the integration variable, In order to connect this sum with scattering information we take into account Eq. (38) and note that the r.h.s of Eq. (4) can then be written as where we have used that the S-matrix is energy diagonal in the continuum, i.e. S(E)ϕ α (E ) = 0 for E E. Merging this result and the previous Eq.(45) we finally get where the energy dependence is in the energy shift. We stress that this formula holds for an equidistant energy grid, which would correspond to a trapezoidal rule in energy. If we write it in terms of momentum variables for the purpose of applying the Gauss-Chebyschev grid defined by Eqs. (13)-(16) we get

Chebyshev-shift
In order to profit from both the use of the Gauss-Chebyschev grid and the previous DeWitt's formula, we introduce the angle φ given by so that the Gauss-Chebyschev grid provides an equidistant angle. Thus, if we use φ as the independent variable we may apply DeWitt's argument mutatis mutandis for a cumulative number N(φ ). Using the analogy between the energy levels of scattering states in a box and the discretization given by a finite grid, and observing that the equidistance happens in the argument of the cosinus function, we prescribe the following formula based on the shift of such an angle, and write: where φ n = π N n − 1 2 , dφ n = π N , and the "distorted" angles Φ n are calculated from Eqs. (13)- (16) replacing p n by P n . Thus, This is the main result of this paper. One might think that using a similar change of variables in the LS equation might alter its convergence properties, since discretization and reparametrization are generally non commutative operations. Note however that, from Eq. (13) one has On the Chebyshev grid one has d p n = w n if dφ n = π/N, so that this change of variables does not modify the original discretized equation in momentum space.

Bound state modifications
The occurrence of a bound state modifies the formulas in the case of the energy shift, since direct application of the differences violates the well known Levinson's theorem [40], which in the continuum becomes with n l the number of bound states. In the discrete case, there appears a discrete momentum scale which requires some reordering of the states [15,17]. In the case of the energy shift it becomes whereP 2 n BS = (P 2 n BS +1 + P 2 n BS −1 )/2. Note that in this prescription only the eigenvalues P 2 n corresponding to momenta p n < p n BS are shifted one position to the left. The cases for the p-shift and φ −shift are similar.

Numerical results
We come to our numerical results and analyze how the scattering phase-shifts obtained from the energy-shift, momentum shift and φ -shift formulas compare with the results obtained from the standard method based on the LS equation. Along the text, we may use the abbreviations p-shift, E-shift, and φ -shift in order to refer to the momentum-shift, energy-shift and angle-shift prescriptions given by Eq. (28), (49) and (52), respectively.
We will see that the φ -shift method prescription is the best and the only one that reproduces almost exactly solution calculated in the continuum, even for a coarse grid. In fact, as we shall show in our numerical study, the method gives reliable predictions even for a grid with a very small number of points. The generalization to any momentum grid amounts to finding the variable that is uniformly distributed along the momentum grid.

Separable models
In our numerical comparisons we use a separable model potential with the structure v l (p , p) = C l g l (p )g l (p) , where C l is positive or negative for repulsive and attractive interactions respectively. For this form of potential the LS equation, Eq. (6) is solved by the ansatz T l (p , p, E) = g l (p )g l (p)T l (E) , which inserted in Eqs. (7)-(8) yields In practice we take the toy Gaussian potential g α (p) = e −p 2 /L 2 α proposed in [15,17] for NN scattering in the Taking these values, we may then proceed to check the phase-shift determined by the p-, Eand φ -shift formulas, which only generates them on grid points.

4.2.
Dependence on the momentum grid and comparison with the standard method In this section we study numerically how our φ -shift results, calculated in a finite momentum grid, differ from the exact solution in the continuum (green, smooth line in all figures) and confront these results with the standard method based on the LS equation. Figure 1 shows the phase shifts obtained by solving the LS equation, Eq. (22) (in orange), and by the φ -shift prescription Eq. (52) (blue dots). Both of them are compared with the numerical fit (smooth, green line) which represents the exact solution. Each column of Fig. 1 corresponds to the same calculation, but using a different number of grid points, N =10, 20, 50, respectively. It is impressive how the φ -shift formula reproduces exactly the exact solution in all cases, even in the grid with least points, N =10. Conversely, we observe how the LS method converges to the continuum as the number of points increases and a comparatively larger number of grid points is required to reproduce the exact solution. The LS results are worse in the case of the 1 S 0 wave, where the phase shifts undergo larger changes is a relatively small interval of momentum and thus the number of points describing the curve becomes more important.

Comparison with different spectral-shift prescritpions
In this section we compare the φ -shift prescription with the p-shift and E-shift ones. Figure 2 shows, for both channels, the resulting phase shifts obtained from the different formulas, as indicated by a label in the corner. All of them are compared with the exact solution.
In all cases, we may represent the results as a function of the interacting momentum P n , or as a function of the free momentum p n . Every graphic in Figure 2 shows both curves and Figure 3 shows a selected interval of the 1 S 0 channel, where the difference between lines is more visible. Here δ (P n ) is represented by a darker line with round markers, while δ (p n ) is represented by a lighter line with square markers. Of course, the phase shifts have the same values but are horizontally displaced from each other by the momentum shift, as it can be observed in the figures. As it turns out, the formulation as a function of the interacting momentum lies closer to the exact solution in all cases. This appears reasonable if one notes that in the case of the p-shift formula, the phase shift must be a function of the interacting momentum by construction. Although there is nothing in DeWitt's argument [13] that suggests that the distorted momentum should be used as the independent variable, we note that in the relativistic case and for very light masses the E-shift prescription converges to the p-shift one 3 . This suggests to 3 In order to check this, one can replace the energy by the relativistic formula E = p 2 + m 2 in Eq. (48), and take the limit m → ∞.
consider the interacting momentum as the independent variable in all cases.
It is remarkable that in the φ -shift case (first column of Figures 2 and 3, in blue), both lines, the one depicted as a function of p n and the one as a function of P n , precisely overlap the exact solution and there is no difference among both criteria. Indeed, one can barely see the green line in this case, and the darker blue line is totally covered by the lighter-blue line.
The fact that the E-shift or the p-shift prescriptions produce worse results than the φ -shift one is reasonable. These prescriptions assume equal-distance separation of energy levels and momentum levels, respectively, while in our Gauss-Chebyshev grid, this separation occurs in the Chebyshev angle. The E-shift and p-shift prescriptions appear thus to be a (non exact) but approximate formula. Nevertheless, both of them still turn out to be a considerably good approximation, since the quality of the results for the studied N = 20 case are by far better than those obtained through the standard LS equation.

Conclusions
In this letter we have presented a new method which solves the scattering problem by diagonalizing the corresponding Hamiltonian in a momentum grid. This guarantees that the phase-shifts are invariant under unitary transformations on the finite grid, unlike the usual solutions based on the Lippmann-Schwinger equation.
We presented and studied the predictive power of the momentum-shift and energy-shift method for calculating phase shifts. We have proposed a new prescription based on an argument that holds for almost any momentum grid. The prescription requires to find the variable that holds an equidistant space between points along the momentum grid. In our case, the chosen grid is a Gauss-Chebyshev quadrature and the equal spacing occurs in the angle φ = π N (n − 1/2). Having identified this quantity, we follow DeWitt's reasoning [13] and calculate phase-shift from the φ -shift produced by the interaction. We observe that this prescription yields remarkable good results, even for a grid with a small number of points. We have observed, furthermore, that, in contrast to what is formulated by DeWitt, it is better to consider the phase shifts as a function of the distorted momentum -following the momentum-shift prescription [11] -, instead of as a function of the free momentum.
Theses result suggest that the new φ -shift method offers a reliable tool for numerical calculations of phase shifts in a discretized momentum grid which turns out to surpass in precision other conventional formulas. Extensions to coupled channels and relativistic systems are straightforward and will be analyzed in detail elsewhere.