The one-electron Schrödinger equation in prolate spheroidal coordinates . 35

Consider a diatomic molecule within the Born-Oppenheimer approximation. Assume that the one-electron Hamiltonian ˆh possesses rotational symmetry with respect to the internuclear axis and can therefore be written as

h =ˆ −1

2∇²+V(ξ, η). (2.16)

The explicit form of V(ξ, η) which depends on the interaction between an electron and ionic cores is not further specified at this point. Solving the one-electron Schrödinger equation (OESE) ˆhψ=ψyields the orbitalsψ. Due to rotational symmetry the orbitals

ψ can be specified by the quantum numbers λ= 0,1, . . . and m=±λ. If the function V(ξ, η) fulfills the conditionV(ξ,−η) =V(ξ, η) (for example, for a homonuclear molecule with equivalent cores), the Hamiltonian ˆh also possesses inversion symmetry with respect to the midpoint between the nuclei,

ξ→ξ, η→ −η, φ→φ+π. (2.17) The eigenvalue℘of the inversion operator is equal to 1 forgeradeand to−1 forungerade states. The complete set of orbital quantum numbers is specified by γ ≡ {λ, m,[℘]}

(square brackets are used for optional quantum numbers). Since the energy does not depend onm,|γ| ≡ {λ,[℘]} is introduced. For example, in the case of H⁺₂ orbitals with symmetryσ_g and π_u correspond to|γ|={0,1} and {1,−1}, respectively.

Using the rotational symmetry of the OESE,

ˆhψ_ν^γ=^|γ|_ν ψ^γ_ν, (2.18) the orbitalsψ^γν can be written as

ψ_ν^γ(ξ, η, φ) = 1

√2πψ^|γ|_ν (ξ, η)e^imφ. (2.19) If the potentialV(ξ, η) has the form

V(ξ, η) =− 2

R²(ξ²−η²){Z_ξ(ξ) +Zη(η)} (2.20) with functionsZξ(ξ) and Zη(η) depending on a single variable, the OESE is separable.

As an example, consider the Hamiltonian ˆh that describes the motion of an electron in the Coulomb field of two nuclei with chargeZ_A and Z_B,

V =−ZA

r_A −ZB

r_B. (2.21)

Then, substituting (2.8) into (2.21) and transforming to the form (2.20), one obtains Z_ξ(ξ) = (Z_B+Z_A)R ξ, Z_η(η) = (Z_B−Z_A)R η. (2.22) In case of separability of the OESE the orbitals can be written as

ψ^γ_ν(ξ, η, φ) = N_ν^|γ|Ξ^|γ|_ν (ξ) Π^|γ|_ν (η) e^imφ, (2.23)

whereN_ν^|γ| is a normalization factor. Substituting (2.11), (2.16), (2.20), and (2.23), the equation (2.18) can be transformed into the set of eigenvalue equations for Ξ^|γ|_ν (ξ) and Π^|γ|_ν (η)

− ∂

∂ξ

(ξ²−1)∂

∂ξ

+ λ²

ξ²−1 −Z_ξ(ξ)−R² 2 ^|γ|_ν ξ²

Ξ^|γ|_ν (ξ) =A^ξ(^|γ|_ν )Ξ^|γ|_ν (ξ)

−∂

∂η

(1−η²) ∂

∂η

+ λ²

1−η² −Zη(η) +R² 2 ^|γ|_ν η²

Π^|γ|_ν (η) =A^η(^|γ|_ν )Π^|γ|_ν (η)

(2.24)

where the eigenvaluesA^ξ(^|γ|ν ) andA^η(^|γ|ν ) connect Eqs.(2.24) via the condition

A^ξ(^|γ|_ν ) +A^η(^|γ|_ν ) = 0. (2.25)

In order to solve the OESE two sets ofB splines are introduced. The first set comprises n_ξB splinesB_r(ξ)(r= 1..n_ξ) of orderk_ξfor variableξ. The knot sequence{ξ_i^b}is chosen according to

1 =ξ₁^b =...=ξ^b_kξ < ξ_k^b_ξ+1< ... < ξ_n^b_ξ+1 =...=ξ_n^b_ξ+kξ =ξ_max, (2.26) where the basis-set parameter ξ_max defines the size of the elliptical “box” in which the OESE is solved. Similarly, the second set comprises n_η B splines B_r(η)(r = 1..n_η) of order kη for variableη, but now using the knot sequence{η^b_i}according to

−1 =η₁^b =...=η_k^b_η < η^b_kη+1 < ... < η_n^b_η+1=...=η_n^b_η+kη = 1. (2.27) Instead of using the B splines directly, the (symmetry dependent) functions

X_α^|γ|(ξ) = (ξ²−1)^λ2 B_α(ξ) and Y_β^|γ|(η) = (1−η²)^λ2 B_β(η) (2.28) are introduced to correctly describe the square-root singularities at ξ = 1 and η =

±1. With the aid of all possible products of Xα^|γ| and Y_β^|γ|, the three-dimensional basis functions

U_αβ^γ (ξ, η, φ) = 1

√2πX_α^|γ|(ξ)Y_β^|γ|(η)e^imφ. (2.29) are constructed. Similar to Sec. 1.3, the wave functions are enforced to vanish at the box boundaryξmax. Since onlyBnξ(ξ) is nonzero atξmax, its coefficients should be zero.

This is equivalent to the simple omission of Xn^|γ|ξ(ξ) and the use of only ˜n_ξ = n_ξ −1 functions Xα^|γ|(ξ), α= 1..˜n_ξ.

While it is in general required to use ˜nη =nη basis functionsY_β^|γ|(η), Hamiltonians with inversion symmetry allow the reduction of the number of basis functions to ˜nη =nη/2.

With the choice of the knot sequence{η_i^b} according to η_r^b=−η^b_n

η+kη+1−r for r= 1.. n_η+k_η (2.30) with even value of nη, the identity Bβ(η) = Bnη+1−β(−η) is fulfilled. Instead of the definition given in Eq. (2.28) the symmetry-adapted functions are defined as

Y_β^|γ|(η) = (1−η²)^λ2{B_β(η) + (−1)^λ℘B_nη+1−β(η)}. (2.31) This leads to correct inversion symmetry ofU_αβ^γ in Eq.(2.29).

After replacement of the two indices α = 1..˜n_ξ and β = 1..˜nη by the single index i= (α−1)˜n_η+β, the OESE is projected onto a set of ˜N = ˜n_ξn˜_η linear independent basis functions{U_i^γ}. This yields the generalized matrix eigenvalue problem

h¯

|γ|C~^|γ|=^|γ|S

¯

|γ|C~^|γ| (2.32)

with

h¯

|γ|

ij = Z Z Z

V

dv U_i^γ∗ˆhU_j^γ and S

¯

|γ|

ij = Z Z Z

V

dv U_i^γ∗U_j^γ (2.33) where dv is given by Eq. (2.9) and V is the integration volume corresponding to the ellipsoid bounded by ξ_max. Due to the local support of B splines the matrices h

¯

|γ|

ij

andS

¯

|γ|

ij are banded. This property is efficiently considered when the diagonalization is performed e.g. with LAPACK subroutine DSBGVX that also accounts for the symmetry of the matrices. The integration over angle φ in Eq. (2.33) can easily be performed.

The remaining two-dimensional integrals have to be solved numerically (for example, by Gaussian quadrature). If the potential V(ξ, η) is of the form (2.20), it is, however, possible to write the integrals as linear combinations of products of one-dimensional integrals. Each of them is a sum of integrals over intervals between neighbor knot points. In the case of a Coulomb field (2.22) the integrands are polynomials at most of order 2(k_ξ+λ) or 2(kη +λ) and Gaussian quadrature of order k_ξ+λor kη +λyields exact results.

The diagonalization subroutine provides a set of eigenvalues ^|γ|ν , sorted in ascending order (ν = 1..N˜), and (optionally) the corresponding coefficients C_ν,i^|γ|. The orbitals

(2.18) are thus given as

ψ^γ_ν(ξ, η, φ) =

N˜

X

i=1

C_ν,i^|γ|U_i^γ(ξ, η, φ). (2.34)

Even if the potential V(ξ, η) allows separability, the described procedure does not pro-vide orbitals of the form given in Eq.(2.23), because every coefficient in Eq.(2.34) belongs to a pair of basis functions with respect toξ andη. This does not only lead to a longer expansion than necessary (˜n_ξn˜η instead of ˜n_ξ+ ˜nη), but, even more importantly, does not allow the significant simplification in the evaluation of one- and two-electron inte-grals that is discussed for the latter ones below (see Eq.(2.51)). Therefore, the following procedure for obtaining the eigenvectors is implemented for the case of separable Hamil-tonians. Only the eigenvalues ^|γ|_ν of Eq.(2.32) are calculated as described above. Then for every eigenvalue ^|γ|ν the Eqs.(2.24) are projected onto the sets of basis functions {X_α^|γ|} and {Y_β^|γ|}. The resulting two generalized matrix eigenvalue problems are again solved (including calculation of the eigenvectors) using subroutine DSBGVX. Out of the obtained eigenvalue sets {A^ξ(^|γ|ν )} and {A^η(^|γ|ν )} that pair A^ξαν(^|γ|ν ) and A^η_β

ν(^|γ|ν ) is chosen, whose sum is close to zero (see Eq.(2.25). The eigenvectors {C_ν,α^|γ|} and {C_ν,β^|γ|} belonging to this pair of eigenvalues define the orbital with energy ^|γ|_ν in the form of Eq.(2.23) with

Ξ^|γ|_ν (ξ) =

˜ nξ

X

α=1

C_να^|γ| X_α^|γ|(ξ) and Π^|γ|_ν (η) =

˜ nη

X

β=1

C_νβ^|γ|Y_β^|γ|(η). (2.35)

Note, that if the functionsY_β^|γ|(η) are defined according to Eq. (2.31) the function Π^|γ|ν (η) satisfies

Π^|γ|_ν (−η) = (−1)^λ℘Π^|γ|_ν (η). (2.36) It is useful to specify the orbitals ψ_ν^γ in Eq. (2.23) by the numbersN_ξ and N_η of nodes of the functions Ξ^|γ|_ν (ξ) and Π^|γ|ν (η). In fact, the set{N_ξ, Nη, m}can be used to uniquely identify an orbital, alternative to the set {ν, γ}, and can be also treated as the set of quantum numbers (see Fig. 2.2 for an example). The physical meaning of Nξ and Nη

can be well understood in the united-atom-limit case. Consider the nuclei with changes Z_A=Z_B=Z/2 and R→0. Applying Eq. (2.7) to Eq. (2.24) it can be derived that

Ξ^|γ|_ν (ξ) ≈ R_nl(r), n =N_ξ+N_η+λ+ 1,

Π^|γ|_ν (η)e^imφ ≈ Y_lm(ˆr), l =N_η+λ, (2.37)

1 10

provided both Ξ^|γ|_ν and Π^|γ|_ν are properly normalized. HereR_nl(r) is the radial wavefunc-tion of hydrogenic atom with the charge Z, with the principal quantum number nand the angular momentum l, while Y_lm(ˆr) is a spherical harmonic. Thus, for R → 0 the orbital (2.23) with the set{N_ξ, N_η, m} is equivalent to the hydrogenic orbital with the quantum numbers{n, l, m}. (Because of the box boundary condition this is not always the case for numerically calculated orbitals).

Using these relations it is possible to select configurations for the subsequent CI calcu-lations in a more intuitive way (as is discussed below in Sec. 3.1). Even in the case of non-separability (as long as it is weak) it is possible to characterize the orbitals by two numbers that are related to the nodal structure. Since Eqs.(2.24) are the Sturm-Liouville problem, the number of nodes can be determined from the indicesαν andβν. For Ξ^|γ|_ν (ξ) one findsN_ν^ξ=α_ν−1. For Π^|γ|_ν (η) one findsN_ν^η =β_ν−1, if the basis functions Y_β^|γ|(η) are chosen according to Eq.(2.28), andNν^η = 2β_ν −(3 + (−1)^λ℘)/2, if the Y_β^|γ|(η) are chosen according to Eq.(2.31).

For the purpose of illustration, Fig. 2.3 shows electronic potential energiesE(R) of the H⁺₂ and HeH⁺⁺ molecules obtained with the present method. In order to present the

0 2 4 6 8 10 12 14 16

Scaled internuclear distance √− R (a.u.)

Figure 2.3.: Low-lying bound states of H⁺₂ (a,b) and HeH⁺⁺(c,d) molecules. a,c) Molecu-lar potential energyE(R)as a function of the internuclear distanceR. b,d) Effective principle quantum number n^∗ = (ZA+ZB)/√

2Ebind, where Ebind being the binding energy, versus the square root of the internuclear distance R.

obtained energies in a wide region of internuclear separations, Figs. 2.3 b and 2.3 d show the effective principle quantum number n^∗ = (Z_A+Z_B)/√

2E_bind which depends on the binding energy E_bind = E(R) −(ZAZB)/R of the molecular ion. In addition, a scaled abscissa axis is used. Such a representation allows one to demonstrate how the molecular energies converge to both the united-atom limit (where the system reduces to a one-electron atom with the nuclear chargeZ_A+Z_B) and the dissociation limit (where the system reduces to two non-interacting atoms).

Since all orbitals are enforced to be zero at ξmax, box discretization is employed in the present approach. Consequently, all orbitals are normalized to one, independently on whether they represent bound or continuum solutions. In one-electron calculations involving the continuum (for example, photoionization of H⁺₂) it is, however, necessary to adopt energy normalization for the continuum orbitals. This can be achieved either

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 Energy (a.u.)

10^-18 10^-16 10^-14 10^-12 10^-10 10^-8 10^-6 10^-4

Oscillator strength f (a.u.) ¹

3 5 7 9

Figure 2.4.: Oscillator strengths from the ground state 1σg of H⁺₂ to computedσu orbitals with differentN_η.

by fitting the orbitals (close to the box boundary) to the known asymptotic behavior or using the density of states.

Since Nη is a good quantum number, it can be used to distinguish between different continuum channels, instead of the angular momentumlin the case of a single-electron atom. This possibility is demonstrated in Fig. 2.4 where the oscillator strengthsf,

f_ij = 2(E_j−E_i)|hi|z|ji|², (2.38) for transitions from the ground state 1σg of H⁺₂ to computed σu orbitals with different Nη are shown. The electronic-structure calculation was performed for the internuclear distanceR= 2.0a₀ using the box sizeξ_max= 700 and 500B splines for variableη. Due to symmetry,Nηmay possess only an odd value, and the oscillator strengthsf to orbitals withN_η = 1 and 3 are evidently stronger (Fig. 2.4) compared to those with higher N_η. As one can see, the oscillator strengths f belonging to the same Nη can be joined by a smooth line (in fact, the density of points is so high that individual points are often indistinguishable), and thus the interpolation procedure can be efficiently performed to

calculate the contribution of continuum channels at an arbitrary energy 0 < E < E_cr where the value of the critical energy Ecr = 2 a.u. can be inferred from the figure.

Substituting relation (2.7) in Eq. (1.29) yields the following expression for an estimation of the critical energy

Ecr≈3 2s

Rξ_max 2

(2.39) which yields Ecr = 1.5 a.u. for the discussed case. In order to calculate the density of states one has to consider separately orbitals belonging to different N_η. However, if the orbitals are used in a CI calculation, they should not be energy renormalized, even if one is interested in continuum properties of the two-electron system. In fact, the orbitals should in this case only be seen as basis functions. Therefore, even the non-physical pseudo-states necessarily obtained in a box-discretized calculation (due to the infinite Rydberg series represented by a finite basis) are important ingredients for a subsequent CI calculation.

2.4. Two-electron basis set and configuration-interaction approach

As a basis set for solving the two-electron Schrödinger equation, H(1,ˆ 2)Ψ(1,2) =

h(1) + ˆˆ h(2) + 1 r12

+VAB(R)

Ψ(1,2) =EΨ(1,2), (2.40) linear combinations of the orbitals [given either in Eqs. (2.23) and (2.35) or in Eqs. (2.34) and (2.29)] are used. The potential V_AB(R) is either the Coulomb interaction Z_AZ_B/R between the nuclei or describes the interaction between two cores. Let Γ specify the full set of quantum numbers, Γ≡ {Λ, M, S,[P,P^Σ]}, where Λ = 0,1, . . . is the absolute value of the component of the total angular momentum along the internuclear axis (Σ,Π, . . .), M =±Λ and the total spin S= 0,1. The optional quantum numbersP andP^Σ specify the parity with respect to inversion symmetry (gerade or ungerade, cf. Eq. (2.17)) and, for Σ states only, with respect to a reflection at a plane through the nuclei,

ξ_i →ξ_i, η_i→η_i, φ_i → −φ_i. (2.41) The last transformation is equivalent to a complex conjugation of the wave function, if the normalization of the latter is chosen in such a way, that only its angular part is complex.

Normalized and fully symmetry-adapted configurations Υ^Γ_i are used in the CI calcula-tions. The configurations are build with the aid of products of two orbitals |νγ ν¯γ¯i ≡ ψν^γ(ξ₁, η₁, φ₁)ψ_ν^¯γ_¯(ξ₂, η₂, φ₂) with ¯γ ≡ {λ,¯ m,¯ [ ¯℘]}. For a given configuration Υ^Γ_i ≡[νγ¯ν¯γ]

the orbitalsνγand ¯νγ¯satisfy the conditionsm+ ¯m=M and ℘℘¯=P. For non-Σ states (Λ6= 0) particle-exchange symmetry is considered using

[νγν¯¯γ]_Λ6=0 = |νγν¯¯γi+ (−1)^S|¯νγ νγi¯

p2(1 +δνν¯δγ¯γ) . (2.42) For Σ states (Λ = 0) inclusion of the additional reflection symmetry leads to (note the misprint in [61])

[νγν¯¯γ]_Λ=0 = |νγν¯γi¯ + (−1)^S|¯νγ νγi¯

+P^Σ |νγν¯γ¯i+ (−1)^S|¯ν¯γ νγi∗

q

4(1 +δ_ν¯νδ_|γ||¯γ|)(1 +δ_0mδ_{0 ¯m})

(2.43)

where it was used that two different orbitals with the same ν and |γ| can be obtained from each other by complex conjugation. Using the relation

|νγν¯γ¯i^∗ ≡ |ν{λ, m,[℘]}ν{¯ λ,¯ m,¯ [ ¯℘]}i^∗

=|ν{λ,−m,[℘]}ν{¯ ¯λ,−m,¯ [ ¯℘]}i ≡ ||νγ^∗ ν¯γ¯^∗i, (2.44) configuration (2.43) can finally be represented as a linear combination of products of two orbitals. Clearly, the choice ofνγ and ¯νγ¯ has always to be done in such a way that a double occurrence of the same configuration is prevented. The two-electron problem

H(1,ˆ 2)Ψ^Γ_µ(1,2) =E_µ^ΓΨ^Γ_µ(1,2) (2.45) is solved by an expansion of the wave functions as linear combinations ofN_Γ configura-tions,

Ψ^Γ_µ(1,2) =

NΓ

X

i=1

C_µi^ΓΥ^Γ_i(1,2). (2.46) Due to the orthonormality of the set{Υ^Γ_i}the ordinary matrix eigenvalue problem

H¯

ΓC^Γ=E^ΓC^Γ with H

¯

Γ ij =

Z

V1

Z

V2

dv₁dv₂Υ^Γ∗_i (1,2) ˆH(1,2)Υ^Γ_j(1,2) (2.47) is obtained. The coefficients{C_µi^Γ}and energies E_µ^Γ are calculated using LAPACK sub-routine DSPEVX. The CI matrix elementsH

¯

Γ

ij comprise one- and two-electron integrals.

The evaluation of the one-electron integrals is trivial, because the orbitals are solutions

of the one-electron Hamiltonian. Thus one finds

hνγν¯γ|¯ˆh(1) + ˆh(2)|ν⁰γ⁰ν¯⁰¯γ⁰i= (^|γ|_ν +^|¯_ν¯^γ|)δ_νν⁰δ_γγ⁰δ_ν¯¯ν⁰δ_γ¯¯γ⁰. (2.48) Using the von Neumann expansion (2.12) the two-electron integrals can be written after some further simplifications [62] as

hνγν¯¯γ| 1 In case of separability of the OESE the integral functionF_νν^|γ||γ0`^0|ρ(ξ) reduces to a product of one-dimensional integrals,

For Hamiltonians with inversion symmetry the functions F_νν^|γ||γ0`<sup>0|ρ</sup>(ξ) are equal to zero, if the condition ℘℘<sup>0</sup> = (−1)<sup>`+1</sup> is fulfilled. Indeed, from P<sub>`</sub><sup>ρ</sup>(−η) = (−1)<sup>`+ρ</sup>P_`^ρ(η) and (2.36) follows

P<sub>`</sub><sup>ρ</sup>(−η)Π<sup>|γ|</sup><sub>ν</sub> (−η)Π<sup>|γ</sup><sub>ν</sub>0<sup>0|</sup>(−η) = (−1)<sup>`+λ+λ0+ρ</sup>℘℘⁰P_`^ρ(η)Π^|γ|_ν (η)Π^|γ_ν0^0|(η), (2.52) where (−1)^λ+λ0+ρ= 1 for anymandm⁰, as can be easily shown from the definition ofρ.

Thus, in (2.51) the integrals with respect to ηare equal to zero for℘℘⁰ = (−1)<sup>`+1</sup>, since in this case their integrands are odd functions with respect to η. If functionF<sub>νν</sub><sup>|γ||γ</sup>0`^0|ρ(ξ) is non-zero, the ratioF_νν^|γ||γ0`<sup>0|ρ</sup>(ξ)/(ξ−1)<sup>t</sup>witht= (ρ+λ+λ<sup>0</sup>)/2 + 1 remains finite even in the limitξ→1. Since (ξ<sup>2</sup>−1)[ ˆP<sub>`</sub>^ρ(ξ)]^{2 ξ→1}−→B(ξ−1)^ρ+1 whereB is a non-zero constant,

one obtains

F_νν^|γ||γ0`<sup>0|ρ</sup>(ξ)F<sub>ν¯</sub><sup>|¯γ||¯</sup><sub>ν¯</sub>0`^γ0|ρ(ξ) (ξ²−1)[ ˆP_`^ρ(ξ)]²

−→ξ→1D(ξ−1)^{λ+λ0+¯λ+¯λ0+1}, with D= const. (2.53)

Since (ξ²−1)[ ˆP<sub>`</sub><sup>ρ</sup>(ξ)]<sup>2</sup> 6= 0 at any point ξ ∈(1, ξmax], the integrand has no singularities withinξ∈[1, ξ<sub>max</sub>] and Gaussian quadrature can be used to calculate this integral. It is important to mention that the variety ofF<sub>νν</sub><sup>|γ||γ</sup>0`^0|ρ(ξ) functions is much smaller than the number of two-electron integrals. The storage of the values of the functionsF_νν^|γ||γ0`^0|ρ(ξ) at Gaussian quadrature roots allows the simultaneous calculation of all two-electron integrals. This speeds up the calculations significantly.

Im Dokument Ionization of molecular hydrogen in ultrashort intense laser pulses (Seite 47-58)