Solution of the Rovibrational Schrödinger Equation of a Molecule Using the Volterra Integral Equation

By using the Rayleigh-Schrödinger perturbation theory the rovibrational wave function is expanded in terms of the series of functions φ0, φ1, φ2, . . . φn, where φ0 is the pure vibrational wave function and φι are the rotational harmonics. By replacing the Schrödinger differential equation by the Volterra integral equation the two canonical functions α0 and β0 are well defined for a given potential function. These functions allow the determination of (i) the values of the functions φι at any points; (ii) the eigenvalues of the eigenvalue equations of the functions φ0, φ1, φ2, . . . φn which are, respectively, the vibrational energy Ev , the rotational constant Bv, and the large order centrifugal distortion constants Dv,Hv, Lv . . . .. Based on these canonical functions and in the Born-Oppenheimer approximation these constants can be obtained with accurate estimates for the low and high excited electronic state and for any values of the vibrational and rotational quantumnumbers v and J even near dissociation. As application, the calculations have been done for the potential energy curves:Morse, Lenard Jones, Reidberg-Klein-Rees (RKR), ab initio, SimonParr-Finlin, Kratzer, and Dunhum with a variable step for the empirical potentials. A program is available for these calculations free of charge with the corresponding author.


Introduction
For the gas-phase of molecules, the rovibrational spectroscopy remains the most accurate and reliable source for understanding the molecular structure and the development of quantum mechanics.The theoretical determination of the rovibrational constants enables the prediction of line positions, which can be used in guiding experimental investigations and greatly facilitating the detection of unknown molecules.The fundamental atomic and molecular data are very important for interpreting observations and modelling the atmospheres of planetary and stellar objects.In the last two decades hundreds of extrasolar planets, cool stars, the atmospheres of exoplanets, and brown dwarfs have been detected where their other proposal uses molecular states.
Since diatomic molecules exhibit very long decoherence times, the molecular electronic states (electronic, vibrational, and/or rotational) are chosen recently to act as the quantum bits (qubits).A femtosecond laser pulse can be shaped to produce high fidelity binary shaped laser pulses to act as a quantum logic gate on the chosen qubits.By using the vibrational states of a linear molecule of N atoms as the qubit basis the number of vibrational degrees of freedom is 3N -6.A huge quantity of spectroscopic data is required to implement experiments using shaped laser pulses for molecular quantum computing.However, these data are still not available in literature.
The quantum mechanical canonical function method is a powerful technique to compute the very large rotation-vibration constants.The vibration rotation energy of an electronic state of a diatomic molecule is commonly represented by where  = J(J + 1), v and J are, respectively, the vibrational and rotational quantum numbers, E v is the pure vibrational energy, B v the rotational constant, and D v , H v , L v + . . .are the centrifugal distortion constants (CDC) related to

The Theory
2.1.Vibration-Rotation Canonical Functions.In the Born-Oppenheimer approximation, the vibration-rotation motion of a diatomic molecule is described by the wave function Ψ V and the energy  V that are, respectively, the eigenfunction and the eigenvalue of the radial Schrödinger equation [18] where  = 2/ℏ 2 ,  is the reduced mass, ℏ is the Planck constant, v and J are the vibrational and rotational quantum numbers, respectively, and r is the internuclear distance.This equation can be simply represented by where  =  −   (  ) is the value of r at the equilibrium, and with  = ( + 1).Equation ( 1) is equivalent to the Voltera integral equation [19] Ψ in the sense that any solution of ( 4) is solution of (1).By inserting Ψ V () by its expression many times in the integral, one can find: where () and () are two canonical functions [20,21] defined by with the well determined initial values The initial value Ψ  V (0) for the unnormalized wave function Ψ V () can be deduced from () and () by using (5), on one hand, and on the other hand the boundary conditions [5] We find 2.2.The Rotational Schrödinger Equations.In the Rayleigh-Schrödinger perturbation theory (RSPT), the eigenvalue and the eigenfunction of (1) are given, respectively, by where  0 =  V is the pure vibrational energy,  1 =  V is the rotational constant,  2 = − V ,  3 =  V , . . .are the centrifugal distortion constants,  0 is the pure vibrational wave function, and   is the rotational correction.Therefore, the energy factor  V () can be written as where From this expression of  V (), the functions () and () are given by By replacing  V (),  V () by their expressions ( 14) and (15) in ( 5) we find A set of Schrödinger equations is obtained by replacing the wave function Ψ V () by its expression (16) in (2) to obtain a set of Schrödinger equations where the first of these equations is the radial Schrödinger equation of pure vibration.All the other equations are called the rotational Schrödinger equations.

Analytic Expressions of the
where  V () and  V () are the pure vibration canonical functions defined by ( 5) and (7) in which we replace  V () by  V () (i.e., we take J=0).

Numerical Method
2.4.1.Calculation of the Vibrational Wave Function Φ 0 ().For one electronic state and for a given potential function, the vibrational wave function is given by The determination of  0 () requires the calculation of  V (),  V (), and   0 ().
(1) Calculation of  V () and  V ().For a given potential, and on one interval,   = [  ,  +1 ] has the polynomial form The canonical functions () and () are particular solutions of (17-0), because () is expanded in polynomial [8]; () and () also can be expanded as By representing () and () by the same function y(r) for a given potential U(r) and energy E, the function y(r) is given by By using (17-0), we obtain the following recursion relation: The initial values ( 1 ) and   (  ) are given by where y(0) = 1, y  (0) = 0 for  function and y(0) = 0, y  (0) = 1 for  function Therefore, the canonical functions () and () are well determined at any point r.

Diatomic Centrifugal Distortion Constants (CDC).
The rotational equations (17-0), (17-1), (17-2), and (17-n) are all of the form We multiply this equation by  0 and integrate between  0 and ∞.By making use of (17-0), (17-1), (17-2), and (17-n), we obtain Then we make use of the boundary conditions for  0 and z (at  ≈ ∞) on one hand, and of (8) on the other hand, and we find and similarly, for the other boundary condition (at  ≈ 0), The continuity condition for s(r) implies the equality of   ( 0 ) given by (44); therefore This equation gives for the successive value of s These equations give simple expression of  1 ,  2 , . . .,   in terms of  0 ,  1 , . . .,  −1 where Once  0 =  V is given, the determination of

Results and Discussion
3.1.Calculated Data.In order to test the accuracy and the validity of the present method for the calculation of the different rovibrational parameters we used as examples for the empirical potentials Morse potential [14,15], Lenard Jones potential [4,8,22], Dunhum Potential [23], and Kratzer potential [24].For the investigation of the ab initio potential [17], we employed the state averaged Complete Active Space Self Consistent Field (CASSCF) followed by Multireference Configuration Interaction (MRCI) method with Davidson correction (+Q), single and double excitations.In these entire calculations, we used the computational chemistry program MOLPRO [25] with the advantage of the graphical user interface GABEDIT [26].The CASSCF configuration space was used as the reference in the MRCI calculations.The potential energy U(r) is determined in terms of the internuclear distance r with a step equal to 0.02 Å.The Schrödinger equation is solved by using the cubic spline interpolation between every two consecutive points.A Rydberg-Klein-Rees (RKR) potential is obtained experimentally and defined by their turning points [27][28][29].The calculation has been done by an interpolation between every two consecutive points also by using the cubic spline interpolation.If the number of the turning points is limited we extrapolate the potential function curve to the right by [30] U(x) = a/r 6 +b/r 8 +c/r 10 and to the left by U(x) = d/r 15 where the constants a, b, c can be obtained from literature (if available) and d is calculated using the used potential energy curve.The comparison of our calculated values of the vibrational energy levels E v , by using the empirical Morse potential for the molecule CO, with those given in [14,15] shows, respectively, excellent agreements with the relative differences 0.00% < ΔE v /E v < 1.25% and 0.00% < ΔE v /E v < 0.03% (Table 1).Similar excellent agreement ΔE v /E v = 0.00% is obtained by comparing our calculated values of the different vibrational levels E v with those given by Martin et al. [16] of the RKR potential of the molecule I 2 (Table 2).Our method is also tested in case of an ab initio potential for the molecule ZnBr [17] (Table 3).Recently [31], we calculated, by using the present method, the rovibrational parameters E v , B v , D v of the molecule LiSr by using an ab initio potential.The comparison of the calculated values of E v and B v with those available in literature shows the good agreement with the relative differences 0.1% < ΔE v /E v < 7.5% for E v and acceptable agreement 10.9% < ΔB v /B v < 12.9% for B v .However, the comparison of our calculated values of B v , by the present method, with those given in literature shows a very good agreement with the relative difference 1.1% < ΔB v /B v < 2.8% for the molecule NaSr [31].

Advantage of the Present Method.
The first explicit analytical expressions for the calculation of the distortion constants were derived by Albritton et al. [2].The conventional Rayleigh-Schrödinger perturbation approach makes use of the wavefunction expansion given in (11) that leads to the expressions for e l , e 2 , e 3 , . . .as first-and higherorder perturbation energies (10).These expressions of the distortion constants were considered fairly good for low-lying vibrational levels only; for higher levels they failed mainly because of the omission of the contribution from continuum levels.
A modification of the perturbation method enabled Hutson [3] to derive simple expressions for e l , e 2 , e 3 , . . . in terms of the internuclear distance r,  o ,  1 ,  2 and to elaborate an efficient routine to solve the diatomic distortion constants problem by using the Numerov method described by Cooley [32].By using this method, Tellinghuisen [4] pointed out a "pesky problem of instability" in the computation of  1 ,  2 in the Hutson approach.Therefore, he described a modification of the Hutson method in order to improve it by eliminating "nuisance instabilities in the nonclassical regions." The advantages of the method given in the present work are as follows: (i) By using the variable step method, for an empirical potential, the number of mesh steps is largely reduced.For ab initio and RKR potentials, the step is the distance between two consecutive points.In the present method we step out to the left and to right by starting from an arbitrary origin.The calculation is stopped when the wave function tends to zero (see (42)) which reduces the range of integration.
(ii) The instability problem and its solution presented by Tellinghuisen [4] are both implied by the perturbation formulation as used by Hutson [3] (which indirectly implies an initial value problem for  1 ,  2 . ..).The present work eliminates this instability problem by giving exact analytical initial values for  1 ,  2 . . .(see (32) and (41)) at an arbitrary origin.
(iii) The present method is free from any conventional problems related to those of Hutson and Tellinghuisen; it does not imply any trial value or any "matching" problem (or initial values assumptions) on any iteration; it is equally easy for low or high levels.b The second and the c third entries are, respectively, for [14,15].
(iv) The computing "effort" to calculate e l is roughly the same as that to obtain e 2 , e 3 , . ... For each constant two new integrals, I n and R n (see (49)), are required.The effort necessary to calculate B v , D v , H v , L v , M v . . . is roughly equivalent to the average required to solve the vibrational eigenvalue.
(v) Calculation of the energy eigenvalue E v near the dissociation energy is another advantage.The spectroscopic constants used in the present work for Morse potential of the molecule CO [14,15] are  e = 2169.81358cm −1 and  e x e = 13.28831cm −1 .The dissociation energy for this potential is given by D e =  e 2 /4 e x e = 88575.80cm −1 .We present in Table SF7   b The second entry is for [16].v = 41 to v = 81).Our calculated value for E v = 88575.55cm −1 for v = 81 is in excellent agreement with that calculated using the formula D e =  e 2 /4 e x e .We are now considering the application of this method as a verification procedure, for additional future and previously characterized molecules.The calculated values of the rovibrational energy levels E vJ of the potential energy curves Morse, RKR, and ab initio for different values of v and J are given in Tables (SF1, SF2, SF3) in Supplementary Material file.The accuracy of our calculated values of the rotational constants B v and the centrifugal distortion constants D v , H v , L v for the six considered potentials by using the present software is presented in Tables (SF1, SF4, SF5-SF6) in Supplementary Material file.In these tables, we compare the calculated values 1 ,  2 , . . . is reduced to that of simple definite integrals  0 ,  1 ,  2 , . . .and  0 ,  1 ,  2 , . . .depending on  0 ,  1 ,  2 , . . . .[19].

Table 1 :
The energy eigenvalues E v , the rotational constants B v , and the centrifugal distortion constants D v , H v , L v of the pure vibrational energy levels of the Morse potential of CO molecule.

Table 1 :
Continued The first entry is for the present work. a (Supplementary Material file) the vibrational energy levels E v , the rotational constants B v , and the centrifugal distortion constants D v , H v , L v of the of the Morse potential for the molecule CO near dissociation (from

Table 2 :
The energy eigenvalues E v , the rotational constants B v , and the centrifugal distortion constants D v , H v , L v of the pure vibrational energy levels of the RKR potential of the electronic state XO + g of the molecule I2 .vE v cm −1 B v × 10 2 cm −1 D v × 10 9 cm −1 -H v × 10 16 cm −1 -L v × 10 22 cm −1The first entry is for the present work. a

Table 3 :
The eigenvalues E v , the rotational constants B v , and the centrifugal distortion constants D v , H v , L v of the pure vibrational levels of the ab initio potential of the electronic state (2) 2 Σ + of ZnBr molecule.v E v B v × 10 2 cm −1 D v × 10 8 cm −1 H v × 10 14 cm −1 L v × 10 19 cm −1