Complete solution of the Marchenko equation for a simple model system

An example of full solution of the inverse scattering problem on the half line (from 0 to ∞) is presented. For this purpose, a simple analytically solvable model system (Morse potential) is used, which is expected to be a reasonable approximation to a real potential. First one calculates all spectral characteristics for the fixed model system. This way one gets all the necessary input data (otherwise unobtainable) to implement powerful methods of the inverse scattering theory. In this paper, the multi-step procedure to solve the Marchenko integral equation is described in full detail. Several important analytic properties of the Morse potential are unveiled. For example, a simple analytic algorithm to calculate the phase shift is derived.


INTRODUCTION
Strict mathematical criteria for the unique solution of the one-dimensional inverse scattering problem have been formulated long ago, thanks to the outstanding contributions by many researchers [1][2][3][4][5][6][7][8][9][10][11] (see [12] for a historical overview).However, in spite of the perfect consistency and mathematical beauty of the concept, these rigorous criteria have rarely been used in practice, because the necessary input data are very difficult to obtain.Indeed, one has to know the full energy spectrum of the bound states E n < 0 (n = 0, 1, 2, . . ., N − 1) and the full energy dependence (from 0 to ∞) of the phase shift δ (E) for the scattering states (E > 0).Even if the mentioned obligatory requirements were fulfilled (which is never the case), this would only be sufficient to construct an N-parameter family of phase-equivalent and isospectral potentials.One has to somehow fix N additional parameters, the so-called norming constants, in order to ascertain the potential uniquely.
On the other hand, the shape of a real potential is not arbitrary, but often it looks like the curve in Fig. 1.In this figure, a comparison is made between an ab initio potential for Ne 2 molecule [13] and its Morse approximant, calculated according to the formula [14] V (r) taking D = 42.153K [13,15] and R e = 3.0895 Å [13].The Morse curve in Fig. 1 was calculated with r 0 ≡ √ h2 /(2mD) = 0.2388 Å, α = 0.4879/r 0 , and ∆r ≡ r − R e .It has only one bound state: E 0 /D = −0.5716,and practically the same value has been determined for the single discrete level of the Ne 2 molecule, both theoretically [13] and experimentally [15].Figure 1 was presented for illustrative purposes, to continue the above discussion.Namely, in principle, the unknown norming constants for a real potential can be treated as variational parameters.One may try to guess their correct values by fitting to Eq. ( 1).Moreover, one might even think about using this simple model potential to examine different solution schemes of the inverse scattering theory, and this is the main idea developed in this paper.The Morse potential has many unique analytic properties, which makes it especially suitable for modelling.In particular, for this model system one can ascertain all the necessary spectral characteristics to any desired accuracy.Thereafter, one can use them as input data, instead of getting these data from real experiments.
The idea to invert an inverse problem is, of course, tautological but not at all meaningless.Nowadays people tend to consider the inverse scattering a 'pedagogical' problem of little scientific interest.However, such an ignorant view is not based on the real understanding of the complexity of the problem, but is a mere belief.The implementation of the methods of the inverse scattering theory is not at all a trivial task.On the contrary, this is a complex and computationally very demanding multi-step procedure that has to be performed with utmost accuracy.In this paper, we are going to describe the full solution of the Marchenko integral equation.To this end, as already mentioned, a very simple model system has been chosen, to make the procedure understandable for a wide audience of potential readers.This can be considered as the first step towards the real implementation of the method.
Performing an applicability test of the Marchenko method is not the only issue addressed in this paper.In the author's opinion, the model system itself deserves special attention.Useful properties of the Morse potential, as well as the solution to the related Schrödinger equation already given by Philip Morse himself [14], are well known to physicists.Since the potential is shape-invariant, this solution can even be found by pure algebraic means, using the Gendenshtein's recipe [16].However, in this context we are talking about the physically correct linear combination of the two special solutions to the Schrödinger equation.There are lots of systems which can be analysed in exactly the same way (the simplest one is the harmonic oscillator), but within this pattern the feature which makes the Morse potential really unique is overlooked.Namely, for the Morse potential the two linearly independent solutions to the Schrödinger equation can always be easily ascertained analytically (even if D < 0!).This is not the case for a harmonic oscillator or any other popular model system.
Consequently, the approach based on Eq. ( 1) is not just an option among many others, but it is often the most appropriate choice for modelling a real quantum system.Moreover, one may construct a more reliable model consisting of several Morse-type components, still preserving the exact solubility of the Schrödinger equation.Such technically rather complex developments are not analysed in this paper (see, e.g., [17] and [18] for more details), because here we concentrate on solving the inverse not the direct Schrödinger problem.For this reason, the model has been chosen to be as simple as possible.
Nevertheless, as a matter of fact, the scattering properties of the Morse potential have not found much attention so far, and there still is some work to do.In Section 3, a simple analytic formula for calculating the phase shift is derived.Unexpectedly, in addition to this specific result, there is a pure mathematical outcome which is directly related to Eq. (1).Namely, a very accurate algorithm for evaluating the Riemann-Siegel function (see [19]) was obtained, as a kind of bonus for calculating the phase shift.
The remaining parts of the paper are organized as follows.The details of calculating the kernel of the Marchenko equation are described in Section 4 and the solution method itself is under examination in Section 5. Finally, Section 6 concludes the work.

MARCHENKO INTEGRAL EQUATION
The basis for this paper is the integral equation derived by Marchenko [8].The Marchenko method is preferred, because the kernel of Eq. ( 2) is directly related to the main spectral characteristics of the system, while an alternative approach, the Gelfand-Levitan method [5], requires an additional calculation of the Jost function.
The kernel of Eq. ( 2) reads where C ≡ h2 /(2m) and s n are the norming constants for the Jost solution of the related Schrödinger equation: so that f (iγ n , r) → exp(−γ n r) as r → ∞, and If one is able to solve Eq. ( 2), the potential can be easily determined: Equation ( 2) and its kernel (3) are a bit different than originally given by Marchenko [8], but there is no contradiction.Indeed, the original Marchenko equation (see, e.g., [20,21]) can be rewritten as Consequently, taking A 0 (x) ≡ −B 0 (x), one comes to Eqs (2)- (3).Note that the second term in Eq. ( 3) begins with a minus sign, and this is an important nuance, because this sign is incorrectly given as 'plus' in several highly appreciated overviews (see, e.g., [12], p. 79).Now, let us specify the problem and the scheme for its solution.First, one calculates the spectral characteristics for a simple model potential described by Eq. (1).Then one determines the kernel of the Marchenko equation according to Eq. ( 3).Finally, one solves Eq. ( 2) and calculates the related potential according to Eq. ( 6).If the procedure is successful, the initial potential is recovered.Throughout this paper, dimensionless units for the energy and radial coordinate will be used, taking For a Morse potential, the norming constants s n can be easily ascertained analytically.Let us fix, for example, α = 2/3.Then the system has only one bound state: where a ≡ √ D/C/α.Consequently, in our unit system, a = 1/α and γ 0 = 2/3 = α.Therefore, this particular case is especially suitable for illustrative purposes and only this case will be examined further on.(The same approach, with minor adjustments concerning the different values for α and γ 0 , can be applied to the model potential shown in Fig. 1.)To conclusively fix the model, let us assume that the minimum of the potential is located at R e = 2.5r 0 .
Strictly speaking, to exclude an extra bound state near E = 0, we have to set the constraint V (0) = ∞, i.e., the potential jumps to infinity (not to zero!) as r → 0. The problem is that if one takes α = 2/3, then a = 3/2 and, according to Eq. ( 8), E 1 = 0.However, this extreme case is excluded because the coordinate here ranges from 0 to ∞, not from −∞ to ∞ as for a classical Morse oscillator.Therefore, as can be proved by elementary perturbational arguments, the possible additional level is shifted slightly above (not below) the dissociation limit, thus becoming a scattering state.Another nuance is that, in principle, one should form a linear combination of both solutions of the Schrödinger equation, to fix the correct regular solution.Again, it can be shown (see, e.g., [22]) that the second special solution is of next to no importance, so that the physical solution reads as in the 'normal' Morse case.Here N 0 is the related norming constant and a new dimensionless coordinate y ≡ 2a exp [−α (r − R e )] was introduced.The parameter N 0 as well as the norming constant for the Jost solution (4) can be easily determined: The second term in the denominator is negligible, so that, to a high accuracy, s 2 0 = y 2 0 α = 6 exp(10/3).This is the value which we are trying to recover.However, to clarify the sign problem described above, the quantity s 2 0 will be treated as a variational (not definitely positive!)parameter.We will see that its expected value is indeed the correct value.

PHASE SHIFT FOR THE MORSE POTENTIAL
Let us see how to determine the full energy dependence of the phase shift, subjected to a strong constraint (N = 1 in our case), according to the Levinson theorem [1,2].The analytic background of the problem has been described elsewhere [18] (see Eqs ( 23)-( 29) there).There are several formally equivalent approaches to building the regular solution of the Schrödinger equation for a scattering state.For example, we may proceed from Eq. ( 25) of [18].Then, after some analytical work, we get the following special solutions: with and (N = 1 in our case), Γ being the gamma function.
From the regularity requirement

one gets the ratio of the coefficients
On the other hand, since S(a, iβ ; y) → 1 as r → ∞ [18], one gets a formula for the phase shift: and consequently, a general expression for the wave function Ψ = A 1 Ψ 1 + A 2 Ψ 2 (for simplicity we drop the arguments of the S-functions) The last formula is valid for any r, and also in the limit r → ∞.From this one immediately concludes that and most importantly, This is the simple formula mentioned in Section 1. Combined with Eq. ( 13), it enables us to easily ascertain the phase shift for any scattering state and thus, their full dependence on energy (or wave number).The result for the specified Morse potential is shown in Fig. 2 (main graph).Note that the k-scale there is logarithmic, so that the figure involves 18 orders of magnitude in the energy space.Full agreement with the Levinson theorem (11) can be seen as well.Figure 3 demonstrates the behaviour of the scattering wave functions as E → 0. In general, the wave functions can be calculated using Eqs (13), (18), and (19).Note, however, that the S-functions are defined as [17] S(a, c; x) ≡ exp(−x/2)Φ(−a + c + 1/2, 2c where Therefore, if E → 0 and consequently, iβ → 0, it is more appropriate to explicitly use Eq. ( 21).As a result, one gets the following rapidly converging expression: , which was actually used to compose Fig. 3.
As the energy is extremely low, linear coordinate dependence can be seen in a wide region, which is in full agreement with the well-known rule where a 0 is the scattering length.From Eqs ( 18), (19), and ( 22) it follows that and this is indeed the case, while a 0 = −4312.06224.Practically the same scattering length is obtained from the real k-dependence of the phase shift (see the lower left graph in Fig. 2), using Eqs ( 13), (19), and ( 20)- (22).

An algorithm for the Riemann-Siegel function
Usually, to calculate δ (k) there is no need to make any use of Eqs ( 14)- (15).In some special cases, however, the convergence of the series (13) may be slow.Then it may be appropriate to combine Eqs ( 25) and ( 28) from [18].As a result, again after some analytical work, one comes to the following expression: where and we assumed that N = a − 1/2.
Thus, returning to the main subject of this section, we have found another option for calculating the phase shift, based on Eqs ( 14), ( 19), ( 23), (27), and ( 31)-(34): This little excursion once again demonstrates wonderful analytic properties of the Morse potential, providing an instructive example of how purely physical considerations may lead to a pure mathematical result.

SCATTERING PART OF THE KERNEL
As the necessary spectral data are now available, we can start the solution of the inverse problem.The first step is to perform the Fourier transform to fix the scattering part of the kernel (3): It can be shown that A s (−x) = [ f (x) − g(x)] /π, so that the inverse Fourier transform would result in the following formulas: Calculations according to Eqs (36), (37) represent the most challenging part of the overall procedure, because the functions f (x) and g(x) must be ascertained in a wide range with an extreme accuracy.It means, in particular, that the fast Fourier transform techniques are useless for our purposes.Fortunately, one can use an asymptotic formula [22] where the coefficients a 1 , a 3 , a 5 , . . .are directly related to the potential and its derivatives (for example, Thus the asymptotic part of the integrals (37) can be calculated analytically.In addition, one can make use of Filon's quadrature formulas (see [26], p. 890), which are perfectly suitable if x is large.Filon's formulas have been used for the range k ∈ [20, 100] when x ≥ 20, while for x < 20 a more common (but very accurate) approach was used.Namely, the domain was divided into sufficiently small intervals where the 64-point Gauss-Legendre quadrature formula is appropriate for the purpose.
The described procedure ensures the correct asymptotic behaviour of both functions defined in Eq. (37).Namely, where b = 7.252681534782e-4, c = 2.315574387346e-4 are some characteristic constants.Quite remarkably, these formulas hold to better than 10 −10 accuracy for x ≥ 30.The second terms in Eq. ( 40) are not important for the following analysis, since they mutually compensate each other.On the other hand, without these additional terms one cannot correctly perform the inverse Fourier transform according to Eq. ( 38).In fact, there is no need for this inverse transform as well, but it is recommended, to be sure that the calculated functions f (x) and g(x) are reliable.Figure 4 demonstrates that this is indeed the case.
Figures 5-7 give more detailed information about the solution procedure.(Although the curves in these figures may look like solid lines, they are assemblies of dots which represent the actual results of calculations.)Looking at Figs 5-7, one can understand why the utmost accuracy is needed in performing the Fourier transform.This is most clearly seen in the lowest graph of Fig. 7 (the dashed line there shows A(0)).Indeed, if the Fourier transform is not sufficiently correct, the important information about the scattering properties of the system (shown in Figs 5,6 and in the two upper graphs of Fig. 7) may easily be lost against the background of the bound states' contribution to the kernel, which dominates at small values of the argument.Thus the whole procedure would become meaningless.
As the result of the Fourier transform according to Eq. (37), one gets two smooth functions that can be very accurately (with absolute error less than 10 −6 ) approximated by piecewise rational functions whose parameters are given in Tables 1 and 2. Using these parameters and applying Eq. (40) for larger arguments, one can reproduce all curves shown in Figs 5-7.

SOLUTION OF THE MARCHENKO EQUATION
Having completed the preliminary work, there remains the last step -the actual solution of the Marchenko equation (2).To some surprise, this is technically much easier than determining the kernel (3).For any fixed r, let us define T (x) ≡ A(r, r + x), t ≡ r + x, and s ≡ r + y.Equation (2) can then be rewritten as where F(y) ≡ A 0 (2r + x + y) T (y) , while X i and W i are the nodes and weights of an appropriate quadrature formula, respectively.Equation (41) has been solved using the Nyström method (see, e.g., [27], p. 782), which is perfectly suitable for our purposes.As the asymptotic behaviour of the kernel is known, a good idea is to split the integral into two parts and discretize them differently.Namely, for I 1 ≡ ∫ R 0 F(y)dy the 64-point Gauss-Legendre quadrature formula can be successfully used.More specifically, the family of isospectral potentials has been calculated with R = 15, and splitting the range [0, R] into two subdomains, so that there are 128 nodes within this region.
To discretize the asymptotic integral I 2 ≡ ∫ ∞ R F(y)dy, a useful trick is to change the variable [28, p. 43]: y = R+α 0 (1−u)/(1+u), where α 0 > 0 is an arbitrary constant.Correspondingly, the domain of integration reduces to [−1, 1] , where one can again use the 64-point Gauss-Legendre quadrature formula.Thus where 2 , and x i , w i are the nodes and weights of the Gauss-Legendre formula, respectively.
The computational procedure itself is simple.Namely, using the discretized version of Eq. (41) at the specified node points x n , one gets a system of linear algebraic equations for the quantities T n ≡ T (x n ).Such a system can be easily solved.In this paper, the Householder reduction algorithm [29] was used for this purpose, because it is numerically more stable than Gaussian elimination.Having found the solution, one uses Eq. (41) once again, taking x = 0 to get T (0) = A(r, r).Finally, one calculates the potential according to Eq. ( 6).
The results of solving the Marchenko equation are shown in Figs 8 and 9.The curves there were calculated with α 0 = ∆ • (1 + x 0 )/(1 − x 0 ) and ∆ = 10000 (this parameter characterizes the width of the domain).As can be seen, the discretization of the Marchenko equation with only 192 nodes related to the relevant quadrature formulas is sufficient to accurately solve the inverse problem.In addition, one can convince himself/herself that the sign of the second term in Eq. ( 3) is indeed 'minus', while s 2 0 = y 2 0 α = 6 exp(10/3), as expected.Figure 8 demonstrates the dependence of the solution A(r, r + x) on both arguments, r and x.In this case, the exact theoretical value for the norming constant (s 2 0 = y 2 0 α) was used to fix the kernel.Actually, to calculate the potential according to Eq. ( 6), one only needs the solution at x = 0, which is shown as a background contour (dash-dot line).
In Fig. 9, a family of isospectral potentials can be seen.It was obtained by varying the norming constant s 0 in a wide range (from 0 to 100), while the scattering part of the kernel remains exactly the same for all Table 1.The fitting parameters for the curves shown in Fig. 5 Table 2.The fitting parameters for the curves shown in Fig. 6 these curves.Looking at Fig. 9, one can imagine that if the real norming constant is not known, a criterion of 'reasonableness' could be used to fix it.For example, assuming that the potential has only one extremum point (minimum) and only one inflection point (as the Morse potential), the correct norming constant can be determined with reasonable accuracy.

CONCLUSION
In this paper, a detailed description of the full solution procedure of the Marchenko integral equation was given.To this end, a complete set of the necessary input data was used, which was obtained by accurately solving the direct Schrödinger problem for a model system.Such a tautological combination of the direct and inverse scattering problems may seem useless for practical purposes, but one has to bear in mind that otherwise the full set of necessary data cannot be obtained at all.Unfortunately, the methods of the inverse scattering theory are still not implemented as planned.These methods have mainly been used to construct families of isospectral potentials [30,31], and they have many interesting applications in the soliton theory (see, e.g., [32,33]), but their capabilities have not been fully exploited.Hopefully, the results of this work clearly demonstrate the real power of the Marchenko method, albeit the necessary input data were obtained   computationally, not experimentally.On the other hand, the model system used for this purpose, the wellknown Morse potential, continues to surprise us by its undiscovered analytical properties.For example, in this paper a simple formula for calculating the phase shift was derived, along with the new analytic algorithm for calculating the Riemann-Siegel function, which is a pure mathematical result.How could the results of this work be useful for further studies?An interesting option is to combine Eq. ( 2) with the Marchenko differential equation (see, e.g., [12], p. 78), assuming that the real potential can be expressed as a sum V (r) = V 0 (r) + ∆V , where V 0 (r) is a known model potential whose spectral characteristics are known as well.Thus according to Eq. ( 6), V 0 (r) = −2 [A 0 (r, r)] ′ , where A 0 (r,t) is the corresponding solution to Eq. ( 2).Then, assuming that V (r) = −2 [A(r, r)] ′ and A(r,t) = A 0 (r,t) + ∆A, one gets the equation If V 0 (r) is sufficiently close to the real potential, the right side of Eq. (43) would contain two small quantities, ∆V and ∆A, which, in principle, can be determined self-consistently along with solving this equation.

Fig. 1 .
Fig. 1.An illustration to the physical background of the proposed approach.

Fig. 2 .
Fig.2.Phase shift for the potential Eq. (1) depending on the wave number (parameters are specified in Section 2).Lower graphs demonstrate the asymptotic behaviour of δ (k).

Fig. 3 .
Fig.3.Wave functions for the scattering states at extremely small energies.The same set of wave functions is shown in all graphs and their insets, but note that the coordinate ranges are essentially different.

Fig. 7 .
Fig. 7. Coordinate dependence of the scattering part of the kernel according to Eq. (3) (two upper graphs) and the full kernel for the examined Marchenko equation (bottom graph).The dashed line shows A(0).