An approach to calculate Franck-Condon factors involving anharmonic potentials using harmonic oscillator bases

In this contribution a discrete variable representation (DVR) approach to calculate Franck-Condon factors (FCFs) is presented. The method is illustrated using a harmonic oscillator basis. The advantage of this approach is that it is possible to calculate FCFs associated with general potentials based on the knowledge of the harmonic oscillator theory. First the coordinate and momentum representations are introduced. Based on this concept, the DVR is proposed, which allows the Franck-Condon factors to be calculated in a simple form. This method takes advantage of the momentum as generator of translations to include the displacement of the potentials. The calculation of the FCFs is exemplified for transitions involving the ground state X1Σg+ and the excited state a 1Π g of the diatomic molecule 14N2 both modelled with Morse potentials.


Introduction
Nowadays spectroscopy is a fundamental tool in physics and chemistry [1] and quantum mechanics the theory to understand the process involved in absorption and emission of light.This phenomenon may take place between ro-vibrational states of the same electronic state but also between different electronic states.In the latter case the dipole transition intensities between the vibronic states ψ v and v y¢ ¢ are proportional, as a first approximation, to the square of the overlap integral which assumes slow varying electronic transition moment and separability of electronic, rotational and vibrational parts of the wave function in accordance with the Born-Oppenheimer and rigid rotor approximations [2][3][4][5][6].Hence the greater the overlap of the vibrational wave functions the larger the transition moment.It is worth stressing the fact that vibronic transitions yield information that cannot be obtained from microwave, infrared or Raman vibrational spectroscopy since the selection rules for a given electronic state are not applied.The integrals (1) are known as the Franck-Condon factors (FCFs) [7][8][9] associated with the Franck-Condon Principle (FCP), a fundamental concept to explain vibronic spectra.The physical explanation of the FCP is included in textbooks [10][11][12][13][14], and widely discussed in several works with chemical interest [15][16][17][18][19].In this contributions we present a simple approach to calculate the overlaps (1) associated with anharmonic potentials using a basis of harmonic oscillator functions.Our method uses a matrix representation in a harmonic oscillator basis as in other cases [20].We believe the material of this contribution is suitable for both advance-level and graduate students interested in applications of the harmonic oscillator.
Given the importance of the overlaps (1), considerable attention has been paid to their calculation.Besides the fact that a harmonic potential is the simplest approximation to be considered in many problems, the harmonic oscillator describes the molecular vibrational spectra through the concept of normal modes [6,21].In this venue the simplest situation corresponds to the case of transitions between two electronic states of a diatomic molecule.Closed formulas for FCFs involving harmonic oscillators functions have been obtained for one oscillator [22][23][24][25][26][27].
Although the harmonic oscillator approximation is good for the low lying region of the spectrum, as the energy increases anharmonic effects start manifesting and consequently alternative potentials have to be considered [28].One of the most useful model corresponds to the Morse potential because it takes into account the salient features of a realistic interatomic potential, but also because it admits exact solutions of the Schrödinger equation [29,30].Taking advantage of this fact, an analytic formula for FCFs involving the Morse potential has been obtained in terms of binomial series containing polygamma function [31].Because of the sums involved, also approximate approaches have been proposed, some of them can be found in [32][33][34][35][36][37][38][39][40][41].The Morse and other anharmonic potentials have been used in molecular physics but also in condensed matter [42,43].Special attention should be paid to the work of R J Le Roy [44], where a program to solve the radial Schrödinger equation as well as 1D potentials is provided.This program include the calculation of FCFs.
As a complement to previous pedagogical contributions involving harmonic potentials [17,19,45,46], in this work we present an approach to obtain the FCFs for the Morse wave functions using a basis of harmonic oscillator functions.Although the harmonic oscillator basis cannot be apply for any system it has been proved to be an invaluable tool in a great variety of applications [47].To accomplish this goal the following three ingredients are needed: (a) the expansion of an eigenstate as a linear combination of a complete orthonormal basis set, (b) the concept of a discrete variable representation and (c) the use of the momentum as a displacement operator.These concepts are the key to establish the calculation of the FCFs, not only for Morse potentials, but also to any potential where the harmonic oscillator basis is suitable to be used to describe the eigenfunctions.The concept of a complete basis is included in any standard course on quantum mechanics, but the DVR approach has been historically presented in the framework of orthogonal polynomials defining a discrete mesh using quadratures [48][49][50][51][52][53].This is a quite complicate machinery beyond the standard undergraduate and even graduate courses.However we show how to introduce the DVR concept just considering the matrix elements of the coordinate and momentum for the harmonic oscillator in a familiar way to students [54].The DVR concept turns out to be a quite powerful approach which allows the Hamiltonian to be expressed in terms of diagonal matrices.Here the student will learn the advantages of changing the representation to deal with a given problem.The last concept regarding the momentum as a displacement operator will be introduced in terms of the concept of Taylor series expansion, and it will be used to introduce the displacement of the potentials for the calculation of the FCFs.
In the framework of the DVR approach it is relatively simple to obtain the solutions for different potentials.This fact together with appropriate choice of units, allows the calculation of FCFs in terms of a simple internal product.In order to illustrate our method, we consider the calculation of FCFs between the ground X g 1 S + and the excited state a 1 Π g of 14 N 2 modelled with Morse potentials, although more general potentials and bases may be considered [54,55].
This article is organized as follows.In section 2 the harmonic oscillator is revisited in order to establish the definition of the DVR basis.A general approach to obtain the solutions of the Schrödinger equation is also included.Section 3 is devoted to obtain the solutions associated with the Morse potential.In section 3 the method to obtain the FCFs is described.Finally, the summary and conclusions are presented.

Harmonic oscillator
In this section we introduce the concept of discrete variable representations (DVR) for the coordinates and momenta in the context of the harmonic oscillator system.We start presenting the basic concepts needed to establish our approach.The Hamiltonian for the harmonic oscillator for a particle of reduced mass μ and frequency ω is given by and eigenvalues [10][11][12] ⎛ ⎝ ⎞ ⎠ This treatment is said to be in configuration space.An alternative approach is provided by the factorization method which consists in expressing the coordinate and momentum in terms of the operators a † and a in the following form [56] ( The substitution of (5) and (6) into (2) gives rise to the Hamiltonian This approach is said to be in Fock space or second quantization.The operators a † (a) has the following action over the eigenstates (8) Using these expressions it is possible to obtain the matrix elements of the coordinate and momentum: which can be obtained without using the expressions (9), but through the recurrence relations for the Hermite polynomials.In practice it is convenient to work with dimensionless quantities and at the end introduce the units of interest.Hence we define while for a general Hamiltonian The harmonic oscillator states |n〉 constitute an infinite complete space . In practice, however, it is not possible to work with infinite spaces, and consequently we may consider the finite Ndimensional subspace Constrained to the subspace (14), we can carry out the digonalization of the matrices defined by (10) and (11) to obtain the eigenvectors where we assume dimensionless units in accordance with (12).The matrices T and W with matrix elements , respectively, are unitary and consequently T −1 = T † , with similar property for W. Notice that the eigenvectors |q i 〉 and |p i 〉 are characterized by their eigenvalues and given in terms of the eigenstates |n〉, which when projected over the position |x〉, provide the eigenstates in terms of the coordinate representation Since the states |q i 〉 are eigenstates of the coordinate operator and |p i 〉 of the momentum operator, the eigenvectors |q i 〉 and |p i 〉 are referred to be in the position and momentum representations, respectively.On the other hand, since the states (8) are eigenvectors of the Hamiltonian, they correspond to the energy representation.Coordinate and momentum representations have been used to calculate the FCFs, using the Fourier Grid Hamiltonian Method, closely related with the DVR approach [57,58] The sets of eigenvectors {|q i 〉; |p i 〉} constitute the discretization of the coordinate and momentum, respectively.In figure 1 the states ( )  q x i given by (18) has been plotted when taken a subspace (14) of dimension N = 10.We can see that they are highly localized around their eigenvalues which coincide with the localized states defined using orthogonal polynomials defining quadratures [52].
The discretized eigenvectors have the remarkable property that an arbitrary function ( ) with similar results for the matrix elements of any function ( ˆ)  F p in the momentum representation In both cases the diagonal matrices are obtained by evaluating the functions ( )  G x and ( )  F p at their corresponding eigenvalues { }   x p ; i i associated with the eigenvectors (15) and (16), respectively.This feature has an enormous importance from the practical point of view, since we can compute in simple form the matrix elements of any operator in the energy representation using the coefficients involved in the eigenvectors (15) and (16).
Let us now consider the solutions of the Schrödinger equation associated with the general potential V(x).First we express the Hamiltonian in units of ÿω and length in λ 0 units We are interested in the calculation of the matrix elements of the potential ( )  V x in the energy representation: To achieve this goal we just introduce the change of basis from the energy to the coordinate representation to obtain, in matrix form: where the diagonal matrix Λ (V ) corresponds to (20): Following a similar approach for the kinetic energy, the Hamiltonian in the energy representation using the transformation coefficients T and W, takes the form defining the coordinate representation.The maximum of the state ( )  q x i is located at the i-th eigenvalue of the coordinate diagonalization.
where Λ (p) is also a diagonal matrix with elements in the momentum representation.The diagonal matrices (26) and (28) are constructed by simple substitution of the variables by their corresponding eigenvalues.
In this analysis we preferred length units λ 0 and energy units ÿω.However in dealing with general potentials the frequency ω of the harmonic oscillator may have the role of a variational parameter, as we next discuss.

Anharmonic potentials
Here we face the problem of obtaining the solutions of the Schrödinger equation associated with anharmonic potentials through the matrix representation (27).As an example we shall consider Morse-like potentials, which are defined by where D is the depth of the potential, x 0 the equilibrium position and β is related to the range of the potential.The importance of this potential is that it presents the salient physical features for a diatomic molecules: anharmonicity and continuum.An additional advantage is that the solutions of the Schrödinger equation are well known for both the discrete and continuum part of the spectrum [30].In particular we are interested in the bound states.
The Hamiltonian associated with the Morse potential is given by x N e y L y , 3 1 where ( ) L y v s are the associated Laguerre functions, the argument y is related to the physical displacement coordinate x by y = (2j + 1)e − β x , while N j v is the normalization constant The eigenvalues are We now proceed to obtain the Morse functions (31) as an expansion of the harmonic oscillator basis using the DVR approach presented in the previous section.The goal is to pave the way to tackle the problem of the calculation of the Franck-Condon factors involving two Morse potentials.We shall take advantage that the exact solutions are known, a fact that provides a convergence test just by direct comparison of the wave functions.We start considering the Morse Hamiltonian in energy units ÿω e and length units Before applying the DVR approach we notice that the transformation coefficients T and W previously defined are given in units of harmonic oscillator length λ 0 .This is important because as we have already mentioned, the frequency ω is in practice a variational parameter.Hence it is necessary to transform p to  p and x to  x.To achieve this goal it is convenient to introduce the parameter s = ω/ω e , in terms of which we have

= =
and consequently We are now ready to apply the DVR approach using (27).Since we plan to calculate the FCFs for the molecule N 2 involving its ground state X g 1 S + and the excited state a 1 Π g , we shall consider their fits to Morse potentials.Their corresponding Morse spectroscopic parameters are displayed in table 1.Notice that j stands for the number of bound states.In figure 2 both potentials are depicted together with their vibrational levels.
The Schrödinger equation is solved by diagonalizing the matrix representation (27).The total number of bound states for the states X g 1 S + and a 1 Π g are 82 and 60, respectively, but we decided to pay attention up to v = 20 since they represent the region where the FCFs will be calculated.As we know the parameter s, corresponding to the ratio ω/ω e , has to be determined given a space dimension N. In figure 3 we plot the root mean square deviation rms of the energies for both potentials as a function of the parameter s, taking into account all the energies with a basis dimension N = 125.We can see that the ωʼs are close to the frequencies of the Morse potentials, but not equal.The optimum parameters are included in table 1, while the corresponding potentials together with their eigenvalues are displayed in table 2. As noticed the exact energies are practically reproduced.Of course, to obtain similar results in the higher energy region a larger basis dimension is needed.Another aspect to take into account is that convergence in energies does not necessarily implies convergence in wave functions.To be sure about the quality of the wave functions we can make a direct comparison with the exact wave functions when possible, or by successive comparison of the eigenstates for different basis dimensions [54].
For both potentials a basis dimension of N = 125 was taken.We wonder about the s-dependence of the basis dimension.Indeed, as the basis dimension increases, the s-dependence is lost and we are free to chose s = 1.This behavior is also manifested when we restrict ourselves to the subspace of 20 bound states.
An additional consideration to take into account is that convergence depends on the location x 0 of the minimum of the potential.In figure 4 the error Δ defined as the difference of the exact and calculated wave functions are displayed for three values of x 0 .As the distance to the origin increases the error becomes larger.This is explained by that fact that the harmonic oscillator basis is located at the origin, and consequently as the equilibrium distance x 0 increases more functions are needed to reach the same accuracy.This problem will be avoided in the next section by introducing the displacement operator.Once we know how to obtain the wave functions we are able to proceed to calculate the FCFs.

Franck-Condon factors
In the previous section the solutions associated with two Morse potential were obtained using the DVR approach.Both were independently obtained.Here we aim at the calculation of the FCFs, which implies that both descriptions has to be connected in order to simplify the calculations.This connection consists in selecting common units of energy and length for both potentials, a trick widely used [17,19].Let us call MP-I and MP-II the potentials involved in the dipole transitions.For potential MP-I we have where we have added the superindex '1' to emphasize that the Hamiltonian refers to units of potential I.Here the following definitions are given ¯( ) For the second potential, we have a similar expression (40), but referred to units of the MP-II.We are interested, however, in taking the same units taken for the MP-I.The result is  ˆ .
Notice that at this stage the displacement has not been taken into account.Both Hamiltonians have the same representation (27) with the appropriate corresponding parameters.The diagonalization of the matrix representation for the potentials (40) and (42) provides the wave functions in approximated form: Notice that the basis dimensions may be different, but the harmonic oscillator functions are the same because of the common units.Taking advantage of the orthonormality Δ is defined as the difference between the exact and calculated wave functions.As the equilibrium distances increases, the differences become larger, as expected.
Table 2. Comparison between the exact energies obtained with (33) and the approximate results provided by the diagonalization of the Hamiltonian (27): Δ = E exact − E calc .Potential I is associated with the ground state X g 1 S + , while potential II with the excited state a 1 Π g .In both cases the basis dimension was N = 125.Here the adjustment parameters s = ω/ω e displayed in table 1 , the FCFs takes the simple form and in vectorial form

»
The sum in (47) involves the minimum dimension length, a fact that helps to simplify the calculations.The expression (48) is restricted to the case of potentials with minima localized at the same equilibrium positions.We now proceed to include the displacement of the second potential.First we consider the displacement operator  T a over a scalar function f (x) defined as If we expand in Taylor series the displaced functions around α = 0, we obtain and identifying the exponential function e − α( d/ dx) , we obtain where we have introduced the momentum operator ˆ( ) . This result suggests to introduce the operator ˆ( ) and in this way take into account the relative displacement between the potentials.Indeed, in terms of the operator ˆ( ) D a , the general expression for the FCFs takes the form ò y ay » and using the expansions (44) and (45) we get ( ) Here it is convenient to identify the integral as where the diagonal matrix ( ) with units referred to the MP-I.Inserting (57) into (54), we finally have for the FCFs:

L »
Here the matrix C j v 2 2 is a column vector involving the components of the eigenstate | v j 2 2 y ñ.In our DVR treatment, however, the full set of eigenvectors are simultaneously obtained, and consequently it is possible to obtain all the FCFs of interest in matrix form.If we denote with C j the matrix including the first 20 eigenstates which are well described with N = 125, the complete set of FCFs are obtained through

L »
Using the expression (59) we have calculated the FCFs between the ground and the excited state of N 2 .In particular we focus on the series FCF( j 1 , 0; j 2 , v 2 ), because reference data are available.In table 3 the results up to v 2 = 19 are given, while in figure 5 the same results are graphically displayed.
For comparison we present in figure 6 the differences of our results and the ones calculated using other models.In particular the FCFs calculated in [31] (with differences in blue) may be considered as the exact results, since the exact wave functions associated with both potentials with the same set of parameters are used.Two additional comparisons are displayed both from calculations using RKR potentials [62].The deviations in green are taken from [60], while the deviations displayed in red were taken from the results provided by [61].In order to appreciate the accuracy in quantitative form, in table 4 the percentages of deviation with respect our calculation are displayed.It is clear that our calculation are quite close to the ones provided by [61].As we already Table 3. Squared of the Franck-Condon factors associated with the potentials displayed in figure 2 for the 14 N 2 molecule.For the potential associated with the state X g 1 S + it was considered the vibrational state with υ 1 = 0.For this calculation the adjustment parameter was taken as s 1 = 0.4 and a basis dimension N = 125.remarked, equation (60) provides the full set of FCFs.As an example of the efficiency of our approach, in figure 7 we have displayed in 3D the distribution of all the FCFs associated with quantum numbers v 1 , v 2 = 0, 1, 2,...19.Finally, we summarize the source of errors in our description.First of all we should remark that our description in given in the framework of the Born-Oppenheimer and rigid rotor approximations, which are errors attached to any approach.In addition we consider only the FCFs, which is only an approximation to the transition intensities.The source of errors belonging to our approach lies on the basis as well as its dimension.The latter can be diminished by choosing the appropriate basis to a given potential with the correct basis dimension to assure convergence.In this work we have chosen the harmonic oscillator basis not because is the best to describe the Morse wave functions, but to illustrate the DVR approach in the framework of a well known quantum mechanical system without recurring to the theory of orthogonal polynomials.As a Supplementary Material we include a code in Mathematica where all the calculations can be followed step by step.

Summary and conclusions
In this contribution we have presented the DVR approach using the harmonic oscillator basis and applied it to the Morse potential.The approach has been introduced using the theory of harmonic oscillator.First the matrix elements of the coordinate and momentum are obtained.The diagonalization of the matrices gives rise to the coordinate and momentum representations, with the remarkable feature that the matrix representation of any function of the coordinate is diagonal in the coordinate representation, with analogous behavior for the momentum.In particular, for the coordinate representation, the functions q i (x) turn out to be strongly localized at the corresponding i-th eigenvalue.Based on the properties (20) and (21), it was shown that the matrix representation of the Hamiltonian for any potential can be expressed in a simple form in terms of diagonal matrices together with the transformation coefficients connecting the energy representation.
The DVR approach was applied to the Morse potential, which presents the main features of a real diatomic molecular potential.Because of the interest in the calculation of the FCFs, two cases were studied in detail: a potential with the minimum located at the origin and a second one displaced potential.It was shown that a fast Table 4. Percentage deviations between the FCFs calculated with our model and using other methods.While in [31] the potentials are also modeled with Morse potentials with FCFs calculated between exact Morse functions, in [60,61], the RKR method is used to model the potentials.convergence is obtained for the potential located at the origin, while for the displaced potential a larger basis was needed.To avoid this problem the displacement operator given in terms of the momentum was introduced.For the calculation of the FCFs a crucial ingredient was the introduction of the same units for both potentials, a fact that reduces the FCFs calculation to a dot product between the eigenvectors.Finally using the momentum representation it was possible to obtain the FCFs embraced in a matrix involving the eigenvectors for the displaced potential together with the transformation coefficients connecting the momentum representation.This is a simple calculation which can be carried out as an exercise in courses of either quantum mechanics, once the harmonic oscillator has been presented and the concept of FCF is given.In addition, our approach shows the advantage of counting on both momentum and coordinates representations.The concept of a representation is fundamental in quantum mechanics, and we have shown its importance using the harmonic oscillator system.
In this contribution we have applied a DVR approach based on the harmonic oscillators basis to describe Morse potentials, but the method may be applied for other potential of interest.However, the harmonic oscillator functions may not work as a suitable bases for certain potentials, like the Lennard-Jones potential [55] for instance.In such cases alternative DVR approaches can be used.In particular, algebraic DVR approaches based on Morse and Pöschl-Teller functions are available [54].The selection of the method is based on the symmetry of the potential to be described.Another aspect of our approach is that it can easly extended to incorporate the non-Condon factors, as explained in [39].

Figure 2 .
Figure 2. Potentials associated with the ground and excited states of 14 N 2 which are considered for the calculation of the FCFs.

Figure 3 .
Figure 3. Plots of rms vs s = ω/ω e involving all the bound sates with N = 125 for potentials associated to the the ground (Potential I) and excited (Potential II ) states of 14 N 2 .

Figure 4 .
Figure 4. Plots of the difference |Δ| vs the position x for different equilibrium distances x 0 , considering a basis with N = 125.Δ is defined as the difference between the exact and calculated wave functions.As the equilibrium distances increases, the differences become larger, as expected.

Figure 6 .
Figure 6.Differences between the FCFs calculated with our model and using other methods.While in[31] (blue color) the potentials are also also modeled with Morse potentials with FCFs calculated between exact Morse functions, in[60,61] (green and red colors) RKR methods are used to model the potentials.

Figure 7 .
Figure 7. Complete set of FCFs provided by equation (60) up to v 1 = v 2 = 19 quantum numbers.Two perspectives of the same results: in (a) the yellow left front of the figure coincides with the plot of figure 5, while in (b) the same figure but rotated is displayed.

Table 1 .
Spectroscopic parameters for the ground state and the excited state of 14 N 2 .The relative displacement between the potentials is denoted by x 0 .The parameter κ are calculated through the relation κ = w e /x e w e .
were used.
Graphic representation of the squared FCFs presented in table3.