Generalized Sturmian Functions in prolate spheroidal coordinates
GGeneralized Sturmian Functions in prolate spheroidalcoordinates
D. M. Mitnik , F.A. L´opez , and L. U. Ancarani Instituto de Astronom´ıa y F´ısica del Espacio (IAFE), CONICET-UBA,C.C. 67, Suc. 28, (C1428EGA) Buenos Aires, Argentina. Universit´e de Lorraine, CNRS, LPCT, 57000 Metz, France.June 12, 2020
Abstract
With the aim of describing bound and continuum states for diatomic molecules,we develop and implement a spectral method that makes use of Generalized SturmianFunctions (GSF) in prolate spheroidal coordinates. In order to master all computa-tional issues, we apply here the method to one–electron molecular ions and compareit with benchmark data for both ground and excited states. We actually propose twodifferent computational schemes to solve the two coupled differential equations.The first one is an iterative 1d procedure in which one solves alternately the angularand the radial equations, the latter yielding the state energy. The second, nameddirect 2 d method, consists in representing the Hamiltonian matrix in a two–dimensionalGSF basis set, and its further diagonalization. Both spectral schemes are timewisecomputationally efficient since the basis elements are such that no derivatives have tobe calculated numerically. Moreover, very accurate results are obtained with minimalbasis sets. This is related on one side to the use of the natural coordinate system and, onthe other, to the intrinsic good property of all GSF basis elements that are constructedas to obey appropriate physical boundary conditions. The present implementation forbound states paves the way for the study of continuum states involved in ionization ofone or two–electron diatomic targets. a r X i v : . [ phy s i c s . a t m - c l u s ] J un INTRODUCTION
The molecular ion H +2 , as well as the isotopic forms such as HD + or D +2 , and other one–electron diatomics such as HHe +2 or HLi +3 , are the simplest molecular quantum three-bodyproblem with Coulomb interactions. H +2 , in particular, has been largely studied since theearly days of quantum mechanics , and is presented in standard molecular physics books asit allows one to understand why molecules form. On top of being important in astrophysics(it is involved in many reaction chains leading to the production of polyatomic molecules),the molecular ion H +2 also serves as benchmark to test any new molecular approach andnumerical method.In the fixed–nuclei approximation, it is well known that prolate spheroidal coordinatesmake the Schr¨odinger equation separable . Aside from the simple azimuthal angle depen-dence due to axial symmetry, the wavefunction depends on two variables, one angular andone radial (actually quasi–angular and quasi–radial). The H +2 bound structure can be foundby solving a system of two coupled ordinary differential equations, one for each of thesetwo variables. An analytical solution exists formally but involves two not so tractableexpansions and therefrom complicated energy equations (see, e.g. , and references therein).In practice, therefore, the energies are found numerically. This is why a wide variety ofmethods, including iterative methods, have been proposed and applied to solve the coupledequations. For continuum states, necessary for example to describe ionization processes fromdiatomic molecules, the energy is known and fixed. However, these non–L states are muchmore difficult to build as they oscillate up to infinity. Some recent investigations dedicatedto their description in prolate spheroidal coordinates include Ref. . Approximate single ordouble continuum wavefunctions borrowed from the atomic literature have been extendedto the two–center case and employed to study ionization processes . Other approachesconsist in extending well established atomic numerical techniques to the diatomic molecularcase, using (see, e.g. , ) or not using (see, e.g. , ) prolate spheroidal coordinates.In the last decade, a spectral method named GSF has been developed and implementedto study the structure of and scattering processes on atomic systems . The method usescomplete and orthogonal basis sets of Generalized Sturmian Functions (GSF) with appro-priate boundary conditions. Negative energy GSFs allow one to study bound states. Thehelium atom, the simplest atomic quantum three–body problem with Coulomb interactions,served as a benchmark to put the method on solid grounds, by studying in details con-vergence issues, the integrals involved and the adequate choice of optimal parameters andnumerical packages (see and references therein). While the aim of the GSF method wasnot to compete with well established structure codes, it proved to be very accurate at areduced computational cost because of intrinsic GSF properties in particular the adequate,and unique, asymptotic decay of all basis elements.After bound states, the GSF approach was rapidly implemented for continuum states forwhich the good properties of positive energy GSFs demonstrated the power of the method.Indeed, for continuum states, the correct asymptotic behavior is crucial in any scatteringcalculation as shown in applications to one and two–electron atomic systems (see, e.g. , ).The method was first presented in spherical coordinates, then extended to hypersphericalcoordinates but limited to atomic systems. An extension to molecules with a heavy cen-tral nucleus has been proposed in a one–center GSF approach and applied to ionization2rocesses . Nothing, however, has been proposed to deal with diatomic molecules.The purpose of this manuscript is to develop and implement a GSF method in prolatespheroidal coordinates, thus combing the two advantages of (i) using the natural coordinatesfor diatomic systems and (ii) the power of a spectral method together with the intrinsicallygood GSF properties. The long term aim is to be able to describe accurately single ordouble ionization of diatomic molecules treated as a two–electron system. The developmentwill follow a path similar to the one adopted for the atomic case. We will first considerbound one–electron molecules before moving to the continuum part of the spectrum. Bystudying benchmark one–electron molecular ions, such as the H +2 , we wish to validate thenew computational procedure and code, check thoroughly all convergence and precisionissues, and test the robustness with respect to the variation of the internuclear distance.We actually present here two different computational methods. In the first one, weadopt an iterative approach, solving alternately the separated Schr¨odinger equations forthe angular part and for the radial part. This iterative 1d procedure, which is repeateduntil convergence, presents the novelty of using GSF with appropriate boundary conditions.Because of such property the approach results to be computationally efficient as only smallbasis are needed to obtain very good energy levels. It is also efficient in computing timebecause the GSF basis elements already solve the Hamiltonian differential operator so thatno derivative calculation is needed at each iteration. The second method, called here the direct d method, consists in representing the Hamiltonian matrix in a two–dimensional GSFbasis set, and its further diagonalization. On top of the same advantages as the first method,the 2 d spectral approach demonstrates its full power by providing accurately many statessimultaneously, and this with very small basis.The remainder of this paper is as follows. In Sec. 2 we provide the theoretical frameworkof the proposed GSF method in prolate spheroidal coordinates. Then in Sec. 3 we apply itto the ground and first three excited states of symmetric (H +2 ) and asymmetric (HHe +2 andHLi +3 ) molecular ions. The successful comparison with benchmark data from the literatureallows us to validate the method for bound states. As indicated in the Conclusion (Sec. 5), thenext step will be to study continuum states for which positive energy GSF, with appropriateboundary conditions, will be used.Atomic units ( (cid:126) = m e = e = 1) are assumed throughout. Consider a diatomic molecular system consisting of one electron and two nuclei of arbitrarycharges Z and Z placed at a fixed distance R along a line defining the z axis; let r denotethe distance of the electron from nucleus 1 and r from nucleus 2. To simplify we neglectany nuclei finite mass effect.In prolate spheroidal coordinates, defined by ξ ≡ r + r R ; η ≡ r − r R ; φ ≡ arctan (cid:16) yx (cid:17) (1)where 1 ≤ ξ < ∞ , − ≤ η ≤ ≤ φ ≤ π , the Schr¨odinger equation for the electron3eads (cid:26) − R ( ξ − η ) (cid:20) ∂∂ξ ( ξ − ∂∂ξ + ∂∂η (1 − η ) ∂∂η ++ ξ − η ( ξ − − η ) ∂ ∂φ (cid:21) + V ( η, ξ ) (cid:27) ψ ( ξ, η, φ ) = E ψ ( ξ, η, φ ) , (2)with the electron-nuclei potential given by V ( ξ, η ) = − Z r − Z r = − R ( Z + Z ) ξ − ( Z − Z ) η ( ξ − η ) . (3)In the fixed–nuclei approximation, the internuclear distance R enters as a parameter, andthe nuclei repulsive potential energy 1 /R may be simply added. Equation (2) is separablein these coordinates, meaning that the solution is expressed as a product of three functions ψ ( ξ, η, φ ) = U ( ξ )Λ( η )Φ( φ ) . (4)The azimuthal function Φ is easily separated, and must fulfill the equationd Φd φ + m Φ = 0 , (5)whose solutions are Φ( φ ) = 1 √ π e imφ , (6)with m = 0 , ± , ± , ± , · · · . Because of the axial symmetry of the potential, m is a goodquantum number.Upon elimination of the azimuthal dependence, and defining p = − R E , a = R ( Z − Z )and a = R ( Z + Z ), the ensuing equation reads (cid:26) ∂∂ξ (cid:20) ( ξ − ∂∂ξ (cid:21) + a ξ − p ξ − m ξ − ∂∂η (cid:20) (1 − η ) ∂∂η (cid:21) − a η + p η − m − η (cid:27) U ( ξ )Λ( η ) = 0 . and is also separable. Denoting the separation constant as A , one obtains a system of twonon–trivial ordinary differential equations, a “radial” equation for U ( ξ ) and an “angular”equation for Λ( η ), (cid:20) ∂∂ξ (cid:20)(cid:0) ξ − (cid:1) ∂∂ξ (cid:21) + a ξ − p ξ − m ξ − A (cid:21) U ( ξ ) = 0 , (8a) (cid:20) ∂∂η (cid:20)(cid:0) − η (cid:1) ∂∂η (cid:21) − a η + p η − m − η − A (cid:21) Λ( η ) = 0 , (8b)which are coupled through both the scaled energy p and the coupling constant A . Stateswith different m values are not coupled, so that they can be considered independently.4n this work, we propose two different methods using a spectral approach based on GSFin prolate spheroidal coordinates. In the first – named hereafter “iterative 1 d method” – wesolve, alternately, the one–dimensional radial equation (8a), assuming a fixed scaled energy p , and solving an eigenvalue equation for the separation constant A . Then, we use thisconstant as a fixed value in the one–dimensional angular equation (8b), obtaining a newenergy p . The process is repeated until convergence is achieved. In this iterative procedure,both equations are solved by using adequate GSF basis sets and are converted into eigenvalueproblems. The main advantage of our GSF approach resides in the fact that the principalpart of these two equations (in particular, the derivatives) are already dealt with by thebasis functions; as a consequence, derivative calculations are not required at every iterationstep. In the second method, we construct a basis set composed of products of the angularand radial GSF. This two–dimensional basis is used to represent the Hamiltonian, which isdiagonalized in order to solve the whole Schr¨odinger equation (2). In this way, we obtain theeigenvalues (energies) and eigenvectors (solutions) of many states at the same time. Thismethod, here referred to as the “direct 2 d method”, while possessing the same advantagesrelated to GSF is computationally even more efficient. d method We search the solution of Eq. (8b), for a given m , as an expansion in Sturmian functionsΛ( η ) = (cid:88) j c j S j ( η ) , (9)the angular basis set being generated by solving the Sturmian equation (cid:20) ∂∂η (cid:20)(cid:0) − η (cid:1) ∂∂η (cid:21) − m − η (cid:21) S j ( η ) = − β j S j ( η ) , (10)with boundary conditions S j (1) = 1 and S j ( −
1) = ( − j for m = 0 and S j (1) = S j ( −
1) = 0for m (cid:54) = 0. The solutions are actually the well known associated Legendre polynomials , S j ( η ) = P mj ( η ), and correspond to eigenvalues β j = j ( j + 1). Figure 1 shows the first 9elements S j ( η ) for m = 0.With expansion (9) and making use of Eq. (10), the angular equation (8b) becomes (cid:88) j c j (cid:2) − β j − a η + p η (cid:3) S j ( η ) = A (cid:88) j c j S j ( η ) . (11)Multiplying from the left by S i ( η ) and integrating over the angular domain [ − , M c = A B c . (12)The matrices involve the elements[ M k ] ij = (cid:90) − S i ( η ) η k S j ( η ) dη (13)5 η -1-0.500.51 S j ( η ) Figure 1: First 9 angular Sturmian basis elements S j ( η ) for m = 0.which can be evaluated analytically using known properties of the Legendre polynomials .Those of interest here are given by (cid:2) M (cid:3) ij = 22 i + 1 ( i + m )!( i − m )! δ ij (14a) (cid:2) M (cid:3) ij = 22 i + 1 ( i + m )!( i − m )! 12 j + 1 [( j − m + 1) δ i,j +1 + ( j + m ) δ i,j − ] (14b) (cid:2) M (cid:3) ij = 22 i + 1 ( i + m )!( i − m )! 12 j + 1 (cid:20) ( j + 1 − m )( j + 2 − m )2 j + 3 δ i,j +2 + (cid:18) ( j + 1 − m )( j + 1 + m )2 j + 3 + ( j + m )( j − m )2 j − (cid:19) δ i,j + ( j − m )( j + m )2 j − δ i,j − (cid:21) , (14c)and are calculated only once, at the first iteration. The elements of the matrices M and B are given by [ M ] ij = − j ( j + 1) [ M ] ij − a [ M ] ij + p [ M ] ij (15a)[ B ] ij = [ M ] ij . (15b)Assuming a given energy value p , the angular part reduces to solving the generalizedeigenvalues problem (12), i.e. , finding the eigenvalue A (the separation constant) and theeigenvector c (the coefficients of expansion (9)). At each iteration, the matrix M is easilyrecalculated with the new energy value p . 6 .1.2 Radial equation Once the A eigenvalue is obtained from the angular equation, the scaled energy p is to befound from solving the radial equation (8a). Setting U ( ξ ) = ( ξ − | m | / f ( ξ ) removes thesingular term m / ( ξ −
1) from the differential equation. A first boundary condition islim ξ →∞ f ( ξ ) = e − pξ . (16)At the other end, when the electron is exactly in the center of the molecular system ( ξ = 1),we have for m = 0 lim ξ → f ( ξ ) = ξ − A e p ξ − a ξ . (17)because the radial equation (8a) reduces to df ( ξ ) dξ = (cid:18) p ξ − a − A ξ (cid:19) f ( ξ ) , (18)while for m (cid:54) = 0, we just require f ( ξ ) to be a regular function, so that U ( ξ ) vanishes.Similarly to the angular part, we propose an expansion U ( ξ ) = ( ξ − | m | / (cid:88) j d j S j ( ξ ) , (19)on a basis of Generalized Sturmian Functions S j ( ξ ) generated by the Sturmian equation (cid:20) ∂∂ξ (cid:20)(cid:0) ξ − (cid:1) ∂∂ξ (cid:21) + 2 ξ | m | ∂∂ξ + a ξ − p s ξ (cid:21) S j ( ξ ) = α j V s ( ξ ) S j ( ξ ) , (20)with eigenvalues α j . In Eq. (20), p s = − R E s is a parameter that can be set freely, with theenergy E s > E s < E s is arbitrary, choosing it close to that ofa desired state will make the GSF basis more efficient from a convergence point of view. V s ,known as generating potential, must be a short–range potential so that the basis elements S j ( ξ ) have an asymptotic behavior similar to (16). Moreover, since we wish S j ( ξ ) to possessalso the same ξ → U ( ξ ), the generating potentialmust obey the relation lim ξ → α j V s ( ξ ) = − A + p − p s . (21)It turns out that is convenient to choose a function nearly constant at ξ = 1, in order tostabilize the iterations. In the present work, the generating potential is chosen to be V s = 12 [1 − tanh( α ( ξ − γ ))] , (22)where the parameters α and γ determine the shape of the potential as illustrated by Figure2. For a given value of α , a larger parameter γ extends the range of the potential (for α = 1, γ approximately represents the range). On the other hand, for a fixed value of γ (solid anddotted curves), higher α parameters correspond to steeper potentials.7 ξ (a.u.) V s α = 1 γ =5α = 1 γ = 7α = 0.5 γ = 5 Figure 2: Generating potential V s , used to generate the radial GSF S j ( ξ ).At ξ → ∞ we could impose on S j ( ξ ) the boundary condition (16), but requiring simplythe basis function to vanish was found to be sufficient. On the other hand, imposing oneach element condition (17) at ξ → m = 0. We generate theSturmian functions by solving the radial equation (20) with a finite difference method . Asit is one–dimensional, we can afford using a numerical grid having a considerable number ofpoints. The first 9 basis elements for m = 0, generated with α = 1 . γ = 5, are shown inFigure 3. As j increases, these functions display an increasing number of nodes. Featuringone of the main GSF properties, all elements behave asymptotically in a unique manner,here in the same exponential manner e − p s ξ as ξ → ∞ .With expansion (19), and making use of (20), the radial equation (8a) takes the form (cid:88) j d j (cid:2) α j V s ( ξ ) + A + m + | m | (cid:3) S j ( ξ ) = (cid:88) j d j ( p − p s ) ξ S j ( ξ ) . (23)Multiplying from the left by S i and integrating over the domain [1 , ∞ [, we obtain anothergeneralized eigenvalues equation N d = λ C d (24)where the eigenvalues are λ = p − p s , and thus the corresponding energies through p = − R E/
2. Let us define the elements (cid:2) N k (cid:3) ij = (cid:90) ∞ S i ( ξ ) ξ k S j ( ξ ) dξ (25a)[ G ] ij = (cid:90) ∞ S i ( ξ ) V s ( ξ ) S j ( ξ ) dξ , (25b)8 ξ -1.5-1-0.500.511.5 S j ( ξ ) Figure 3: First 9 radial basis elements S j ( ξ ) for m = 0.that are calculated, numerically, only once. The matrices N and C have for elements[ N ] ij = ( A + m + | m | ) [ N ] ij + α j [ G ] ij (26a)[ C ] ij = [ N ] ij . (26b)Here A is a fixed parameter obtained from the previous step, when solving the angularpart. The solutions of (24) provide both the eigenvalues λ and the eigenvectors made of thecoefficients d j of the radial expansion (19).This iterative method has a significant advantage. The Hamiltonian is separated intotwo coupled equations, and both of them are one–dimensional reducing significantly thecomputational cost. Moreover, the use of expansions on GSF basis greatly simplifies thetask since each basis element already solves a substantial part of the equations, in particularthe differential operators. As a consequence, it is not necessary to solve numerically thedifferential equations at each step. Computationally, one only solves – iteratively – twogeneralized eigenvalue problems. There is, however, a drawback in this methodology: eachmolecular state requires a new basis set. This means that, from all the eigenvalues A and p resulting from the calculations, we must select only those corresponding to the eigenvectorshaving the right number of nodes. For each one of the molecular states, a different iterationprocedure is thus needed. This difficulty is avoided in the alternative method presentedhereafter. d method We propose now a method in which equation (2) is solved directly. As before, we first removethe azimuthal part and write ψ ( ξ, η, φ ) = Ψ( ξ, η )Φ( φ ) (27)9ith Ψ( ξ, η ) solution of the two-dimensional equation (cid:26) ∂∂ξ (cid:20) ( ξ − ∂∂ξ (cid:21) + a ξ − p ξ − m ξ − ∂∂η (cid:20) (1 − η ) ∂∂η (cid:21) − a η + p η − m − η (cid:27) Ψ( ξ, η ) = 0 . (28)This time we propose to expand the solution Ψ( ξ, η ) over a two-dimensional basis S ij ( ξ, η ) ψ ( ξ, η ) = (cid:88) ij a ij S ij ( ξ, η ) = ( ξ − | m | / (cid:88) ij a ij S i ( ξ ) S j ( η ) (29)where the one–dimensional Sturmian functions are obtained with the same methodologydescribed above, i.e. , from equations (10) and (20).Upon substitution of expansion (29), the two–dimensional equation (28) becomes (cid:88) ij a ij (cid:26) α i V s ( ξ ) + m + | m | + p s ξ − a η − β j (cid:27) S ij ( ξ, η )= (cid:88) ij a ij p ( ξ − η ) S ij ( ξ, η ) . (30)A matrix system is constructed by multiplying from the left by a basis element S i (cid:48) j (cid:48) ( ξ, η ) andintegrating over both ξ and η variables (note here the absence of the volume element ξ − η in spheroidal prolate coordinates). We obtain a generalized eigenvalues problem P a = λ D a , (31)in which the matrices P and D are given by[ P ] i (cid:48) j (cid:48) ,ij = α i [ G ] ii (cid:48) (cid:2) M (cid:3) j (cid:48) j + p s (cid:2) N (cid:3) i (cid:48) i [ M ] j (cid:48) j − a (cid:2) N (cid:3) i (cid:48) i [ M ] j (cid:48) j + ( m + | m | − β j ) (cid:2) N (cid:3) i (cid:48) i [ M ] j (cid:48) j (32a)[ D ] i (cid:48) j (cid:48) ,ij = (cid:2) N (cid:3) i (cid:48) i [ M ] j (cid:48) j − (cid:2) N (cid:3) i (cid:48) i [ M ] j (cid:48) j . (32b)We solve this eigenvalue problem, obtaining a solution matrix a ; each column consists ofthe coefficients vector (cid:126)a n , which expands that solution corresponding to the molecular statehaving eigenenergy λ n = p n . To be more specific, if the basis size is N , we have a = a a a . . . a N a a a . . . a N . . . . . . . . . . . . . . .a a a . . . a N a a a . . . a N . . . . . . . . . . . . . . .a ij a ij a ij . . . a Nij . . . . . . . . . . . . . . . λ = p p . . .p N (33)This direct methodology avoids iterations. Moreover, it allows us to obtain the solutionsfor many molecular states simultaneously. Since the matrices are two–dimensional, at firstsight the method seems computationally costly. However, all integrations leading to the ma-trix elements of P and D are separable and reduce to products of one–dimensional integralsas given by (32a) and (32b). 10 RESULTS
We present now the results of our calculations and make a comparison with the data providedin the literature. We start by applying the GSF iterative 1 d method for both the groundand some m = 0 excited states of the hydrogen molecular ion H +2 for which Z = Z = 1and thus a = 0. Next we consider asymmetric (heteronuclear) molecular ions with Z (cid:54) = Z .Finally, for H +2 , we will show how the GSF direct 2 d method yields the ground and severalexcited states in a single run. d method for the ground state of H +2 The best values of the energy E for the ground state 1 σ g , and the corresponding separationconstant A from the work of Scott et al. are used here as a benchmark to analyze theconvergence issues of our Sturmian method. We assume here an internuclear distance R = 2a.u., thus fixing the values of a and a . In order to solve the angular equation (8b), an initial value for the energy E (more precisely, p = 2 . ξ coordinate, only even elements S j ( η ) are included in the expansion. Asshown through Table 1, convergence towards the benchmark result A of Ref. is reachedwith just 4 elements. Number of basis elements A A in Eq. (8b) for fixed energy E = 1 . c j that allow usto construct the ground state angular solution (9) which is shown in Figure 4. The excitedstates will be discussed in the next section. 11 η -101 Λ ( η ) σ g σ u σ g σ u Figure 4: The angular Λ( η ) solutions for the four lower energy states of H +2 . Once the angular equation is solved, we turn to the radial equation. In contrast to theangular part, here the approach is completely numerical. On one hand we have to generatethe basis set S i ( ξ ) and, on the other, the matrix elements of the corresponding eigenvalueproblem (24) must be calculated numerically.The basis elements are generated by solving the Sturmian equation (20) with a numericalmethod described previously . It is based on a predictor–corrector algorithm, propagatingthe solution from the origin to some defined matching point (this is the outgoing solution),and from an effective infinite towards this point (the inward solution). The inward functionis normalized, in such a way that both solutions coincide at the matching point. If thederivatives disagree at this point, the eigenvalue is adjusted and the procedure starts again,until convergence. This algorithm, allows one to produce very accurate solutions for atomicsystems, even with a reasonably small ( ∼ . However, wenoticed that it was hard to obtain the radial solutions of Eq. (8a), even when a large numberof points was included in the numerical grid. In fact, to solve this equation appropriately,the crucial aspect resides in the fulfillment of the boundary conditions (17) at the origin.We endorsed this conclusion, trying to solve the radial equation with other methods, andusing different mathematical softwares, obtaining very different results for different numericalgrids. We even tried to solve the equation fixing the energy value to E = − . ,but the converged solutions yielded eigenvalues A too far from the correct value.We also tried to use standard diagonalization routines from lapack to solve equation(8a) directly. However, within this approach it is not simple to introduce explicitly theboundary conditions in contrast to our GSF expansion approach for which it is straightfor-ward. Thus, our method allows us to obtain very accurate results, even with a very fewnumber of points in the numerical grid. Nevertheless, since all the required integrals are12ne–dimensional, we used a significant number of points ( ∼ ), regardless of whether itwas necessary.Having solved the Sturmian equation and generated the basis set S i ( ξ ), we can proceedto analyze convergence issues for the expansion (19) of the function U ( ξ ). In Table 2 thebasis size dependence of the energy value E , obtained by fixing the separation constant A = 0 . E s chosen to generate the Sturmian basis (20). Although it is anarbitrary energy, it is convenient to choose its value close to the true state energy. Assigning,for the first iteration, the arbitrary value E = − E = − . i.e. , we set E s = − . E (a.u.) ˜ E (a.u.)1 -1.0 -1.13 -1.1 -1.10266 -1.1024 -1.10263409 -1.1026 -1.1026346Reference -1.1026342Table 2: Convergence of the energy E in Eq. (8a) for fixed A = 0 . E obtained withan improved (recalculated) basis.Once the eigenvalues equation is solved, the eigenvectors of (24) provide the expansioncoefficients d i , which build the radial function U ( ξ ) through (19). The converged result isshown in Figure 5; the excited states will be discussed in the next section.The product of the angular and radial solutions Λ( η ) U ( ξ ) gives, up to the azimuthaldependence, the wavefunction which is best visualized by converting the prolates ( ξ, η, φ )into cartesian coordinates ( x, y, z ) through x = R (cid:112) (1 − η )( ξ −
1) cos( φ ) (34a) y = R (cid:112) (1 − η )( ξ −
1) sin( φ ) (34b) z = R ηξ . (34c)In the top left panel of Figure 6 we show the obtained ψ σ g for a fixed angle φ (for m = 0states, the results are symmetric respect to rotations over the z axis, and therefore, there isno dependence on the angle φ ). 13 ξ -0.500.511.5 U ( ξ ) σ g σ u σ g σ u Figure 5: The radial U ( ξ ) solutions for the four lower energy states of H +2 . In the ground state results presented above we have fixed, adopting the Born–Oppenheimerapproximation, the internuclear distance at R = 2 a.u. Calculations can be easily repeatedby varying R , and in each case, one obtains the total energy E tot ( R ) = E ( R ) + 1 R . (35)The radial Sturmian functions should be generated through Eq. (20) in which one modifies a = R ( Z + Z ) for each R . This option can be taken but we found it convenient to use aunique basis generated with a given value a s = R s ( Z + Z ); except for very high internucleardistances R , we simply took R s = 2 a.u.. In so doing, the use of the Sturmian equationfor the radial Schr¨odinger equation (8a) leads to a slightly modified Eq. (23) and thus thesupplementary matrix element ( a − a s ) [ N ] ij must be added to matrix N defined by (26a).Figure 7 presents the calculated total energy as a function of the internuclear distance. Theinset allows one to see a clear minimum at R = 1 . E tot = − . .We challenged our computational method with energy calculations considering very smallinternuclear distances R for which, in general, many numerical instabilities and errors arise.The energy values displayed in Table 3 demonstrate that our Sturmian method remainsrobust for decreasing distances R , even in the limit R →
0, for which the solution correspondsto the atomic ion He + with energy E He + = − Z / − −4−3−2−10123 z −4 −3 −2 −1 0 1 2 3 ψ g x −4−3−2−10123 z −4 −3 −2 −1 0 1 2 3 ψ g x −4−3−2−10123 z −4 −3 −2 −1 0 1 2 3 ψ g x −4−3−2−10123 z −4 −3 −2 −1 0 1 2 3 ψ g Figure 6: Wavefunction ψ σ g converted to cartesian coordinates, for the H +2 ground state,calculated at four different internuclear distances, moving from the molecular ion H +2 to theatomic ion He + : (top left) R = 2 . R = 1 . R = 0 . R = 0 .
008 a.u. To better appreciate the evolution we have renormalized thewavefunctions. 15 .5 1 1.5 2 2.5 3 3.5 4
R (a.u.) E t o t Figure 7: Total energy of the H +2 ground state as a function of the internuclear distance R .R (a.u.) E (a.u.)2 -1.10263401 -1.45178230.4 -1.8007540.1 -1.97825520.025 -1.99841130.008 -1.9998307He + -2.0Table 3: Ground state energy of the system H + H + e − , as a function of the internucleardistance R . d method for some excited states of H +2 By modifying the way the GSF basis functions are constructed, the GSF spectral methodallows one to obtain not only the ground state but also excited and continuum states. Tostart with, let us look at the first excited state 1 σ u . For the generation of the radial basis, it isnecessary to choose an arbitrary Sturmian energy as an external parameter. In a first, crude,approach we take the same energy obtained for the ground state calculation ( E s = − . p s = 1 . A . This value is introduced as a parameter into the radial equation, whosesolutions produce a new scaled energy value p . As shown in Table 4, a very precise result withsix significant figures is obtained after only eight iteration steps. However, as we discussedfor the ground state, we can make the whole calculation even better, choosing the Sturmian16nergy value from the last convergence step ( p s = 1 . E s = − . .Iteration p E (a.u.) ˜ p ˜ E (a.u.)0 1.485015 -1.10263 1.154791 -0.6667712 1.175548 -0.690957 1.155444 -0.6675254 1.155869 -0.668017 1.155451 -0.6675346 1.154793 -0.6667738 1.154791 -0.666771Reference p and energy E for the first H +2 excited state 1 σ u . The fifth columncorresponds to the energy ˜ E obtained with an improved (recalculated) basis.The same procedure is repeated for the generation of other excited states, such as 2 σ g and 2 σ u . In Table 5 the energy results obtained with our iterative method are displayed andcompare very favorably with the results obtained by Bian . Note that the latter coincide,up to the eighth decimal with those of Madsen and Peek .State A E (a.u.) ˜ E (a.u.) E (a.u.) Ref. σ g σ u -1.8689 -0.667 -0.6675338 -0.667534392 σ g σ u -1.69179 -0.25 -0.25535 -0.25541317Table 5: Parameter A and energies E of the lowest energy states of H +2 calculated with ouriterative GSF method. The fourth column corresponds to the energy ˜ E obtained with animproved (recalculated) basis. The last column reports the energy values found by Bian .The radial U ( ξ ) and the angular Λ( η ) solutions of the four lowest states of H +2 are shown,respectively, in Figures 4 and 5. The total wavefunctions for the excited states 1 σ u , 2 σ g and2 σ u are shown in Figure 8 as a function of the cartesian coordinates ( x, z ). We recall thatthe density is invariant under rotations around the z axis.17 −4−3−2−10123 z −4 −3 −2 −1 0 1 2 3 ψ u −101 x −6−5−4−3−2−1012345 z −6 −5 −4 −3 −2 −1 0 1 2 3 4 5 ψ g −101 x −6−5−4−3−2−1012345 z −6 −5 −4 −3 −2 −1 0 1 2 3 4 5 ψ u −101 Figure 8: Wavefunctions for the first excited states 1 σ u , 2 σ g and 2 σ u of H +2 .18 .3 Iterative d method for the asymmetric molecular ions HHe +2 and HLi +3 We apply now our GSF approach to other monoelectronic diatomic systems, such as theHHe +2 and HLi +3 molecular ions. For these heteronuclear ions, Z (cid:54) = Z and thus a (cid:54) = 0.In order to compare with other sample results published in the literature, we have keptthe internuclear distance fixed at R = 4 a.u. (for HHe +2 the equilibrium value is around R = 3 .
89 a.u.).The ground state wavefunctions of the heteronuclear molecular ions are shown in Figure 9.The distribution of the electron density is now clearly asymmetric, the logical shift towardsthe nucleus with larger charge being more evident as the Coulomb attraction increases.The shape of the wavefunction acquires more and more an atomic–like form centered onthe heavier nucleus with only relatively small values close to the hydrogen nucleus. Thesefeatures will obviously strongly depend on the internuclear distance, here fixed at R = 4a.u.. x −3−2−1012 z −3 −2 −1 0 1 2 ψ x −3−2−1012 z −3 −2 −1 0 1 2 ψ Figure 9: Ground state wavefunctions 1 σ of HHe +2 (left) and HLi +3 (right), assuming aninternuclear distance R = 4 a.u..Table 6 displays the calculated ground state energies, whose absolute value increasesapproximately as Z / Z the charge of the heavier nucleus. The efficiency of ourmethod can be appreciated by giving a few numbers of other methods. The results given19y Avery et al. were calculated with 10 basis elements (Coulomb Sturmian functions) foreach nucleus. Kereselidze et al. used 10 basis functions per nucleus (Coulomb Sturmian inprolate spheroidal coordinates). Xue–Bin Bian employed an imaginary–time–propagationmethod based on a Crank–Nicolson scheme to solve the separate equations, using 20 B–splines of order 7 to solve the radial equation, and 80 B–splines of order 7 for the angularpart. Campos et al. used 22 functions per coordinate. The aim of our calculation here wasnot to obtain very high accuracies, but rather to demonstrate that our simple and versatilemethod is computationally more efficient when compared to other approaches. If desired,we can achieve even better energy accuracies by improving the employed numerical methods(number of points or the finite differences order) or by tuning the generating potential as tooptimize the GSF basis set. E σ g H +2 (a.u.) E σ HHe +2 (a.u.) E σ HLi +3 (a.u.)This work -1.1026340 -2.2506056 -4.7501126Avery -1.10220 - -4.75011Kereselidze -1.102614 - -4.750111Bian -1.1026342 -2.2506054 -Campos - -2.2506054 -4.7501118Table 6: Ground state energies for the monoelectronic molecular ions: H +2 , assuming aninternuclear distance R = 2 a.u., and HHe +2 and HLi +3 , assuming an internuclear distance R = 4 a.u. . d method for the ground and excited statesof H +2Although we found excellent results with the iterative method, we wish to exploit the fulladvantages of the spectral method which allows one to obtain many states in one shot. Thedirect diagonalization of a 2 d Hamiltonian is generally very costly from the computationalpoint of view. Within the finite differences framework, and taking into account that everycoordinate has to be represented by hundreds of points, the matrix becomes huge and isintractable. A spectral method can reduce significantly the size of the Hamiltonian matrixto diagonalize, but computationally it still represents a hard task. Within the GSF method,the size of the matrices are reduced even more, since the appropriate physical behavior isexplicitly introduced in the basis set. In this manner, the numerical treatment is optimized.The use of expansion (29) on a two–dimensional basis S ij ( ξ, η ) transforms the Schr¨odingerequation into an equation (30) where all the derivatives have been removed and replaced bysimple expressions. Moreover, since the basis functions are optimized, the size of the basisis very small. For example, in our calculations, we introduced only 18 functions (3 angular S j ( η ) and 6 radial S i ( ξ )). Finally, the direct diagonalization of this small matrix produces,as a result, 18 states simultaneously without the need to perform separate iterations for eachstate. 20e have applied our GSF direct 2 d method to the benchmark ion H +2 , again taking R = 2 a.u.. With only one diagonalization we obtained the energy values displayed in Table7. They compare very well with the results of Madsen and Peek , following their statesnotation. We should point out that our aim here was to produce all the levels at the sametime without a focus on a single state. To generate the Sturmian basis we chose here theenergy value E s =-0.2 a.u. which is clearly quite different from the ground state energy; itis an acceptable compromise that leads to a good precision for the whole set of presentedmolecular states. The table shows that it is possible to obtain excellent results, in particularfor the lower states, at a rather small computational cost. If one wishes to improve theenergy accuracy for one particular state, a different Sturmian energy E s closer to this stateenergy should be chosen, as was shown in the 1 d method. Since the generation of a newSturmian basis requires one–dimensional calculations and the 2 d matrix only has a few dozenof elements, this further optimization procedure is rather inexpensive.State E (a.u.) E (a.u.) Ref. Sσ g -1.102630 -1.102634212 P σ u -0.66753431 -0.667534392 Sσ g -0.360863 -0.360864883 P σ u -0.25541312 -0.255413173 Dσ g -0.2357775 -0.235777633 Sσ g -0.1776 -0.177681054 P σ u -0.133 -0.13731293Table 7: Energies of seven energy states of H +2 , obtained with the GSF direct 2 d method.The third column indicates the results of Madsen and Peek . Both were obtained for afixed internuclear distance R = 2 a.u.. The spectral method, based on Generalized Sturmian Functions, has been here extended, toallow its use in prolate spheroidal coordinates which should provide, in principle, the mosteffective framework to treat diatomic molecular systems. We developed and implemented twodifferent numerical methods for the calculation of the molecular structure of monoelectronicmolecular ions.The first one consists in separating the Schr¨odinger equation in one angular and oneradial equations, coupled through the energy and a coupling parameter. The equations aresolved alternately, fixing the energy and the coupling parameter in each case, and after a fewiterations, these parameters converged to the final values. The advantage of using GSF istwofold. On the one hand, it allows one the replacement of most of the Hamiltonian calcu-lations by a simple expression thus substantially reducing the complexity of the calculationat any iteration step. On the other hand, the GSF method is based in the valuable propertythat the right boundary conditions are enforced onto the basis functions. Therefore, the21ize of the basis is minimal, turning the method in a very efficient procedure that producesground and excited states of high quality.The second method also uses GSF, and the angular and radial basis sets are generatedin the same way as in the first one. Then, a two–dimensional basis set is constructed, andthe Schr¨odinger equation solution becomes a 2 d generalized eigenvalues problem. Since thebasis elements have the correct boundary conditions, the size of the basis is very small, andthe diagonalization is not a costly procedure. This direct 2 d method does not require anyiteration and a single calculation yields – simultaneously – many molecular states. Verygood results can be obtained already with small basis size.As a first step towards the extension of the GSF method to diatomic molecules, we havepresented here an investigation of molecular ions having only one electron. We calculatedthe ground and excited states of the molecular hydrogen ion H +2 , in excellent agreementwith benchmark results (7 significant figures in the case of the ground state). We alsostudied heteronuclear molecular ions, like HHe +2 and HLi +3 , with again excellent results.The method proved to be robust over a wide range of internuclear distances R , including inthe notoriously difficult atomic limit.The whole numerical investigation gives us confidence in our implementation of the GSFmethod in prolate spheroidal coordinates, as to contemplate exploring the continuous partof of the spectrum. As demonstrated for atomic systems, the advantages of the GSF spectralmethod are more evident in the treatment of collision problems. In this case, the continuumSturmian basis elements are generated with a positive energy parameter E s and one imposesappropriate scattering boundary conditions. As a consequence, the basis needs to solve theSchr¨odinger equation only in the interaction region. Scattering problems involving one ortwo electrons in the continuum can then be treated efficiently with compact bases .The same arguments apply to diatomic molecular systems, and we plan to extend the presentinvestigation in prolate spheroidal coordinates to scattering problems such as single or doubleionization by photon or electron impact. First we will examine the single continuum bystudying the single photoionization of the benchmark one–electron molecular ion H +2 ; then,we will move to the more challenging two–electron correlated case, by investigating singleand double ionization processes on H and on quasi two–electron targets like N as done forexample in Ref. . DM gratefully acknowledge the financial support from the following Argentine institutions:Consejo Nacional de Investigaciones Cient´ıficas y T´ecnicas (CONICET), PIP 11220130100607,Agencia Nacional de Promoci´on Cient´ıfica y Tecnol´ogica (ANPCyT) PICT–2017–2945, andUniversidad de Buenos Aires UBACyT 20020170100727BA.22
REFERENCESReferences [1] Burrau, Ø., Kgl. Danske, Videnskab. Selskab. Mat. Fys. Medd., , 14 (1927).[2] Hylleraas, E. A., Z. Phys., , 739 (1931).[3] Jaff´e, G., Z. Phys., , 535 (1934).[4] Bransden B. H., and Joachain, C. J., The Physics of Atoms and Molecules (LongamnScientific and Technical: Harlow, UK, 1983).[5] Carrington, A., McNab I. R., and Montgomerie, C.A., J. Phys. B, , 3551 (1989).[6] Bian, X. B., Phys. Rev. A, , 033403 (2014).[7] Tamaz Kereselidze, Irakli Noselidze and Alexander Devdariani, J. Phys. B: At. Mol.Opt. Phys., , 105003 (2019).[8] Serov V. V., Joulakian B. B., Pavlov D. V., Puzynin I. V., and Vinitsky S. I., Phys.Rev. A, , 062708 (2002).[9] Chuluunbaatar O., Joulakian B. B., Tsookhuu K., and Vinitsky S. I., J. Phys. B: At.Mol. Opt. Phys., , 2607 (2004).[10] Serov V. V., Joulakian B. B., Derbov V. L., and Vinitsky S. I., J. Phys. B: At. Mol.Opt. Phys., , 2765 (2005).[11] Chuluunbaatar O., Joulakian B. B., Puzynin I. V., Tsookhuu Kh., and Vinitsky S. I.,J. Phys. B: At. Mol. Opt. Phys., , 015204 (2008).[12] Tao L., McCurdy C. W., and Rescigno T. N., Phys. Rev. A, , 012719 (2009).[13] Serov V. V., and Joulakian B. B., Phys. Rev. A, , 062713 (2009).[14] Foster M., Colgan J., Al–Hagan O., Peacher J. L., Madison D. H., and Pindzola M. S.,Phys. Rev. A, , 062707 (2007).[15] Mitnik, D. M., Colavecchia, F. D., Gasaneo, G., and Randazzo, J. M., Comp. Phys.Comm., , 1145 (2011).[16] Gasaneo, G., Ancarani, L. U., Mitnik, D. M., Randazzo, J. M., Frapiccini, A. L., andColavecchia, F. D., Adv. Quantum Chem., , 153 (2013).[17] Randazzo, J. M., Ancarani, L. U., Gasaneo, G., Frapiccini, A. L., and Colavecchia, F.D., Phys. Rev. A, , 042520 (2010).[18] Randazzo, J. M., Mitnik, D., Gasaneo, G., Ancarani, L. U., and Colavecchia, F.D., Eur.Phys. J. D, , 189 (2015). 2319] Ambrosio, M. J., Mitnik, D. M., Dorn, A., Ancarani, L. U., and Gasaneo, G., Phys.Rev. A, , 032705 (2016).[20] Ambrosio, M. J., Ancarani, L. U., Gomez, A. I., Gaggioli, E. L., Mitnik, D. M., andGasaneo, G, Eur. Phys. J. D, , 127 (2017).[21] Granados–Castro, C., “Application of Generalized Sturmian Basis Functions to Molec-ular Systems”, PhD thesis, Universit´e de Lorraine, (2016).[22] Granados–Castro, C. M., Ancarani, L. U., Gasaneo, G., and Mitnik, D. M., Adv. Quan-tum Chem., , 3 (2016).[23] Granados–Castro, C. M., and Ancarani, L. U., Eur. Phys. J. D, , 65 (2017).[24] Ali, E., Granados, C., Sakaamini, A., Harvey, M., Ancarani, L. U., Murray, A. J.,Dogan, M., Ning C., Colgan, J., and Madison, D., J. Chem. Phys., , 194302 (2019).[25] Edmonds, A. R. Angular Momentum in Quantum Mechanics (Princeton UniversityPress: Princeton, NJ, 1957).[26] Scott, T. C., Aubert-Fr´econ, M., and Grotendorst, J., Chem. Phys., , 851 (1970).[29] Madsen, M. M., and Peek, J. M., At. Data and Nucl. Data Tables, , 171 (1971).[30] Avery, J., and Avery, J., J. Phys. Chem. A, , 14565 (2009).[31] Kereselidze, T., Chkadua, G., and Defrance, P., Molec. Phys., textbf113:22, 3471 (2015).[32] Campos, J. A., Nascimento, D. L., Cavalcante, D. T., Fonseca, A. L. A., and Nunes, A.O. C., Int. J. Quantum Chem., , 2587 (2006).[33] Chuluunbaatar O., Gusev A. A., and Joulakian B. B., J. Phys. B: At. Mol. Opt. Phys., , 015205 (2012).[34] Bulychev A. A., Chuluunbaatar O., Gusev A. A., and Joulakian B., J. Phys. B: At.Mol. Opt. Phys.,46