# Exact solution of the Schrödinger equation with a Lennard-Jones potential

###### Abstract

The Schrödinger equation with a Lennard-Jones potential is solved by using a procedure that treats in a rigorous way the irregular singularities at the origin and at infinity. Global solutions are obtained thanks to the computation of the connection factors between Floquet and Thomé solutions. The energies of the bound states result as zeros of a function defined by a convergent series whose successive terms are calculated by means of recurrence relations. The procedure gives also the wave functions expressed either as a linear combination of two Laurent expansions, at moderate distances, or as an asymptotic expansion, near the singular points. A table of the critical intensities of the potential, for which a new bound state (of zero energy) appears, is also given.

## 1 Introduction

The interaction between two atoms is frequently represented by means of a Lennard-Jones potential,

(1) |

alone or with addition of some corrections. In this expression is the reduced mass of the system of two atoms, is the equilibrium distance (minimum of ) and is a dimensionless parameter accounting for the intensity of the interaction. Both and are empirically adjusted for each particular kind of interacting atoms. Other classical interatomic potentials, like the Morse, Rydberg or Buckingham ones, can be simulated, as shown by Lim [1], by one of the Lennard-Jones type.

Given a diatomic system and assumed a certain potential to represent the interaction, one is interested, from a theoretical point of view, mainly on the determination of its spectrum of energies, to be compared with the experimentally observed bound states. Nevertheless, in many cases, one needs to know also the corresponding wave functions in order to compute the expected values of quantities that may be obtained in the experiment. A large variety of algebraic methods are discussed in the monographs by Fernández and Castro [2] and by Fernández [3]. References to later developments can be found in recently published papers [4, 5, 6, 7]. Numerical methods have been developed, among others, by Simos and collaborators [8, 9, 10, 11]. An extensive bibliography concerning those methods can be found in Section 2 of a recent paper [12]. Except for a few familiar potentials, for which the differential equation can be solved exactly [13], those methods provide only with approximate values of the energies and wave functions. This may be sufficient in most of cases. However, due to the strong singularity at the origin of the Schrödinger equation with a Lennard-Jones potential, those approximate methods cannot represent faithfully the behaviour of the wave function in the neighbourhood of the origin. This fact, besides of being unsatisfactory from a mathematical point of view, may constitute a serious inconvenient for the computation of the expected values of certain operators.

The purpose of this paper is to call the attention of users of the Lennard-Jones potential towards a method of solution of the Schrödinger equation that is able to give the correct behaviour of the wave function in the neighbourhood of the origin and the infinity, the two singular points of the differential equation. The method is exact, free of approximations, although errors due to the computational procedure are unavoidable. But these errors can be reduced by increasing the number of digits carried along the calculations.

We present, in the next Section, fundamental sets of solutions of the Schrödinger equation that serve as a basis to express the physical solution. The requirement of a regular behaviour of this solution at the singular points establishes a condition, in terms of the connection factors, to be fulfilled by the energies of the bound states. The procedure to determine the connection factors is explained in Section 3. The energies of the bound states in a potential of intensity are shown in Figure 1. Expressions of the corresponding wave functions are given in Section 4. As increases, new bound states appear. We denote as critical those values of for which a state of zero energy exists. In Section 5, a method is suggested to find those critical intensities, which are reported in Table 5. Section 6 contains some pertinent comments. Finally, we recall, in an Appendix, a procedure to solve the nontrivial problem of finding the Floquet solutions.

## 2 Solutions of the Schrödinger equation

For a given energy and angular momentum , the Schrödinger equation for the reduced radial wave function, , of a particle of mass in the potential , given in Eq. (1), reads

(2) |

As usual, we will express the solutions of this differential equation in terms of dimensionless radial variable, , and energy parameter, , defined by

(3) |

For the radial wave function in terms of the new variable we will use

(4) |

Then, the Schrödinger equation becomes

(5) |

This differential equation presents two irregular singular points: one of rank 5 at the origin, an another of rank 1 at infinity. The physical solution must be regular at both singular points. To express this solution, we find convenient to consider three different fundamental systems of solutions.

### 2.1 Floquet solutions

Except for certain particular values of the parameters and , that we exclude from this discussion, there are two independent Floquet or multiplicative solutions expressed as Laurent power series of the form

(6) |

The indices are not uniquely defined. They admit addition of any integer (with an adequate relabeling of the coefficients). In the general case, the indices and the coefficients may be complex. The requirement that be a solution of (5) gives the recurrence relation

(7) |

The solution of this difference equation is not trivial. It can be treated as a nonlinear eigenvalue problem. In Appendix A we show an implementation of the Newton method to determine the indices and the coefficients .

### 2.2 Thomé solutions for large values of

There are two other independent solutions characterized by their behaviour for , namely

(8) |

It can be easily checked, by taking

(9) |

and coefficients given by (omitting the second subindex, )

(10) |

that the right hand side of Eq. (8) is a solution of the differential equation (5). In fact, it is a formal solution, as the series is an asymptotic one that does not converge in general. The two values of the subindex in Eq. (8) correspond to the two possible values of the right hand side of Eq. (9). In the case of negative energies, we adopt the convention

(11) |

Accordingly, is physically acceptable, as it vanishes at infinity, whereas diverges and, therefore, should be eliminated from the physical solution. In the case (not to be considered in this paper) of positive energies, both and are oscillating solutions and correspond to incoming and outgoing waves.

### 2.3 Thomé solutions near the origin

In the neighbourhood of the origin, the role analogous to that of and at infinity is played by two other solutions, and , such that, for ,

(12) |

Substitution of these expressions in Eq. (5) gives for the coefficients in the exponents

(13) |

and for the coefficients in the series (omitting the second subindex, )

(14) |

a recurrence relation that allows one to obtain the by starting with

(15) |

The two solutions correspond to the two possible values of the right hand side of the first of Eqs. (13). By convention we take

(16) |

Then, is acceptable, from the physical point of view, whereas should be discarded.

### 2.4 The physical solution

As the solutions and of the differential equation constitute a fundamental system, any solution can be written as a linear combination of them. In particular, the physical solution would be

(17) |

with constants and , to be determined, such that becomes regular at the origin and at infinity. To impose this condition we need to know the behaviour of and at the singular points. In other words, we need to calculate the connection factors defined by

(18) | |||||

(19) |

In terms of them, the behaviour of the physical solution in the neighbourhood of the singular points would be

The regularity of the physical solution at the singular points is guaranteed if and are chosen in such a way that

(20) |

which is possible if and only if

(21) |

For given values of the parameters of the potential, the left hand side of this equation is a function of whose zeros correspond to the values of the energies of the bound states. Equation (21) is, therefore, the quantization condition. Solving it requires to know the connection factors. We present in the next Section our procedure to determine them.

## 3 The connection factors

Let us design by the Wronskian of two functions and ,

(22) |

Then, from Eqs. (18) and (19), one obtains immediately

(23) | |||||

(24) |

All Wronskians in these equations are independent of . Those in the denominators can be calculated directly to obtain

(25) | |||||

(26) |

The calculation of the numerators is not so easy. In a former paper [14] we suggested a procedure that has been used to find the bound states in a spiked harmonic oscillator [15]. For convenience of the reader, we recall here the procedure, adapted to the present problem.

We consider firstly the Wronskians of each one of the Floquet solutions with the two Tomé solutions at infinity, (). Let us introduce the auxiliary functions

(27) |

Obviously,

(28) |

Both sides of this equation obey the first order differential equation

(29) |

A direct computation of the left hand side of Eq. (28), by using the definitions (27) and the expansions (6) and (8), gives the doubly infinite series

(30) |

whose coefficients

(31) |

are solution of the first order difference equation

(32) |

An expansion of the right hand side of Eq. (28), analogous to that in (30), can be obtained by making use of the so-called Heaviside’s exponential series [16]

(33) |

By taking and choosing , one gets an expansion,

(34) |

in series of the same powers of as in (30) with coefficients obeying the same first order difference equation,

(35) |

Both solutions and of the difference equation must be related by a multiplicative constant that, in view of Eq. (28), shold be . Therefore,

(36) |

an expression that, together with Eq. (25), would allow one to calculate the connection factors given by Eq. (23). Nevertheless, the validity of Eq. (36) is subordinate to the fulfilment of the condition , necessary for the validity of Eq. (34). Such condition is satisfied in the case , as, for , . There is no difficulty in computing by substituting, in the second of Eqs. (23),

(37) |

In the case , instead, the above mentioned condition is not satisfied and Eq. (36) is not valid for . In fact, the positive real semiaxis is a Stokes ray for , that should be taken as the average

(38) |

of its values in the sectors separated by the ray. Equivalently, one may define

(39) |

an average of the Wronskians for slightly above and below the positive real semiaxis. The result is

(40) |

This equation provides with the needed value of the numerator in the first of Eqs. (23).

The procedure to calculate the Wronskians, , () of each one of the Floquet solutions with the two Thomé solutions at the origin is analogous to that just described, with the unavoidable differences due to the fact that the singularity at the origin is of rank five, whereas it was of rank one at infinity. The auxiliary functions are now

(41) |

Then,

(42) |

For the left hand side we have the doubly infinite series

(43) |

with coefficients

(44) |

which obey the fifth order difference equation

(45) |

Five independent solutions of this difference equation are constituted by the coefficients of the five Heaviside’s exponential series

(46) |

with

(47) |

Then, analogously to Eqs. (37) and (40), one has

(48) | |||||

(49) |

Now it is immediate to calculate the connection factors and by means of Eq. (24).

## 4 Bound states

By using the above described procedure, we have determined the values of which are solution of Eq. (21) for different intensities of the potential in the range and for five values of the angular momentum, . The results are shown graphically in Figure 1.

Besides the energies of the bound states, our procedure gives also their wave functions. For the values of satisfying Eq. (21), and can be determined, save for a common arbitrary multiplicative constant, by using any one of Eqs. (20). To fix the arbitrary constant, we may impose, for instance, that

(50) |

Then

(51) |

and, in view of Eqs. (17) and (6), the wave function of the bound state becomes

(52) |

being a normalization constant such that

(53) |

For large values of , the series in Eq. (52) converge slowly and are not convenient for the computation of . In this case, it is preferable to use the asymptotic expansion

(54) |

stemming from

(55) |

bearing in mind Eqs. (20) and (50) and the expansion in Eq. (8). For the same reason, one should use the asymptotic expansion

(56) |

in the neighbourhood of the origin.

We have obtained, by way of illustration, the parameters of the four existing bound states in a potential of intensity . Tables 1 to 4 show the values of the energy, the indices of the Floquet solutions, the connection factors, and the coefficients to be substituted in Eq. (52), for each one of those bound states.
For the determination of the indices and the coefficients of the Floquet solutions, we used the Newton iteration method, to be recalled in the Appendix. We benefited from the subroutines `bandec`

and `banbks`

[17, pp. 45–46] to obtain the initial values, and from `ludcmp`

and `lubksb`

[17, pp. 38–39] in the iteration process. Double precision Fortran was used in the computation. The iteration was stopped when the correction in the absolute value of became less than . Usually, two or three iterations were enough. Simultaneously, the coefficients , with , were obtained. (Due to the fact that Eq. (7) relates coefficients with subindexes of the same parity, the ambiguity in the definition of , mentioned in Subsection 2.1, allows one to cancel all coefficients with odd .)
According to the condition (63), to be justified in the Appendix, the indices of the Floquet solutions either are real or, being complex, have opposite imaginary parts. In this case, thanks to the ambiguity in the definition of the , one may choose them to be complex conjugate to each other. Then, , , and are the complex conjugate of, respectively, , , and . Consequently, becomes real.

A word of caution about the computation of the wave function is in order. Our double precision calculations have revealed that Eq. (52), with the series truncated to , allows one to obtain values of with eight correct significant digits whenever roughly , whereas the asymptotic expansions in Eqs. (54) and (56) become useful for and , respectively. Therefore, double precision is not sufficient for a computation of the values of in the whole interval . Quadruple precision calculations, instead, provide with satisfactory results.

angular momentum | |
---|---|

energy | |

EE | |

EE | |

EE | |

EE | |

EE |

angular momentum | |
---|---|

energy | |

EE | |

EE | |

EE | |

EE | |

EE |

angular momentum | |
---|---|

energy | |

EE | |

EE | |

EE | |

EE | |

EE |

angular momentum | |
---|---|

energy | |

EE | |

EE | |

EE | |

EE | |

EE |

## 5 Critical values of the intensity

It may be interesting to know the values of for which a new bound state (of zero energy) appears. Our method of solution of the Schrödinger equation is also applicable in this case, but in a much simpler form. For zero energy, the singular point at infinity is a regular one and the basic Floquet solutions of the general case are replaced by Frobenius solutions whose coefficients can be obtained trivially. The procedure in this case is the same used to obtain the scattering length [18]. In fact, as it is well known, the presence of a new bound state of zero energy is revealed by a pole in the scattering length. We report, in Table 5, some critical values of the intensity for different values of the angular momentum .

7.04314 | 13.29573 | 21.48500 | 31.60949 | 43.66864 | 57.66218 |

46.61703 | 61.64985 | 78.58395 | 97.43067 | 118.19665 | 140.88604 |

121.28583 | 145.10984 | 170.82095 | 198.43005 | 227.94507 | 259.37186 |

231.08863 | 263.70031 | 298.19340 | 334.57660 | 372.85712 | 413.04084 |

376.02780 | 417.42555 | 460.70191 | 505.86378 | 552.91734 | 601.86799 |

## 6 Final comments

We have shown the applicability of our method for obtaining global solutions of the Schrödinger equation in the case of bound states in a (12,6) Lennard-Jones potential. The method can be similarly applied to any other Lennard-Jones-type potential, whatever exponents in the attractive and repulsive terms. The physical solution results as a determined linear combination of the two Floquet solutions and its asymptotic expansion at the singular points is proportional to the respective regular Thomé solutions.

Given a value of the intensity of the potential, a study of the indices of the Floquet solutions reveals that they are real for small energy. They may be taken in the interval , with . As the energy increases, increases and decreases, both approaching the value 1/2 for a certain energy. As , only one multiplicative solution exists: any other independent solution of the Schrödinger equation contains logarithmic terms. Increasing the energy makes both and to become complex, with fixed common real part equal to and opposite imaginary parts increasing with the energy. The physical wave function, however, may be taken real by adjusting the arbitrary global phase.

Special mention deserve the critical values of the intensity discussed in Section 5. Our Table 5 allows one to know immediately the number of states, of each angular momentum, bounded by a potential of given intensity.

## Appendix

We have mentioned in Subsection 2.1 that the computation of the indices and coefficients of the Floquet solutions can be treated as a nonlinear eigenvalue problem, whose solution we are going to consider in this Appendix. Along it we will omit, for brevity, the subindex in and . The condition in Eq. (6) implies that

(57) |

which allows one to truncate the infinite set of equations (7) and to restrict the label to the interval , both and being positive integers large enough to guarantee that the solution of the truncated problem does not deviate significantly from that of the original infinite one. Algorithms to solve finite-order problems have been discussed by Ruhe [19]. Here we recall the Newton iteration method suggested by Naundorf [20]. The procedure consists in moving from an approximate solution, , to another one, , by solving the system of equations

(58) | |||

(59) |

that results, by linearization [20], from (7) and from the truncated normalization condition

Obviously, the values of with or entering in some of Eqs. (58) should be taken equal to zero, in accordance with the truncation done. The iteration process is stopped when the difference between consecutive solutions, and is satisfactory. The resulting values of and may serve as initial values for a new iteration process, with larger values of and , to check the stability of the solution.

Of course, the Newton method just described needs initial values not far from the true solution. The two different values of can be obtained from the two eigenvalues

(60) |

of the circuit matrix [21] for the singular point at . The entries of that matrix can be computed by numerically integrating Eq. (5) on the unit circle, from to , for two independent sets of initial values. If we consider two solutions, and , obeying, for instance, the conditions

then

and

(61) |

The two signs in front of the square root produce two different values for , unless the parameters and in Eq. (5) be such that , in which case only one multiplicative solution appears, any other independent solution containing logarithmic terms. The ambiguity in the real part of due to the multivaluedness of the logarithm in the right hand side of (61) reflects the fact already mentioned that the indices are not uniquely defined. Notice that

(62) |

and, therefore,

(63) |

This may serve as a test for the integration of Eq. (5) on the unit circle.

Although Eq. (61) is exact, the are obtained by numerical integration of a differential equation and are not sufficiently precise. The resulting values of may only be considered as starting values, , for the Newton iteration process. As starting coefficients one may use the solutions of the homogeneous system

(64) |

with the already mentioned truncated normalization condition

(65) |

## Acknowledgments

Financial support from Departamento de Ciencia, Tecnología y Universidad del Gobierno de Aragón (Project E24/1) and Ministerio de Ciencia e Innovación (Project MTM2009-11154) is gratefully acknowledged.

## References

- [1] T.C. Lim, Connection among classical interatomic potential functions. J. Math. Chem. 36, 261–269 (2004).
- [2] F.M. Fernández, E. A. Castro, Algebraic methods in Quantum Chemistry and Physics (CRC Press, Boca Ratón, 1996).
- [3] F.M. Fernández, Introduction to Perturbation Theory in Quantum Mechanics (CRC Press, Boca Ratón, 2001).
- [4] K.J. Oyewumi, K.D. Sen, Exact solutions of the Schrödinger equation for the pseudoharmonic potential: an application to some diatomic molecules. J. Math. Chem. 50, 1039–1050 (2012).
- [5] H. Akcay, R. Sever, Analytical solutions of Schrödinger equation for the diatomic molecular potentials with any angular momentum. J. Math. Chem. 50, 1973–1987 (2012).
- [6] M. Hamzavi, S.M. Ikhdair, K.-E. Thylwe, Equivalence of the empirical shifted Deng-Fan oscillator potential for diatomic molecules. J. Math. Chem. 51, 227–238 (2013).
- [7] K.J. Oyewumi, O.J. Oluwadare, K.D. Sen, O.A. Babalola, Bound state solutions of the Deng-Fan molecular potential with the Pekeris type approximation using the Nikiforov–Uvarov (N–U) method. J. Math. Chem. 51, 976–991 (2013).
- [8] T.E. Simos, J. Vigo-Aguiar, A symmetric high order method with minimal phase-lag for the numerical solution of the Schrödinger equation. Int. J. Modern Phys. C 12, 1035–1042 (2001).
- [9] T.E. Simos, J. Vigo-Aguiar, An exponentially-fitted high order method for long-term integration of periodic initial-value problems. Comput. Phys. Commun. 140, 358–365 (2001).
- [10] T.E. Simos, J. Vigo-Aguiar, A dissipative exponentially-fitted method for the numerical solution of the Schrödinger equation and related problems. Comput. Phys. Commun. 152, 274–294 (2003).
- [11] J. Vigo-Aguiar, H. Ramos, Variable stepsize implementation of multistep methods for . J. Comput. Appl. Math. 192, 114–131 (2006).
- [12] T.E. Simos, New high order multiderivative explicit four-step methods with vanished phase-lag and its derivatives for the approximate solution of the Schrödinger equation. Part I: Construction and theoretical analysis. J. Math. Chem. 51, 194–226 (2013).
- [13] S. Flugge, Practical Quantum Mechanics (Springer, New York, 1974).
- [14] F.J. Gómez, J. Sesma, Connection factors in the Schrödinger equation with a polynomial potential. J. Comput. Appl. Math. 207, 291–300 (2007).
- [15] F.J. Gómez, J. Sesma, Spiked oscillators: exact solution. J. Phys. A: Math. Theor. 43, 385302 (2010).
- [16] G.H. Hardy, Divergent series (Clarendon Press, Oxford, 1949).
- [17] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran 77 (Cambridge University Press, Cambridge, 1992).
- [18] F.J. Gómez, J. Sesma, Scattering length for Lennard-Jones potentials. Eur. Phys. J. D 66, 6 (2012).
- [19] A. Ruhe, Algorithms for the nonlinear eigenvalue problem, SIAM J. Numer. Anal. 10, 674–689 (1973).
- [20] F. Naundorf, Ein Verfahren zur Berechnung der charakteristischen Exponenten von linearen Differentialgleichungen zweiter Ordnung mit zwei singulären Stelle, ZAMM 57, 47–49 (1977).
- [21] W. Wasow, Asymptotic expansions for Ordinary Differential Equations (Dover, Mineola, N. Y., 2002).