Quadrupole transitions in the bound rotational-vibrational spectrum of the hydrogen molecular ion

# Quadrupole transitions in the bound rotational-vibrational spectrum of the hydrogen molecular ion

Horacio Olivares Pilón and Daniel Baye Physique Quantique C.P. 165/82, and Physique Nucléaire Théorique et Physique Mathématique, C.P. 229, Université Libre de Bruxelles (ULB), B 1050 Brussels, Belgium
###### Abstract

The three-body Schrödinger equation of the H hydrogen molecular ion with Coulomb potentials is solved in perimetric coordinates using the Lagrange-mesh method. The Lagrange-mesh method is an approximate variational calculation with variational accuracy and the simplicity of a calculation on a mesh. Energies and wave functions of up to four of the lowest vibrational bound or quasibound states for total orbital momenta from 0 to 40 are calculated. The obtained energies have an accuracy varying from about 13 digits for the lowest vibrational state to at least 9 digits for the third vibrational excited state. With the corresponding wave functions, a simple calculation using the associated Gauss quadrature provides accurate quadrupole transition probabilities per time unit between those states over the whole rotational bands. Extensive results are presented with six significant figures.

###### pacs:
31.15.Ag, 31.15.Ac, 02.70.Hm, 02.70.Jn
: J. Phys. B: At. Mol. Opt. Phys.

August 2, 2019

## 1 Introduction

The H bound rotational-vibrational spectrum possesses about 420 bound states corresponding to the electronic configuration as well as about 60 narrow quasibound levels. For this simple three-body system, the Schrödinger equation involving Coulomb potentials can not be solved exactly but it is possible to reach a high accuracy for both energies and wave functions. A calculation of the energies of most rotational-vibrational bound states and some quasi-bound levels has been performed in 1993 with 11-digit accuracy by Moss [1]. He also calculated energies for the few extended bound states. Since then, several calculations have reached a higher accuracy, but always for a limited number of bound states. The energy of the ground state, with total orbital momentum and positive parity, has been improved in a number of papers [2, 3, 4, 5, 6, 7, 8] to reach an accuracy around or beyond 30 digits [7, 8]. The energy of the first excited state has also been improved [6, 8]. The lowest vibrational state has been considered by several authors [5, 6, 7, 8]. The and excited vibrational states have been studied by Hilico et al[3] (see also [9]). Energies of the full ground-state rotational band were determined with 12-digit accuracy by Hesse and Baye [10] and energies up to with higher accuracy were obtained by Yan and Zhang [11]. Some of the quoted works also provide accurate information on mean radii [10, 8] or quadrupole moments [10].

In opposition to this large number of accurate studies, very few studies concern the transition probabilities between the rotational-vibrational bound states of H. There are several reasons for this. As the electric dipole transitions are forbidden by the symmetry the two protons, the more complicated electric quadrupole transitions are the dominant mode of decay. Moreover, few existing studies provide the necessary wave functions. Since the pioneering work of Bates and Poots [12], a systematic study of all transitions for states up to within the Born-Oppenheimer approximation has been published with two significant figures [13]. Here we present accurate E2 transition probabilities without Born-Oppenheimer approximation. They are obtained from three-body wave functions calculated with the Lagrange-mesh method in perimetric coordinates [14, 15, 10], with which the calculation is particularly simple and very precise.

The Lagrange-mesh method is an approximate variational calculation using a basis of Lagrange functions and the associated Gauss quadrature. It has the high accuracy of a variational approximation and the simplicity of a calculation on a mesh. Lagrange functions are orthonormal infinitely differentiable functions that vanish at all points of this mesh, except one. Used as a variational basis in a quantum-mechanical calculation, the Lagrange functions lead to a simple approximation when matrix elements are calculated with the associated Gauss quadrature. The variational equations take the form of mesh equations with a diagonal representation of the potential only depending on values of this potential at the mesh points [16, 17]. The most striking property of the Lagrange-mesh method is that, in spite of its simplicity, the obtained energies and wave functions can be as accurate with the Gauss quadrature approximation as in the original variational method with an exact calculation of the matrix elements [18, 17]. The accuracy of the lowest energies exceeds by far the accuracy of the Gauss quadrature for the individual matrix elements. The Lagrange-mesh method allows very accurate calculations not only in simple quantum-mechanical problems but also in various more complicated applications in atomic [19, 20, 21], molecular [14, 15, 10, 22, 23] and nuclear [24, 25, 26] physics.

In the H case, the Lagrange-mesh method is applied in perimetric coordinates, i.e. three angles describing the orientation of the plane containing the particles and three linear combinations of the distances between them [27]. The dependence on the three Euler angles is treated analytically [15, 10]. The three perimetric coordinates vary from zero to infinity and can easily be discretized on a three-dimensional Lagrange-Laguerre mesh [14]. An additional advantage is that the resulting matrix is rather sparse. The Lagrange-mesh method also provides analytical approximations for the wave functions that lead to very simple expressions for a number of matrix elements when used with the corresponding Gauss-Laguerre quadrature.

In section 2, the basic expressions for the transition probabilities are recalled. The E1 and E2 operators are expressed in perimetric coordinates. Some definitions about Lagrange functions lead to Lagrange-mesh expressions for the transition matrix elements. In section 3, energies are given for the lowest four vibrational states over the full rotational bands and E2 transition probabilities are tabulated. Concluding remarks are presented in section 4.

## 2 Lagrange-mesh calculation of transition probabilities

### 2.1 Oscillator strength and transition probability per time unit

The dimensionless oscillator strength for an electric transition of multipolarity between an initial state and a final state is defined as [28, 29]

 f(λ)i→f=mecℏ(2λ+1)(λ+1)[(2λ+1)!!]2λk2λ−1S(λ)if2Ji+1 (1)

where is the electron mass,

 k=|Ef−Ei|ℏc (2)

is the photon wavenumber and

 S(λ)if=S(λ)fi=∑MiMfμ|%$⟨$γiJiMi|O(λ)μ|γfJfMf% ⟩|2=|⟨γiJi||O(λ)||γfJf⟩|2, (3)

where is a total angular momentum, is its projection, represents the other quantum numbers and the reduced matrix element is defined according to [30]. The transition irreducible tensor operator is given in units of by

 O(λ)μ=∑iZir′λiC(λ)μ(Ω′i) (4)

where is the relative coordinate of particle with respect to the center of mass, is its charge in units of and . Notice that the charge unit is included in the coefficient in . One has

 f(λ)f→i=2Ji+12Jf+1f(λ)i→f. (5)

If atomic units are used, the oscillator strength reads

 f(λ)i→f=(2λ+1)(λ+1)[(2λ+1)!!]2λα2λ−2(Ef−Ei)2λ−1S(λ)if2Ji+1 (6)

where is the fine-structure constant.

The transition probability per time unit for is given in atomic units (the atomic unit of time is s) by

 W(λ)i→f=2(λ+1)(2λ+1)λ[(2λ+1)!!]2α2λ+1(Ei−Ef)2λ+1S(λ)if2Ji+1 (7)

where all quantities are in a.u. For any multipolarity , the transition probability is related to the oscillator strength by

 W(λ)i→f=2α3(Ei−Ef)2f(λ)i→f. (8)

### 2.2 Dipole and quadrupole operators in perimetric coordinates

After elimination of the centre-of-mass motion, the Hamiltonian depends on the two Jacobi coordinates of proton 2 with respect to proton 1 and of the electron with respect to the centre of mass of both protons. These coordinates can be expressed as three Euler angles () defining the orientation of the triangle formed by the three particles, and three internal coordinates describing the shape of this triangle. The and angles correspond to the angular spherical coordinates of vector and the angle is the angular cylindrical coordinate of vector in the relative frame where the -axis is moved along by and rotations [31]. For the internal degrees of freedom we use the perimetric coordinates () defined as linear combinations of interparticle distances [27],

 x=R−re2+re1,y=R+re2−re1,z=−R+re2+re1, (9)

where and are the distances between the electron and protons 1 and 2, respectively. The domains of variation of these six variables are for and , for and for and .

In the body-fixed frame, the radial component of and the polar and axial components of are expressed in perimetric coordinates (9) as [15, 10]

 R = x+y2, (10) ρ = √xyz(x+y+z)(x+y)2, (11) ζ = (x−y)(2z+x+y)4(x+y). (12)

For H, the dipole tensor operator reads in Jacobi coordinates,

 d(1)μ=−(1+meM)r(1)μ (13)

where is the total mass of the molecular ion. This operator changes sign under space reflection (odd operator). It is invariant under the permutation of the protons. The tensor components of can be written as a function of the Euler angles as

 r(1)μ=ζD1μ0(ψ,θ,ϕ)+ρ√2[D1μ1(ψ,θ,ϕ)−D1μ−1(ψ,θ,ϕ)]. (14)

In both terms, the Wigner matrices change sign under space reflection [15] while and remain unchanged. With respect to proton exchange [15], and are both odd while and are both even.

In the Jacobi coordinate system, the E2 tensor operator reads

 Q(2)μ=√32{12[R(1)⊗R(1)](2)μ−γ[r(1)⊗r(1)](2)μ} (15)

where

 γ=1−2meM−m2eM2. (16)

In perimetric coordinates, it becomes

 Q(2)μ=12[R2−γ(2ζ2−ρ2)]D2μ0(ψ,θ,ϕ)−√32γζρ[D2μ1(ψ,θ,ϕ)−D2μ−1(ψ,θ,ϕ)] −√38γρ2[D2μ2(ψ,θ,ϕ)+D2μ−2(ψ,θ,ϕ)]. (17)

Operator is even with respect to parity and to permutation.

### 2.3 Transition matrix elements

The three-body Hamiltonian that we consider involves Coulomb forces between the particles but no spin-dependent forces. Hence the total orbital momentum and parity are good quantum numbers corresponding to constants of motion. The wave functions with orbital momentum and parity are expanded as [10]

 Ψ(Lπ)σM=L∑K=0DLπMK(ψ,θ,ϕ)Φ(Lπ)σK(x,y,z). (18)

In practice, the sum can be truncated at some value . The normalized angular functions are defined for by

 DLπMK(ψ,θ,ϕ) = √2L+14π(1+δK0)−1/2[DLMK(ψ,θ,ϕ) (19) +π(−1)L+KDLM−K(ψ,θ,ϕ)]

where represents a Wigner matrix element. They have parity and change as under permutation of the protons. Hence is symmetric for and antisymmetric for , when and are exchanged. Most bound states belong to the band where dominates and is equal to .

Since the symmetry of the proton spin part is where is the total spin of the protons, physical states (i.e. states antisymmetric with respect to the exchange of the protons) have . In the rovibrational band, states have natural parity, . Even- states have positive parity and . The protons are thus in a singlet state and the total intrinsic spin of the molecule is . Odd- states have negative parity and . The protons are in a triplet state and the total intrinsic spin of the molecule is or . E1 transitions are forbidden because of different proton symmetries when . However, they are possible from the three weakly bound states of the band where . E2 transitions are possible within the rovibrational band for , . Such states have the same parity and symmetry.

Using the property

 ⟨DL′M′K′|Dλμκ|DLMK⟩=8π22L′+1(LλMμ|L′M′)(LλKκ|L′K′), (20)

one can write for ,

 FLL′;λKK′;κ = ⟨DL′π′M′K′|Dλμκ+(−1)κDλμ−κ|DLπMK⟩ (21) = (2L+1)1/2[(2L′+1)(1+δK′0)(1+δK0)]1/2(LλMμ|L′M′) × {(LλKκ|L′K′)+(−1)κ(LλK−κ|L′K′) + π′(−1)L′+K′[(LλKκ|L′−K′)+(−1)κ(LλK−κ|L′−K′)]}.

The reduced matrix elements of the operators between initial and final states and can be written as

 ⟨Ψf(Lπff)σf||Q(2)||Ψi(Lπii)σi⟩=∑KiKf2∑κ=0(1+δκ0)−1FLiLf;2KiKf;κALiLfKiKf;κ (22)

where the perimetric matrix elements

 ALiLfKiKf;κ=⟨Φf(Lπff)σfKf|Aκ|Φi(Lπii)σiKi⟩ (23)

are calculated by integration over the perimetric coordinates with the volume element and, from (17),

 A0=12[R2−γ(2ζ2−ρ2)],  A1=−√32γζρ,  A2=−√38γρ2, (24)

with , and replaced by their expressions (10)-(12). With (3) and (22), the oscillator strength is given explicitly for transitions between natural-parity states by

 S(2)if=(2Li+1)∣∣∣∑KiKf{CLi2LfKi0KfALiLfKiKf;0 +[CLi2LfKi1Kf(1+δKi0)1/2−CLi2LfKi−1Kf(1+δKf0)1/2]ALiLfKiKf;1 +[CLi2LfKi2Kf(1+δKi0)1/2+CLi2LfKi−2Kf(1+δKf0)1/2−CLi2Lf2−11δKi1δKf1]ALiLfKiKf;2}∣∣∣2 (25)

where the Clebsch-Gordan coefficients are written as and the sums over and are truncated at in practice. The remaining calculation of the matrix elements is particularly simple within the Lagrange-mesh method as shown in the next subsection.

### 2.4 Lagrange-mesh method

The three-dimensional Lagrange functions are infinitely differentiable functions defined by

 FKijk(x,y,z)=N−1/2KijkRK(x,y,z)fNxi(x/hx)fNyj(y/hy)fNzk(z/hz). (26)

The one-dimensional Lagrange-Laguerre functions are given by

 fNi(u)=(−1)iu1/2iLN(u)u−uie−u/2 (27)

where is the Laguerre polynomial of degree and is one of its zeros, i.e. . They vanish at all with . Basis (26) is exactly equivalent to the set of functions with to . The mesh points correspond to the zeros , , of Laguerre polynomials of respective degrees , , . Three scale parameters are introduced in (26) in order to fit the different meshes to the size of the actual physical problem. The function is a regularization factor introduced because of the presence of singularities in the Hamiltonian operator when differs from zero [15, 10]. It is equal to 1 when and to otherwise. The normalization factor is defined as

 NKijk=hxhyhz(hxui+hyvj)(hxui+hzwk)(hyvj+hzwk)R2Kijk (28)

where .

The functions satisfy the Lagrange property with respect to the three-dimensional mesh , i.e. they vanish at all mesh points but one,

 FKijk(hxup,hyvq,hzwr)=(NKijkR−2KijkλNxiλNyjλNzk)−1/2δipδjqδkr. (29)

The coefficients , , are the Christoffel numbers which appear as weights in the Gauss-Laguerre quadrature approximation

 ∫∞0∫∞0∫∞0G(u,v,w)dudvdw≈Nx∑i=1Ny∑j=1Nz∑k=1λNxiλNyjλNzkG(ui,vj,wk). (30)

At the Gauss approximation scaled with , , , the three-dimensional Lagrange functions (26) are orthonormal with respect to the perimetric volume element because of the Lagrange property (29).

The functions of equation (18) are expanded in the Lagrange basis as

 Φ(Lπ)σK(x,y,z) = N∑i=1i−δK∑j=1Nz∑k=1C(Lπ)σKijk[2(1+δij)]−1/2 (31) ×[FKijk(x,y,z)+σπ(−1)KFKjik(x,y,z)]

where we use the same number of mesh points and the same scale factor for the two perimetric coordinates and in order to take advantage of the Lagrange conditions (29) when the two coordinates are exchanged. Because of the symmetrization the sum over is limited by the value , where is equal to 0 when and to 1 when .

The three-body Hamiltonian in perimetric coordinates for each good quantum number and its discretization on a Lagrange mesh are given in [15]. The singularities of the three Coulomb terms are automatically regularized in the matrix elements by the volume element so that the Lagrange-mesh method is not affected by those singularities. The potential matrix is diagonal and its elements are the values of the potential at the mesh points. The calculation would be as easy with other form factors for the potentials. The resulting matrix is rather sparse. The remaining problem is to calculate the lowest eigenvalues and the corresponding eigenvectors of a large sparse matrix.

For given , the eigenvalues in increasing order are labeled by a quantum number related to the vibrational excitation in the Born-Oppenheimer picture. The corresponding eigenvectors provide the coefficients appearing in expansion (31).

Let us consider initial and final components (31) with respective coefficients and . Because of the Lagrange property (29), the matrix elements (23) are simply obtained with the Gauss quadrature (30) as

 ALiLfKiKf;κ≈N∑i=1i−δK∑j=1Nz∑k=1Ci(Lπii)σiKiijkCf(Lπff)σfKfijkAκ(hui,hvj,hzwk) (32)

where .

## 3 Energies and E2 transition probabilities

### 3.1 H+2 bound and quasibound energies

The energies of the lowest vibrational bound states for to 35 have been calculated with the present method in [10]. Here we extend those calculations to quasibound states up to and to the first three excited vibrational states. The main reason making this extension possible is a better technique for searching the eigenvalues of a large sparse matrix [32] and faster personal computers.

Since the main aim is to calculate transition matrix elements involving two different wave functions, it is convenient to use a single three-dimensional mesh for all states. An excellent accuracy is obtained when the parameters of the calculation are chosen as , and , . For a given value, the total number of basis states is then 16400 or 15600 depending on the parity of . The size of the matrix is larger by about a factor ( when is limited by . For , like in [10], calculations are performed with . They correspond to a size of 48400. In order to make comparisons with more accurate literature results, we use for the proton mass the benchmark value a.u. The dissociation threshold is then at a.u. or Hartrees. The obtained energies are presented as the first line for each value in Table 1. The accuracy is estimated from the stability of the digits with respect to calculations with mesh points. The error is expected to be at most of a few units on the last displayed digit. Literature results sometimes truncated at the 15th digit are displayed in each second line. Except in the low- or low regions where other references are mentioned, the literature results are the 11-digit energies obtained by Moss [1].

Let us start the discussion with the ground state. The energy of the ground state has been improved with respect to Moss’ result in a number of papers [2, 3, 4, 5, 6, 7, 8] to reach an accuracy around or beyond 30 digits in [7, 8]. Our accuracy is about . For the first vibrational excited state, the accuracy is better than [7]. For the and states, the accuracies are about and , respectively [3]. There is thus about a one-digit loss for each vibrational excitation. When the numbers of mesh points are increased to and , the accuracies of all these states reach about (see the third line of the energies). The energy of the lowest state is known with about 25 digits [7, 8]. Our result has an accuracy better than . Comparisons with references [1, 3] indicate that the accuracies of the vibrational excited states behave like for . Results with an accuracy close to 18 digits are available for -12 and [11]. They show that our error for the lowest vibrational energy remains smaller than for all these states. When comparing the rest of our results with those of Moss [1], one observes that both works agree very well. We are a little more accurate for and a little less accurate for and 3. But in addition, the Lagrange-mesh method provides easy-to-use wave functions. The obtained spectrum is depicted in Fig. 1.

Typical convergence tests are displayed in Table 2 for two quite different sets of initial and final states. For the values and close to those suggested in [10], Table 2 displays initial and final energies and transition probabilities for various choices of and . The transition probability is obtained by restricting (25) to while corresponds to . First one observes that the contributions have an importance smaller than 0.05 % and that the contributions are smaller than 0.001 %. Second one sees that the convergence of the transition probabilities is slower than the convergence of the energies, as expected from the variational principle. A good convergence with respect to is already obtained for . The convergence with respect to is slower. Increasing is more expensive than increasing since the size of the basis increases with . Since the convergence is exponential, one can estimate that the relative accuracy on is about for and still better than for . This means that the wave functions are quite accurate. Further similar tests have been performed for other transitions.

The convergence of the transition probabilities with respect to can be studied by comparison with results from wave functions truncated at and . The relative error when is smaller than 0.3 % for all considered transitions while the error for is smaller than . This is rather similar to truncations with respect to . By extrapolation, we estimate that the relative error on the present transition probabilities obtained with should be smaller than .

The probabilities per second for transitions within a same rotational band, and , are presented in Table 3. They include some transition probabilities involving quasibound states. In accord with the previous discussion, we limit the number of significant figures to six. The probabilities increase slowly with with a maximum around for , for , for , and for , not far from the end of the rotational bands. This is due to a maximum of the energy differences around . The maximum of the transition probabilities is shifted toward higher values by a steady increase of the reduced matrix elements. The transition probabilities behave similarly for the displayed values with a slight decrease when increases.

The probabilities per second for other transitions are displayed in Table 4. The columns correspond to transitions between different vibrational states. For each value, the successive lines correspond to increasing values of , i.e., to for , for , and , respectively. For , the obtained probabilities agree with those of Posen et al[13] and improve them significantly. More than 95 % of the results of Posen et alexactly correspond to our result rounded at two significant figures. In most other cases, the rounding is of 6 units on the third digit rather than at most 5 for a normal rounding. An example of the few ’worst’ cases is the transition probability where our result in Table 4 is while the result of Posen et alis .

The strongest transition from each state occurs in general towards the nearest vibrational state () for . Exceptions can be found between and . For , in the vicinity of and beyond, the transitions are replaced by transitions because the sign of the energy difference changes (see the arrows in Fig. 1 for the transitions). These numbers are indicated in italics in Table 4. For example, the first number in the last line for corresponds to the transition. Hence the transition probabilities strongly vary in this region.

Oscillator strengths are depicted in Figs. 2, 3 and 4. For the transitions with displayed in Fig. 2, they present a strong variation along the band. They also vary strongly with . The transitions present a deep minimum around due to the change of sign of the energy difference (see Fig. 1). Beyond that value, the initial state is lower than the final state and the strengths are negative. Otherwise the strengths slowly increase with the vibrational excitation. The strengths are smaller by more than an order of magnitude. The minimum occurring around is here due to a change of sign of the matrix element appearing in (3). The strengths are smaller by more than an order of magnitude than the strengths. Except near the minimum at also due to a change of sign of the matrix element, they are rather constant along the band.

The oscillator strengths presented in Fig. 3 do not vary much along the bands as expected from the similar vibrational structures of the initial and final states. They slowly decrease for and slowly increase for . They are remarkably flat for . Except near the end of the band, for given , they increase with .

The strengths presented in Fig. 4 behave similarly to the strengths. Minima take place around for , for and for . These minima thus occur now at increasing values with increasing and are all due to a change of sign of the matrix element.