# Complete solution of the Marchenko equation for a simple model system/Martsenko vorrandi taielik lahendus lihtsa mudelsusteemi jaoks.

1. 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-11] (see  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.sub.n] < 0 (n = 0, 1, 2, ..., N - 1) and the full energy dependence (from 0 to [infinity]) of the phase shift [delta](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 Ne2 molecule  and its Morse approximant, calculated according to the formula 

V (r)/D = exp [- 2[alpha] (r - [R.sub.e])/[r.sub.0]] - 2 exp [- [alpha](r - [R.sub.e])/[r.sub.0]], (1)

taking D = 42.153 K [13,15] and [R.sub.e] = 3.0895 [Angstrom] . The Morse curve in Fig. 1 was calculated with [r.sub.0] [equivalent to] [square root of ([[??].sup.2]/(2mD))] = 0.2388 [Angstrom], [alpha] = 0.4879/[r.sub.0], and [DELTA]r [equivalent to] r - [R.sub.e]. It has only one bound state:

[E.sub.0]/D = -0.5716, and practically the same value has been determined for the single discrete level of the [Ne.sub.2] molecule, both theoretically  and experimentally .

[FIGURE 1 OMITTED]

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 Schrodinger equation already given by Philip Morse himself , 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 . However, in this context we are talking about the physically correct linear combination of the two special solutions to the Schrodinger 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 Schrodinger 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 Schrodinger equation. Such technically rather complex developments are not analysed in this paper (see, e.g.,  and  for more details), because here we concentrate on solving the inverse not the direct Schrodinger 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 ) 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.

2. MARCHENKO INTEGRAL EQUATION

The basis for this paper is the integral equation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (2)

derived by Marchenko . 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 , requires an additional calculation of the Jost function.

The kernel of Eq. (2) reads

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (3)

where C = [[??].sup.2]/(2m) and [s.sub.n] are the norming constants for the Jost solution of the related Schrodinger equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (4)

so that f (i[[gamma].sub.n], r) [right arrow] exp(-[[gamma].sub.n]r) as r [right arrow] [infinity], and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (5)

If one is able to solve Eq. (2), the potential can be easily determined:

V (r) = -2C dA(r,r)/dr. (6)

Equation (2) and its kernel (3) are a bit different than originally given by Marchenko , but there is no contradiction. Indeed, the original Marchenko equation (see, e.g., [20,21]) can be rewritten as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (7)

Consequently, taking [A.sub.0](x) [equivalent to] -[B.sub.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., , 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 [r.sub.0] [equivalent to] [raiz cuadrada de ([[??].sup.2]/(2mD))] = 1, D = 1, and consequently, C = 1.

For a Morse potential, the norming constants [s.sub.n] can be easily ascertained analytically. Let us fix, for example, [alpha] = 2/3. Then the system has only one bound state:

[E.sub.n] = -D [(1 - n+1/2/a).sup.2] [right arrow] [E.sub.0] = -D [(1 - 1/2a).sup.2], (8)

where a [equivalent to] [square root of (D/C)]/[alpha]. Consequently, in our unit system, a = 1/[alpha] and [[gamma].sub.0] = 2/3 = [alpha]. 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 a and 70, 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.sub.e] = 2.5[r.sub.0].

Strictly speaking, to exclude an extra bound state near E = 0, we have to set the constraint V(0) = [infinity], i.e., the potential jumps to infinity (not to zero!) as r [right arrow] 0. The problem is that if one takes [alpha] = 2/3, then a = 3/2 and, according to Eq. (8), [E.sub.1] = 0. However, this extreme case is excluded because the coordinate here ranges from 0 to [infinity], not from -[infinity] to [infinity] 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 Schrodinger equation, to fix the correct regular solution. Again, it can be shown (see, e.g., ) that the second special solution is of next to no importance, so that the physical solution reads

[[PSI].sub.0] = [N.sub.0] exp(-y/2)y (9)

as in the 'normal' Morse case. Here N0 is the related norming constant and a new dimensionless coordinate y [equivalent to] 2aexp [-[alpha] (r - [R.sub.e])] was introduced.

The parameter [N.sub.0] as well as the norming constant for the Jost solution (4) can be easily determined:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (10)

The second term in the denominator is negligible, so that, to a high accuracy, [s.sup.0] = [y.sup.2.sub.0][alpha] = 6exp(10/3). This is the value which we are trying to recover. However, to clarify the sign problem described above, the quantity [s.sup.2.sub.0] will be treated as a variational (not definitely positive!) parameter. We will see that its expected value is indeed the correct value.

3. 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

[delta](0) - [delta]([infinity]) = N[pi] (11)

(N = 1 in our case), according to the Levinson theorem [1,2]. The analytic background of the problem has been described elsewhere  (see Eqs (23)-(29) there). There are several formally equivalent approaches to building the regular solution of the Schrodinger equation for a scattering state. For example, we may proceed from Eq. (25) of . Then, after some analytical work, we get the following special solutions:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (12)

with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (13)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (14)

and

B [equivalent to] arg [[GAMMA](2i[beta])/[GAMMA](i[beta])], (15)

(N = 1 in our case), [GAMMA] being the gamma function.

From the regularity requirement [PSI](b) = [A.sub.1] [[PSI].sub.1] (b) + [A.sub.2] [[PSI].sub.2](b) = b, one gets the ratio of the coefficients

[A.sub.1]/[A.sub.2] = - tan [[[alpha].sub.0] + arg ([S.sub.0])], [S.sub.0] [equivalent to] S (a, i[beta] ; [y.sub.0]). (16)

On the other hand, since S(a, i[beta]; y) [right arrow] 1 as r [right arrow] [infinity] , one gets a formula for the phase shift:

tan [[delta](k)] = [A.sub.1]/[A.sub.2] cos [[alpha].sub.0] + sin [[alpha].sub.0]/sin [[alpha].sub.0] - cos [[alpha].sub.0], (17)

and consequently, a general expression for the wave function [PSI] = [A.sub.1][[psi].sub.1] + [A.sub.2][[PSI].sub.2] (for simplicity we drop the arguments of the S-functions)

[PSI](k, r) = [A.sub.2][absolute value of (S)]/cos [[[alpha].sub.0] + arg ([S.sub.0])] sin [kr + arg ([S.sub.0]) - arg (S)].

The last formula is valid for any r, and also in the limit r [right arrow] [infinity]. From this one immediately concludes that [A.sub.2] = - cos [[[alpha].sub.0] + arg ([S.sub.0])] and [A.sub.1] = sin [[a.sub.0] + arg ([S.sub.0])], which means that

[PSI](k, r) = [absolute value of (S)] sin [kr + arg ([S.sub.0]) - arg (S)], (18)

and most importantly,

[delta](k) = arg ([S.sub.0]). (19)

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 2 OMITTED]

Figure 3 demonstrates the behaviour of the scattering wave functions as E [right arrow] 0. In general, the wave functions can be calculated using Eqs (13), (18), and (19). Note, however, that the 5-functions are defined as 

S(a, c; x) [equivalent to] exp(-x/2)[PHI](-a + c + 1/2,2c + 1; x), (20)

where

[PHI](a, c; x) = 1 + ax/1!c + ax a(a + 1)[x.sup.2]/2!c(c+1) + ... (21)

Therefore, if E [right arrow] 0 and consequently, i[beta] [right arrow] 0, it is more appropriate to explicitly use Eq. (21). As a result, one gets the following rapidly converging expression:

S(a, i[beta]; y) = exp(-y/2)j 1 - y + i[beta]y [3 - 0!/[(2!).sup.2] y - 1!/(3!) [y.sup.2] - 2!/(4!) [y.sup.3] - ...]},

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

[delta](k) = N[pi] - arctan ([ka.sub.0]), (22)

where [a.sub.0] is the scattering length. From Eqs (18), (19), and (22) it follows that

[PSI](k, r) [right arrow] sin [kr + [delta](k)] [right arrow] sin (kr) + tan [[delta](k)] [right arrow] k (r - [a.sub.0]),

and this is indeed the case, while [a.sub.0] = -4312.06224.

[FIGURE 3 OMITTED]

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).

3.1. An algorithm for the Riemann-Siegel function

Usually, to calculate [delta](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 . Asa result, again after some analytical work, one comes to the following expression:

tan [[[alpha].sub.0] + arg(S)] = [[alpha].sup.2.sub.N][beta]/[y.sup.2N+1] exp(y)F (N, i[beta]; y), (23)

where [[alpha].sub.N] [equivalent to] [absolute value of ((i[beta] + 1) (i[beta] + 2) ... (i[beta] + N))],

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (24)

and we assumed that N = a - 1/2.

To exploit Eqs (23), (24), one takes y = [y.sub.0], uses Eq. (14), and finally gets the phase shift according to Eq. (19). The only problem is how to ascertain the parameter B defined by Eq. (15). Fortunately, we can apply to the theory of gamma functions (see , Section 1.2):

[GAMMA](1/2 + z)[GAMMA](1/2 - z) = [pi]/cos([pi]z), (25)

[GAMMA](2z) = [2.sup.2z-1]/[square root of [pi]] [GAMMA](z)[GAMMA](1/2 + z) = [2.sup.2z-1] [square root of [pi]]/cos ([pi]z) x [GAMMA](z)/[GAMMA](1/2 - z). (26)

Taking here 2z = 1/2 + i[beta], one easily gets the desired result:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (27)

where

[theta]([beta]) [equivalent to] arg [[GAMMA](1/4 + i[beta]/2)] - [beta]/2 ln [pi] (28)

is the Riemann-Siegel theta function, whose behaviour is well known as [beta] [right arrow] 0 or [beta] [right arrow] [infinity]. Namely ,

[beta] [right arrow] 0: [theta]([beta]) = -[beta]/4 [2[gamma] + [pi] + 2 ln (8[pi])] + [[beta].sup.3]/24 [[[pi].sup.3]+28[zeta](3)] - [[beta].sup.5] [[[pi].sup.5]/96 + 31[zeta](5)/10] + ... (29)

and

[beta] [right arrow] [infinity]: [theta]([beta]) = - [beta]/2 ln (2[pi]/[beta]) - [beta]/2 - [pi]/8 + 1/48[beta] + 7/5760[[beta].sup.3] + 31/80640[[beta].sup.5] + .... (30)

Here [gamma] = 0-5772156649 ... is the Euler-Mascheroni constant and [zeta](3) = 1.2020569032 ..., [zeta](5) = 1.0369277551... are the corresponding values of the Riemann zeta function. In practice, Eq. (29) can be used if [beta] < 0.1, while Eq. (30) is recommended for [beta] > 20. For the intermediate region [beta] [member of] [0.1,20] an innovation can be introduced. Namely, combining Eq. (27) with a useful formula for calculating arg [[GAMMA](i[beta] + 1/2)] , we get

2[theta]([beta]) = [beta] [ln ([square root of (1 + 4[[beta].sup.2]/4[pi]) - 1] - arctan [tanh ([beta][pi]/2)] - I([beta])/2, (31)

where 

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (32)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (33)

The integral I([beta]) can be conveniently evaluated numerically or even analytically :

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (34)

Here [B.sub.n] denotes the nth-order Bernoulli number ([B.sub.1] = 1/6, [B.sub.2] = 1/30, etc.).

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):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (35)

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.

4. 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):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (36)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (37)

It can be shown that [A.sub.s](-x) = [f (x) - g(x)]/[pi], so that the inverse Fourier transform would result in the following formulas:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (38)

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 

[delta](k) = [a.sub.1]/k + [a.sub.3]/[k.sup.3] + [g.sub.5]/[k.sup.5] + ..., k [right arrow] [infinity], (39)

where the coefficients [a.sub.1], [a.sub.3], [a.sub.5], ... are directly related to the potential and its derivatives (for example, [a.sub.1] = -[(2C).sup.-1] [[integral].sup.[infinity].sub.0] V(r)dr). Thus the asymptotic part of the integrals (37) can be calculated analytically. In addition, one can make use of Filon's quadrature formulas (see , p. 890), which are perfectly suitable if x is large. Filon's formulas have been used for the range k [member of] [20,100] when x [greater than or equal to] 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,

f(x) [approximately equal to] [pi]/4 [y.sub.0] exp (- [alpha]x/2) + b exp(-cx), g(x) [approximately equal to] [pi]/4 [y.sub.0] exp (- [alpha]x/2) - bexp(-cx), (40)

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 [greater than or equal to] 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.

[FIGURE 4 OMITTED]

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.sup.-6]) approximated by piecewise rational functions

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

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.

5. 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) [equivalent to] A(r, r + x), t [equivalent to] r + x, and s [equivalent to] r + y. Equation (2) can then be rewritten as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (41)

where F (y) [equivalent to] [A.sub.0] (2r + x + y) T (y), while [X.sub.i] and [W.sub.i] are the nodes and weights of an appropriate quadrature formula, respectively.

[FIGURE 5 OMITTED]

Equation (41) has been solved using the Nystrom method (see, e.g., , 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.sub.1] [equivalent to] [[integral].sup.R.sub.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.sub.2] = [[integral].sup.[infinity].sub.R] F (y)dy, a useful trick is to change the variable [28, p. 43]: y = R + [[alpha].sub.0] (1 - u) / ( 1 + u), where [[alpha].sub.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

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (42)

where [x'.sub.i] [equivalent to] R + [[alpha].sub.0] (1 - [x.sub.i])/(1 + [x.sub.i]), [w'.sub.i] [equivalent to] 2[[alpha].sub.0][w.sub.i]/[([x.sub.i] + 1).sup.2], and [x.sub.i], [w.sub.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.sub.n], one gets a system of linear algebraic equations for the quantities [T.sub.n] [equivalent to] T ([x.sub.n]). Such a system can be easily solved. In this paper, the Householder reduction algorithm  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).

[FIGURE 6 OMITTED]

The results of solving the Marchenko equation are shown in Figs 8 and 9. The curves there were calculated with [[alpha].sub.0] = [DELTA] x (1 + [x.sub.0])/(1 - [x.sub.0]) and [DELTA] = 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.sup.2.sub.0] = [y.sup.2.sub.0][alpha] = 6exp(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.sup.2.sub.0] = [y.sup.2.sub.0][alpha]) 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.sub.0] in a wide range (from 0 to 100), while the scattering part of the kernel remains exactly the same for all 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.

6. 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 Schrodinger 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.

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

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., , p. 78), assuming that the real potential can be expressed as a sum V(r) = [V.sub.0] (r) + AV, where [V.sub.0] (r) is a known model potential whose spectral characteristics are known as well. Thus according to Eq. (6), [V.sub.0](r) = -2 [[A.sub.0](r, r)]f, where [A.sub.0](r, t) is the corresponding solution to Eq. (2). Then, assuming that V(r) = -2 [A(r, r)]f and A(r, t) = [A.sub.0](r, t) + [DELTA]A, one gets the equation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (43)

If [V.sub.0] (r) is sufficiently close to the real potential, the right side of Eq. (43) would contain two small quantities, [DELTA]V and [DELTA]A, which, in principle, can be determined self-consistently along with solving this equation.

doi: 10.3176/proc.2016.3.07

ACKNOWLEDGEMENTS

The author acknowledges support from the Estonian Ministry of Education and Research (target-financed theme IUT2-25) and from ERDF (project 3.2.1101.12-0027) for the research described in this paper. Special thanks to the referees V. Sovkov and V. Spirko for valuable suggestions and comments on the paper. The publication costs of this article were covered by the Estonian Academy of Sciences.

REFERENCES

[1.] Levinson, N. On the uniqueness of the potential in a Schrodinger equation for a given asymptotic phase. K. Danske Vidensk. Selsk. Mat-fys. Medd., 1949, 25, 1-29.

[2.] Levinson, N. The inverse Sturm-Liouville problem. Mat. Tidsskr. B, 1949, 13, 25-30.

[3.] Bargmann, V. Remarks on the determination of a central field of force from the elastic scattering phase shifts. Phys. Rev., 1949, 75, 301-303.

[4.] Bargmann, V. On the connection between phase shifts and scattering potential. Rev. Mod. Phys., 1949, 21, 488-493.

[5.] Gel fand, I. M. and Levitan, B. M. On the determination of a differential equation from its spectral function. Izv. Akad. Nauk SSSR. Ser. Mat., 1951, 15, 309-360 (in Russian) [Am. Math. Soc. Transl. (ser. 2), 1955, 1, 253].

[6.] Jost, R. and Kohn, W. Construction of a potential from a phase shift. Phys. Rev., 1952, 87, 977-992.

[7.] Jost, R. and Kohn, W. Equivalent potentials. Phys. Rev., 1952, 88, 382-385.

[8.] Marchenko, V. A. Certain questions of the theory of second-order differential operators. Dokl. Akad. Nauk SSSR, 1950, 72, 457-460 (in Russian).

[9.] Marchenko, V. A. On the reconstruction of the potential energy from phases of the scattered waves. Dokl. Akad. Nauk SSSR, 1955, 104, 695-698 (in Russian).

[10.] Krein, M. G. Solution of the inverse Sturm-Liouville problem. Dokl. Akad. Nauk SSSR, 1951, 76, 21-24 (in Russian).

[11.] Krein, M. G. On the theory of accelerants and S-matrices of canonical form. Dokl. Akad. Nauk SSSR, 1956, 111, 1167-1180 (in Russian).

[12.] Chadan, K. and Sabatier, P. C. Inverse Problems in Quantum Scattering Theory. 2nd edn, Springer, New York, 1989.

[13.] Hellmann, R., Bich, E., and Vogel, E. Ab initio potential energy curve for the neon atom pair and thermophysical properties of the dilute neon gas. I. Neon-neon interatomic potential and rovibrational spectra. Mol. Phys., 2008,106, 133-140.

[14.] Morse, P. M. Diatomic molecules according to the wave mechanics. II. Vibrational levels. Phys. Rev., 1929, 34, 57-64.

[15.] Wuest, A. and Merkt, F. Determination of the interaction potential of the ground electronic state of Ne2 by high- resolution vacuum ultraviolet laser spectroscopy. J. Chem. Phys., 2003, 118, 8807-8812.

[16.] Gendenshtein, L. Derivation of exact spectra of the Schrodinger equation by means of supersymmetry. Pis ma v ZETF, 1983, 38, 299-302 (in Russian) [JETP Lett., 1983, 38, 356-359].

[17.] Selg, M. Numerically complemented analytic method for solving the time-independent onedimensional Schrodinger equation. Phys. Rev. E, 2001, 64, 056701-12.

[18.] Selg, M. Observable quasi-bound states of the H2 molecule. J. Chem. Phys., 2012,136, 114113-12.

[19.] Edwards, H. M. Riemann's Zeta Function. Chapter 7. Academic Press, New York, 1974.

[20.] Agranovitz, Z. S. and Marchenko, V. A. The Inverse Problem of Scattering Theory. Gordon and Breach, New York, 1963.

[21.] Marchenko, V. A. Sturm-Liouville Operators and Applications. Birkhauser, Basel, 1986.

[22.] Selg, M. Reference potential approach to the inverse problem in quantum mechanics. Mol. Phys., 2006,104, 2671- 2684.

[23.] Bateman, H. and Erdelyi, A. Higher Transcendental Functions. Vol. 1. Mc Graw-Hill, New York, 1953.

[24.] Chebotarev, L. V. Extensions of the Bohr-Sommerfeld formula to double-well potentials. Am. J. Phys., 1998, 66, 1086-1095.

[25.] Selg, M. Visualization of rigorous sum rules for Franck-Condon factors: spectroscopic applications to Xe2. J. Mol. Spectrosc., 2003, 220, 187-200.

[26.] Abramowitz, M. and Stegun, I. A. Handbook of Mathematical Functions. 10th printing, Dover, New York, 1972.

[27.] Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. Numerical Recipes in Fortran 77, The Art of Scientific Computing, Vol. 1. 2nd edn, Cambridge University Press, Cambridge, 1997.

[28.] Delves, L. M. and Mohamed, J. L. Computational Methods for Integral Equations. Cambridge University Press, Cambridge, 1985.

[29.] Householder, A. S. Unitary triangularization of a nonsymmetric matrix. J. ACM, 1958, 5, 339-342.

[30.] Abraham, P. B. and Moses, H. E. Changes in potentials due to changes in the point spectrum: anharmonic oscillators with exact solutions. Phys. Rev. A, 1980, 22, 1333-1340.

[31.] Pursey, D. L. New families of isospectral Hamiltonians. Phys. Rev. D, 1986, 33, 1048-1055.

[32.] Eckhaus, W. and van Harten, A. The Inverse Scattering Transformation and the Theory of Solitons. North-Holland, Amsterdam, 1981.

[33.] Ablowitz, M. J. and Clarkson, P. A. Solitons, Nonlinear Evolution Equations and Inverse Scattering. Cambridge University Press, Cambridge, 1991.

Matti Selg

Institute of Physics of the University of Tartu, Ravila 14c, 50411 Tartu, Estonia; mselg@ut.ee

Received 29 June 2015, revised 1 October 2015, accepted 8 October 2015, available online 27 June 2016
```Table 1. The fitting parameters for the curves shown in Fig. 5

Parameters of the rational fit for the function f(x)

0 [less than or equal to] x   0.2 < x [less than   2.5 < x [less than
[less than or equal to] 0.2    or equal to] 2.5     or equal to] 10

a         -8.0742834752365    -8.07451927770804      0.710177479533
b          -139.70663281        21.09396832161       0.078971528944
c          -1964.01791531      -20.103272203664     -0.276506389451
d          9295.18176575       10.475159850219       0.065047583073
e          -8951.34532192      -2.595277468302     -4.695098041353e-3
f                0              0.196451244097     1.138463611195e-4
g           19.59774821        -0.316168857071      -0.711060782617
h           287.27716444        0.804977555407       0.277492349329
i          -511.76899191       -0.412291982417      -0.053313022152
j          -357.91228265        0.355627589898     6.620633459645e-3
k          -747.01469965       -0.117132459825     -3.875251671984e-4
l                0              0.024943017758     1.561945583680e-5

0 [less than or equal to] x    10 < x [less than
[less than or equal to] 0.2     or equal to] 30

a         -8.0742834752365      -1.318246783662
b          -139.70663281        0.350978767195
c          -1964.01791531       -0.025997992558
d          9295.18176575       8.813724334413e-4
e          -8951.34532192     -1.492795622311e-5
f                0             1.117794694102e-7
g           19.59774821         -0.345296699762
h           287.27716444        0.068087218353
i          -511.76899191      -6.848766446371e-3
j          -357.91228265       4.849331544015e-4
k          -747.01469965      -1.655209945751e-5
l                0             4.455108455684e-7

Table 2. The fitting parameters for the curves shown in Fig. 6

Parameters of the rational fit for the function g(x)

0 [less than or equal to] x    0.25 < x [less than
[less than or equal to] 0.25    or equal to] 2.42

a         3.072175507081882      3.070830178732
b          19665382.109603       -0.08759364482
c         -57958653.686259      -24.726764443843
d          37126643.692319       34.416513599217
e          14117338.25322       -13.367864496127
f         -11218702.142292       1.423672500917
g          6400798.980337        3.338766375967
h          2752650.386055        0.570620534536
i          3804036.247907        1.953451792127
j           -49419.924115        -0.164502127028
k          -268374.263668         0.34471652522
l          1638562.625682        0.011395750684

0 [less than or equal to] x    2.42 < x [less than
[less than or equal to] 0.25    or equal to] 9.2

a         3.072175507081882    -4.755191843506e-5
b          19665382.109603       0.149350091364
c         -57958653.686259       -0.015988794423
d          37126643.692319       -0.032070876143
e          14117338.25322       6.826423264373e-3
f         -11218702.142292     -2.754362667700e-4
g          6400798.980337        -1.199750111544
h          2752650.386055        0.712576159696
i          3804036.247907        -0.216463815475
j           -49419.924115         0.03998655984
k          -268374.263668      -3.895572554791e-3
l          1638562.625682       2.046893098608e-4

0 [less than or equal to] x    9.2 < x [less than
[less than or equal to] 0.25    or equal to] 42

a         3.072175507081882     -0.350846307261
b          19665382.109603       0.087431813456
c         -57958653.686259     -5.295698261860e-3
d          37126643.692319     1.228492194728e-4
e          14117338.25322      -7.497417661209e-7
f         -11218702.142292     -9.379725943814e-9
g          6400798.980337       -0.403835495838
h          2752650.386055        0.082504174267
i          3804036.247907      -8.976469204382e-3
j           -49419.924115      6.288995558866e-4
k          -268374.263668      -2.318297714188e-5
l          1638562.625682      5.388478464471e-7
```
Title Annotation: Printer friendly Cite/link Email Feedback MATHEMATICAL PHYSICS Selg, Matti Proceedings of the Estonian Academy of Sciences Report Sep 1, 2016 6169 Ionic liquids as solvents for making composite materials from cellulose/Tsellulooskomposiitide valmistamine ioonsete vedelike abil. Assemblage of turbulent jet flows through static particulate media/Turbulentsete jugavooluste kogumi liikumine labi osakestega taidetud staatilise... Algorithms Atmospheric physics Differential equations, Partial Functions, Inverse Immunoglobulin A Inverse functions Mathematical research Partial differential equations Scattering (Mathematics)