Molecule H2+

Quantum Mechanics
                       K. Konishi, G. Paffuti

Give a variational computation of the potential energy for (Subscript[H, 2]^(+)) molecule.

Global variables and initialization

ClearAll["Global`*"];
Off[General::spell, General::spell1];

Theory and method of solution

Coordinate system

In atomic units the Schrödinger equation for the ionized molecule H2 is

-(1/2) Δ Ψ  - ( (1/(Subscript[r, 1])) +  (1/(Subscript[r, 2])) )Ψ = E Ψ

where (Subscript[r, 1]) are (Subscript[r, 2]) are electron - nuclei distances.

Show[Graphics[{{RGBColor[1,0,0],Disk[{-1,0},0.05]},
{RGBColor[1,0,0],Disk[{1,0},0.05]},{RGBColor[0,1,0],Disk[{0.2,0.7},0.05]},
Line[{{-1,0},{1,0}}],Line[{{-1,0},{0.2,0.7}}],Line[{{1,0},{0.2,0.7}}],
Text[#[[1]],#[[2]]]&/@{{"R",{0,0.1}},{"(Subscript[r, 1])",{-0.5,0.4}},{"(Subscript[r, 2])",{0.65,0.4}},
{"H",{-1,0.17}},{"H",{1,0.17}}}},AspectRatio->Automatic,Axes->False]]

The equation is most easily solved by using elliptic coordinates. With R the inter-nuclear distance we define

ξ = (((Subscript[r, 1])+ (Subscript[r, 2]))/R); η = (((Subscript[r, 2])-(Subscript[r, 1]))/R);  1 ≤ ξ ≤ ∞;  -1 ≤ η ≤ 1;

In this coordinates the laplacian takes the form

Δ = (4/((R^2)((ξ^2)- (η^2))))((∂_ξ) (ξ^2-1)(∂_ξ) + (∂_η) (1-η^2)(∂_η)) + (1/(R^2)) (1/(((ξ^2)-1)(1-(η^2))))(Subscript[(∂^2), φ])

The system is invariant under rotation around z axis (directed along the two nuclei). For eigenstates of (Subscript[L, z]) we have Ψ = ψ Exp[i m φ]  and the Schrödinger is

-((∂_ξ) (ξ^2-1)(∂_ξ) + (∂_η) (1-η^2)(∂_η))ψ - 2 R ξ ψ +( ((m^2)/((ξ^2)-1))+ ((m^2)/(1-(η^2))))ψ = ((E (R^2))/2)((ξ^(2 ))-(η^2))ψ;

For simplicity we will consider only eigenstates with m=0. The equation is separable, but here we will study the problem from the variational view point. We will look for a solution written as a superposition of products of ξ and η functions.

Variational method

For more information on this method the reader is referred to the notebooks on the variational problem, We assume here some familiarity with this technique.

We will use as a set of basis functions Legendre polynomials in η and Laguerre polynomials (times an exponential) in x, where ξ = x+1. This choice is suggested by the limits in the variables written in ().

ψ = (∑_(i,j)) (Subscript[C, ij])  Exp[-x/2]LaguerreL[i,x] LegendreP[j,η]; 
0≤ i ≤ maxξ; 0≤ j ≤ maxη;  dim = (maxξ+1)* (maxη+1);

The parameter dim is the dimension of the base. Indexes i, j can be joined in a index I with

I = i + maxξ* j + 1 ;   ⇒  1 ≤ I ≤ dim.

As the basis vectors are factorized the matrix elements of ξ and η variables follow the same factorization. For instance, with obvious notation:

(i1 j1| (ξ^2) | i2 j2) = (i1|(ξ^2)|i2) * (j1| j2).

Once the matrices for the relevant operators are written a variational calculation can start. A possible form of the variational procedure is to introduce a scale parameter λ and look for the minimization in λ. As our computation as only a didactic purpose we will not introduce any additional scale and we will limit ourselves to a choice a large  basis.

The matrices

From matrix elements of () it appears that we need the following matrices

matD2x = ((∂_ξ) ψ[i], ((ξ^2)-1) (∂_ξ) ψ[j]); matx2 = ( ψ[i], (ξ^2) ψ[j]); matx1 = ( ψ[i], ξ ψ[j]);
matD2η = ((∂_η) ψ[i], (1-(η^2)) (∂_η) ψ[j]); matNη = ( ψ[i],  ψ[j]);      matη2 = ( ψ[i], (η^2) ψ[j]);

These have to be assembled together to build matrices defined in the tensor product space. This is done below by

(matKin - 2 R matV)ψ = E((R^2)/2 matN)ψ

Details on the matrices

Single variable matrix elements

Using recursion relations for polynomials and a bit of guessing the explicit expression for the  matrices can be given. We just quote the result, leaving to the reader the task of performing checks. In the following formulas indexes i, j run from 0. We write only matrix elements different from zero. The necessary symmetrization is understood, for instance (i1|mat|i2) = (i2|mat|i1).

matD2x

(i|matD2x|i+1) =  1/2 (i+1);              (i|matD2x|i+2) =   -1/2 -3/4 i - 1/4 i^2;
(i|matD2x|i) =  1 + 3/2 i + 1/2 i^2 .

matx2

(i|matx2|i+1) = -6 -10i-4 i^2;     (i|matx2|i+2) =   2 + 3 i + i^2;
 (i|matx2|i) =  5 +10 i + 6 i^2 ;

matx1

(i|matx1|i+1) = -1 -i;                        (i|matx1|i) = 2+2i;

matD2η

(i|matD2η|i) =  2 i (i+1)/(2i+1);

matNη

(i|matNη|i) =  2 /(2i+1);

matη2

(n|matη2|n+2) =  (2/((2n+5)))(-2+n+5 (n^2)+2 (n^3))/(-3-2 n+12 (n^2)+8 (n^3));
(n|matη2|n) =  (2/((2n+1)))(-1+6 (n^2)+4 (n^3))/(-3-2 n+12 (n^2)+8 (n^3)).

Code for matrices

getMatricesη[maxL_]:= Module[{matD2η,matNη,matη2},
matD2η = DiagonalMatrix[Table[2 i (i+1)/(2i+1),{i,0,maxL}]];
matNη = DiagonalMatrix[Table[2 /(2i+1),{i,0,maxL}]];
matη2 = Table[0,{maxL+1},{maxL+1}];
(* il +1 è dovuto al fatto che gli indici cominciano da 1 *)
Table[matη2[[n+1,n+2+1]] = (2/((2n+5)))(-2+n+5 (n^2)+2 (n^3))/(-3-2 n+12 (n^2)+8 (n^3));,{n,0,maxL-2}];
matη2 = matη2 + Transpose[matη2];
matη2 = matη2 +  DiagonalMatrix[Table[(2/((2n+1)))(-1+6 (n^2)+4 (n^3))/(-3-2 n+12 (n^2)+8 (n^3)),{n,0,maxL}]];
{matD2η,matNη,matη2}];

getMatricesξ[maxL_]:= Module[{matD2x,matx1,matx2},
matx1 = matx2=matD2x = Table[0,{maxL+1},{maxL+1}];
(* il +1 è dovuto al fatto che gli indici cominciano da 1 *)
Table[matD2x[[i+1,i+2]]= 1/2 (i+1);,{i,0,maxL-1}];
Table[matD2x[[i+1,i+3]]= -1/2 -3/4 i - 1/4 i^2;,{i,0,maxL-2}];
matD2x = matD2x + Transpose[matD2x];
matD2x = matD2x  + DiagonalMatrix[Table[1 + 3/2 i + 1/2 i^2,{i,0,maxL}]];

Table[matx2[[i+1,i+2]]= -6 -10i-4 i^2;,{i,0,maxL-1}];
Table[matx2[[i+1,i+3]]=2 + 3 i + i^2;,{i,0,maxL-2}];
matx2 = matx2 + Transpose[matx2];
matx2 = matx2  + DiagonalMatrix[Table[5 +10 i + 6 i^2,{i,0,maxL}]];

Table[matx1[[i+1,i+2]]= -1 -i;,{i,0,maxL-1}];
matx1 = matx1 + Transpose[matx1];
matx1 = matx1  + DiagonalMatrix[Table[2 +2 i,{i,0,maxL}]];

{matD2x,matx1,matx2}];

tensorMatrices[maxξ_,maxη_]:= Module[{matD2x,matx1,matx2,matD2η,matNη,matη2,
matKin,matV,matN,dim,maxindξ,maxindη,idξ},
maxindξ = (maxξ+1);
maxindη = (maxη+1);
dim = (maxξ+1)* (maxη+1);
matN = matKin = matV = 0*IdentityMatrix[dim];
idξ = IdentityMatrix[maxindξ];

{matD2x,matx1,matx2}=getMatricesξ[maxξ];
{matD2η,matNη,matη2}= getMatricesη[maxη];

Do[  Do[ Do[Do[
If[ Abs[i1-i2]≤ 2 && Abs[j1-j2]≤ 2,
i=i1+maxindξ*j1+1;
j=i2+maxindξ*j2 +1 ;
matN[[i,j]] = matNη[[j1+1,j2+1]]*matx2[[i1+1,i2+1]] - matη2[[j1+1,j2+1]]idξ[[i1+1,i2+1]];
matKin[[i,j]] = matD2x[[i1+1,i2+1]]*matNη[[j1+1,j2+1]] + 
matD2η[[j1+1,j2+1]]idξ[[i1+1,i2+1]];
matV[[i,j]] = matx1[[i1+1,i2+1]]matNη[[j1+1,j2+1]];];
,
{i1,0,maxξ}],{i2,0,maxξ}],{j1,0,maxη}],{j2,0,maxη}];
SparseArray/@{matKin,matV,matN}]

energy[R_?NumberQ]:= ((Sort[Eigenvalues[{matKin- 2 R matV,R^2/2*matN},8,Method->
{Arnoldi,Shift->-2.0000}]][[{1,2}]])+1/R )*2;

Computation

Example

We define the maximum number of terms in ξ and η variables and write the dimension of the resulting tensor space

maxξ = 15; maxη =15;
dim = (maxξ+1)* (maxη+1)

This define the global matrix variables needed in the Hamiltonian

{matKin,matV,matN}=tensorMatrices[maxξ,maxη];

This compute the energy for R=1.5 in Rydberg. The first two energy levels are given in output.

energy[1.5]

Stability with respect the basis

As a preliminary check let us verify that the result is stable under the increase of the base dimension. We take maxη=maxξ.

Table[{matKin,matV,matN}=tensorMatrices[maxξ,maxξ];
energy[1.25],{maxξ,5,20,5}]

Energy as a function of R

maxξ = 15;
{matKin,matV,matN}=tensorMatrices[maxξ,maxξ];

nuclearDistance = Union[Range[0.5,5,0.5*0.125],Range[1.75,2.25,0.0625]];

eneValues = energy/@ nuclearDistance;

Here is a plot as a function of R for the ground state and the first excited state. The dashed line at E = -1 Ry is the energy for infinite distance nuclei, the double ionized molecule.

ListPlot[{Transpose[{nuclearDistance,eneValues[[All,1]]}],
Transpose[{nuclearDistance,eneValues[[All,2]]}]},Joined-> True,
Epilog-> { Dashing[{0.015,0.015}],Line[{{0,-1},{5,-1}}]}]

In the text it has been shown that the R dependence of the ground state is given to the first approximation by

enAp[R_]:= ((2/R - ((2 E^-R (1+R)+2/R(1-(E^(-2R))(1+R)))/(1+(E^(-R))(1+R+R^2/3))))-1)

Let us compare the variational calculation with this approximate result. The dashed line is the approximate result:

ListPlot[{Transpose[{nuclearDistance,eneValues[[All,1]]}],
Transpose[{nuclearDistance,eneValues[[All,2]]}],
Transpose[{nuclearDistance,enAp[nuclearDistance]}]},
Joined-> True,
PlotStyle-> {{},{},Dashing[0.015]}]

The wave function

Let us fix an R value and compute eigenvalues and eigenvectors:

R=2.0;

maxξ = maxη= 7;
{matKin,matV,matN}=tensorMatrices[maxξ,maxη];

{eigH,eigV}=Eigensystem[{matKin- 2 R matV,R^2/2*matN}];
ordine = Ordering[eigH];
eigH = eigH[[ordine]][[Range[1,2]]];
eigV = eigV[[ordine]][[Range[1,2]]];

Let us call z the coordinate along the nuclear axes, y a coordinate in the perpendicular plane. Let us write the coordinate relations () on the plane (z,y) in cartesian coordinates:

αCart = (1/R)((Sqrt[(((z+(R/2)))^2)+(y^2)]) + (Sqrt[(((z-(R/2)))^2)+(y^2)]))-1;
ηCart = (1/R)(-(Sqrt[(((z+(R/2)))^2)+(y^2)]) + (Sqrt[(((z-(R/2)))^2)+(y^2)]));

The basis vectors

The basis vectors are given by the tensor product of the ξ basis and of the η basis. In Mathematica the tensor multiplication is Outer, (see Help), then our basis can be defined by

fBase[α_,η_]:= Flatten[Outer[Times, LegendreP[Range[0,maxη],η],
 Exp[-α/2]LaguerreL[Range[0,maxξ],α]]];

If φ[[i]] are the basis functions and C[[i]] the eigenvectors computed with Eigensystem the Schrödinger wave function is  (∑_i)C[[i]] φ[[i]], then

ψ[state_,α_,η_]:= eigV[[state]].fBase[α,η];

Plots

Let us compute the wave function on a list of z values to do a plot

zVar = Range[-3,3,0.2];

psiVal1 =ψ[1,(1/R)((Sqrt[(((#+(R/2)))^2)]) + (Sqrt[(((#-(R/2)))^2)]))-1, (1/R)(-(Sqrt[(((#+(R/2)))^2)]) + (Sqrt[(((#-(R/2)))^2)]))]&/@zVar;
psiVal2 =ψ[2,(1/R)((Sqrt[(((#+(R/2)))^2)]) + (Sqrt[(((#-(R/2)))^2)]))-1, (1/R)(-(Sqrt[(((#+(R/2)))^2)]) + (Sqrt[(((#-(R/2)))^2)]))]&/@zVar;

The two lowest states are shown below, the red  dashed line representing the excited state

ListPlot[{Transpose[{zVar,-psiVal1}],Transpose[{zVar,-psiVal2}]},Joined-> True,PlotStyle-> {Black,{Red,Dashing[0.015]}}]

It is apparent the characteristic double-well type shape for the symmetric ground state and the antisymmetric excited state.

A contour plot for  |ψ| can be obtained by expressing ξ-1 and η in cartesian coordinates:

makeContour[{stato_,R_},{zmin_,zmax_,ymax_}]:=
ContourPlot[Abs[ψ[stato,(1/R)((Sqrt[(((z+(R/2)))^2)+(y^2)]) + (Sqrt[(((z-(R/2)))^2)+(y^2)]))-1, (1/R)(-(Sqrt[(((z+(R/2)))^2)+(y^2)]) + (Sqrt[(((z-(R/2)))^2)+(y^2)]))]],{z,zmin,zmax},{y,-ymax,ymax},
ContourShading-> False,ContourLines->True,AspectRatio-> Automatic]

Show[GraphicsGrid[{{makeContour[{1,2.0},{-2,2,1.5}],makeContour[{2,2.0},{-2,2,1.5}]}}],ImageSize-> {8 72,3.5 72}]