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 , where 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 and 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 which are, respectively, the vibrational energy , the rotational constant , and the large order centrifugal distortion constants . 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 quantum numbers 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, Simon-Parr-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.
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 , v and J are, respectively, the vibrational and rotational quantum numbers, is the pure vibrational energy, the rotational constant, and are the centrifugal distortion constants (CDC) related to a potential energy curve U(r) and v. In the conventional approach of the Raleigh-Schrödinger perturbation theory (RSPT)  Albritton et al.  derived the first analytical expressions of the centrifugal distortion constants (CDC). Because of the complexity and the tedious computation of these expressions, Hutson derived, by using the Numerov difference equation, an algorithm  that allows the determination of the constants , , , and only in terms of the vibrational wave function . However, difficulties appeared in this algorithm for some potentials (like the Lennard-Jones potential) in treating high vibrational level (near dissociation for example) . Because of the development in spectroscopic techniques, higher orders are needed . Tellinghuisen  introduced an improvement to Huston algorithm, but it was not sufficient to accede these levels and reach large order of the centrifugal distortion constants.
In this work, we present the canonical function method (superposition of difference equation) which is an advantageous technique with respect to other difference equation methods for solving the vibration-rotation eigenvalue problem . Moreover, with this canonical approach with variable step , the large order of CDC and the calculation of the vibration-rotation energy levels, even near dissociation , for any potential energy curves (either RKR , ab initio [10–12], or empirical ) and for any electronic state can be done by one single and simple routine. A program is available for these calculations free of charge with the corresponding author, with applications on the six potential energy curves: Morse potential, Lenard Jones potential, RKR potential, ab initio potential, Simon-Parr-Finlin (Kratzer) potential, and Dunhum potential. Since the canonical function method is proven to be valid for any type of potential function, the program is structured in different subroutines for the different cited potentials in order to give the possibility to use the same program for other types of potentials that can be added.
2. 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 and the energy that are, respectively, the eigenfunction and the eigenvalue of the radial Schrödinger equation where , μ 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 bywhere is the value of r at the equilibrium, andwith . Equation (1) is equivalent to the Voltera integral equation in the sense that any solution of (4) is solution of (1). By inserting by its expression many times in the integral, one can find:where and are two canonical functions [20, 21] defined bywith the well determined initial valuesThe initial value for the unnormalized wave function can be deduced from and by using (5), on one hand, and on the other hand the boundary conditions 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, bywhere is the pure vibrational energy, is the rotational constant, , are the centrifugal distortion constants, is the pure vibrational wave function, and is the rotational correction.
Therefore, the energy factor can be written aswhere From this expression of , the functions and are given by By replacing , by their expressions and in (5) we findA set of Schrödinger equations is obtained by replacing the wave function by its expression 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.
2.3. Analytic Expressions of the Rotation Harmonics
2.3.1. Pure Vibration
For one electronic state and for a given potential function, the solution of the vibrational Schrödinger equation is given by where and are the pure vibration canonical functions defined by (5) and in which we replace by (i.e., we take J=0).
2.3.2. Calculation of the Rotational Harmonics
A rotational Schrödinger equation , , , is given byIf we multiply this equation by and integrate it between zero and x, it will becomeBy replacing by its expression , , in , , , we obtainwherefor x = 0Thereforefor the derivation of order iandFor the unnormalized wave function, we chose , and by using we obtainBy replacing by its value in we obtainOn the other hand, the rotation harmonics must be vanished at the boundaries ; thus givesand the rotation harmonic is given by This expression is valid without any restriction for the form of a given potential function.
2.4. Numerical Method
2.4.1. Calculation of the Vibrational Wave Function
For one electronic state and for a given potential function, the vibrational wave function is given byThe determination of requires the calculation of , , and .
(1) Calculation of and . For a given potential, and on one interval, has the polynomial formThe canonical functions and are particular solutions of , because is expanded in polynomial ; and also can be expanded asBy representing and by the same function y(r) for a given potential U(r) and energy E, the function y(r) is given byBy using , we obtain the following recursion relation:withThe initial values and are given bywhere y(0) = 1, for function and y(0) = 0, for function
Therefore, the canonical functions and are well determined at any point r.
(2) Calculation of . From equation (5), the wave function is given byBy using the boundary conditions one can findFor the unnormalized wave function , the vibration function is determined at any point r.
2.5. Diatomic Centrifugal Distortion Constants (CDC)
The rotational equations , , , and are all of the formWe multiply this equation by and integrate between and . By making use of , , , and , we obtainThen we make use of the boundary conditions for and z (at ) on one hand, and of on the other hand, and we findand similarly, for the other boundary condition (at ),The continuity condition for s(r) implies the equality of given by ; thereforeThis equation gives for the successive value of s These equations give simple expression of in terms of whereOnce is given, the determination of is reduced to that of simple definite integrals and depending on .
3. 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 , and Kratzer potential . For the investigation of the ab initio potential , 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  with the advantage of the graphical user interface GABEDIT . 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–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  U(x) = a/r6 +b/r8 +c/r10 and to the left by U(x) = d/r15 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 , by using the empirical Morse potential for the molecule CO, with those given in [14, 15] shows, respectively, excellent agreements with the relative differences and (Table 1). Similar excellent agreement = 0.00% is obtained by comparing our calculated values of the different vibrational levels with those given by Martin et al.  of the RKR potential of the molecule I2 (Table 2). Our method is also tested in case of an ab initio potential for the molecule ZnBr  (Table 3). Recently , we calculated, by using the present method, the rovibrational parameters , , of the molecule LiSr by using an ab initio potential. The comparison of the calculated values of and with those available in literature shows the good agreement with the relative differences for and acceptable agreement for . However, the comparison of our calculated values of , by the present method, with those given in literature shows a very good agreement with the relative difference for the molecule NaSr .
3.2. Advantage of the Present Method
The first explicit analytical expressions for the calculation of the distortion constants were derived by Albritton et al. . The conventional Rayleigh-Schrödinger perturbation approach makes use of the wavefunction expansion given in that leads to the expressions for as first- and higher-order perturbation energies . 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  to derive simple expressions for in terms of the internuclear distance r, , , and to elaborate an efficient routine to solve the diatomic distortion constants problem by using the Numerov method described by Cooley . By using this method, Tellinghuisen  pointed out a “pesky problem of instability” in the computation of , 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 ) which reduces the range of integration.
(ii) The instability problem and its solution presented by Tellinghuisen  are both implied by the perturbation formulation as used by Hutson  (which indirectly implies an initial value problem for ). The present work eliminates this instability problem by giving exact analytical initial values for (see and ) 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.
(iv) The computing “effort” to calculate is roughly the same as that to obtain . For each constant two new integrals, and (see ), are required. The effort necessary to calculate is roughly equivalent to the average required to solve the vibrational eigenvalue.
(v) Calculation of the energy eigenvalue 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 = 2169.81358 cm−1 and = 13.28831 cm−1. The dissociation energy for this potential is given by = 88575.80 cm−1. We present in Table SF7 (Supplementary Material file) the vibrational energy levels , the rotational constants , and the centrifugal distortion constants , , of the of the Morse potential for the molecule CO near dissociation (from v = 41 to v = 81). Our calculated value for = 88575.55 cm−1 for v = 81 is in excellent agreement with that calculated using the formula .
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 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 and the centrifugal distortion constants , , 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 of the rovibrational energy levels by using (1) and . This comparison shows the high accuracy of the present technique.
In the present work, the Volterra integral equation is used to solve the rovibrational Schrödinger equation of a diatomic molecule. This method showed that a high order precision is obtained for large order centrifugal distortion constant as . By using the present canonical function approach one can obtain these constants with the high values of vibrational levels even near dissociation by one single and simple routine. This method provides strong evidence for our assumption that the higher order of and the higher vibrational v and rotational level J for any electronic state and any type of potential energy curves (either experimental, empirical, or theoretical) are as accurate as the low-order values.
The data used to support the findings of this study are available from the corresponding author upon request.
For the calculations, a program is available, free of charge, with the corresponding author.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
In the supplementary material there are three main files: Read me: Rovib-1 is a program for calculating the rovibrational energy eigenvalues EvJ, the rotational constant Bv, and large order centrifugal distortion constants Dv, Hv, Lv,.. of bound levels of different types of potential energy curves: empirical, RKR, and ab initio. Input Data: Table of all the data used in the program. Program, input, output files: In this file the reader can find the following: (i) the Rovib-1 program; (ii) examples of input files for the different kinds of potential energy curves; (iii) examples of output files for the different kinds of potential energy curves. (Supplementary Materials)
A. Messiah, Mechanique, vol. 2, Dunod, Paris, 1972.
M. Dagher and H. Kobeissi, “DIRIGE- a program for calculating eigenvalues and initial values of log derivative eigenfunctions for a diatomic molecule,” Computer Physics Communications, vol. 46, p. 445, 1987.View at: Google Scholar
N. El-Kork, N. Abu el kher, F. Korjieh, J. A. Chtay, and M. Korek, “Electronic structure of the polar molecules XF (X: Be, Mg, Ca) with rovibrational and dipole moment calculations,” Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy, vol. 177, pp. 170–196, 2017.View at: Publisher Site | Google Scholar
M. Rasoolzadeh and R. Islampour, “Estimation of vibrational energy levels of diatomic molecules (CN, CO and CS) using Numerov algorithm and comparison with the empirical values,” Australian Journal of Basic and Applied Sciences, vol. 5, no. 12, pp. 2041–2047, 2011.View at: Google Scholar
F. Martin, R. Bacis, S. Churassy, and J. Vergès, “Laser-induced-fluorescence Fourier transform spectrometry of the XOg+ state of I2: Extensive analysis of the BOu+ → XOg+ fluorescence spectrum of 127I2,” Journal of Molecular Spectroscopy, vol. 116, no. 1, pp. 71–100, 1986.View at: Publisher Site | Google Scholar
S. Elmoussaoui, N. El-Kork, and M. Korek, “Electronic structure with dipole moment and ionicity calculations of the low-lying electronic states of the ZnF molecule,” Canadian Journal of Chemistry, vol. 95, no. 1, pp. 22–27, 2016.View at: Google Scholar
G. Herzberg, Spectra of Diatomic Molecules, vol. 56, Von Nostrand, Toronto, 1950.View at: Publisher Site
J. Lennard-Jones, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 198, no. 1, 1949.
H. J. Werner, P. J. Knowles, G. Knizia et al., version .1, a package of ab initio programs, http://www.molpro.net/info/users.
N. El-Kork, I. zeid, H. A. Razzouk, S. Atwani, R. Abou arkoub, and M. Korek, “Electronic structure with dipole moment calculations of the high-lying electronic states of BeH, MgH and SrH molecules,” Journal of Physics Communications, vol. 2, no. 5, p. 055030, 2018.View at: Publisher Site | Google Scholar