Two dimensional problems

Quantum Mechanics
 K. Konishi, G. Paffuti

Study the variational procedure for the two-dimensional problem:

V[x,y] = ((x^2)/2)+((y^2)/2)+(x^2)(y^2).

Apply the method to a regularization of the Henon potential

V[x,y]= α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1(x^3-3 x (y^2))= α1 (r^4) + β1 (r^2) + γ1 (r^3) Cos[3θ].

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

Introduction

General remarks

In this notebook we will present some elementary formulations of the variational principle for two-dimensional Schrödinger equation. In particular we shall cover the following topics

The main functions used in the notebook will be collected in a subsection below named functions.

The notebook consists of three main sections, one for each of the above methods. Each section is endowed with a subsection in which the numerical techniques are presented and the functions briefly commented and  a subsection devoted to the checks. These two parts can be skipped by the reader not interested in the details of Mathematica functions, possibly with the exception of the first section, as for the factorized Hamiltonian one cannot write a general purpose function. The last part of each section is devoted to an application. We will consider a kind of deformed anharmonic oscillator, which is of some importance. In a certain limit it describes one of the classically well known example of chaotic system and for our purposes it presents some interesting symmetry properties.

Before starting let us make some comments on the above points.

1 Spectral method

Let us call x and y the two coordinates describing our problem. The Hilbert space for the two dimensional problem is  (𝕃^2)[ℝ]⊗ (𝕃^2)[ℝ]i.e. the space of normalizable functions of two variables (we assume cartesian coordinates for definiteness and no boundaries).  Let (Subscript[f, i])[x]a complete system in the one dimensional Hilbert space (𝕃^2)[ℝ], every normalizable state can be written in the form

ψ[x,y] = (∑_(i,j))(Subscript[c, ij]) (Subscript[f, i])[x] (Subscript[f, j])[y]

The simplest variational technique, the spectral technique, is to consider a finite approximation of the infinite sum (). In principle this is exactly the same procedure followed in the one dimensional case, the only complication being a double index instead of a single index to identify basis functions. One can enumerate all the pairs of indexes appearing in (), let us call I the index which identify a couple. The sum () can be re-written as

ψ[x,y] = (∑_I) (Subscript[C, I]) (Subscript[F, I])[x,y];    (Subscript[F, I])[x,y] = (Subscript[f, i])[x] (Subscript[f, j])[y]

and this completes the analogy with one dimensional case. We will obtain a matrix formulation of  the Schrödinger equation in the form

(Subscript[T, IJ]) (Subscript[C, j]) + (Subscript[V, IJ])(Subscript[C, J]) = E (Subscript[C, J])

where T denotes the kinetic term and V the potential term.

A particular form of Hamiltonian terms is a factorized form:

(Subscript[A, x]) ⊗ (Subscript[B, y])

This is for instance the case for the kinetic term in cartesian coordinates, where the usual expression

T = (1/(2m))(Subscript[p, x]^2) + (1/(2m))(Subscript[p, y]^2)

has indeed to be correctly written as

T = (1/(2m))(Subscript[p, x]^2) ⊗ 𝟙+ 𝟙⊗(1/(2m))(Subscript[p, y]^2)

From now on we will use the natural units and omit all factors m, ℏ etc.

A potential term with the same kind of structure can be

(Subscript[V, 2])=(1/2)(x^2) ⊗ 𝟙+ 𝟙⊗(1/2)(y^2)

An Hamiltonian as T+ (Subscript[V, 2]) is quite trivial as it is the sum of two completely independent Hamiltonians, and the spectrum is simply

E = (Subscript[E, 1])+ (Subscript[E, 2])

where (Subscript[E, 1])and (Subscript[E, 2])are eigenvalues of the separate systems.

The next, a slightly more complicated case is a potential term of the form, for instance

(Subscript[V, b])=(x^2)(y^2)

which is of the form (). If the Hamiltonian can be written as a sum of terms of this kind the computation of matrix elements in () greatly simplifies as they can be written as products of one-dimensional matrix elements:

∫𝕕x 𝕕y (Subscript[f, (Subscript[i, 2])])[x] (Subscript[f, j2])[y] A[x]B[y]  (Subscript[f, (Subscript[i, 1])])[x] (Subscript[f, j1])[y]= ∫𝕕x  (Subscript[f, (Subscript[i, 2])])[x]  A[x]  (Subscript[f, (Subscript[i, 1])])[x]∫𝕕y (Subscript[f, j2])[y]B[y](Subscript[f, j1])[y]

This will be our first kind of generalization.

Let us note that the sometimes the classification of factorized terms depends on the coordinates: for instance 1/r is factorized in polar coordinates (r,θ) while it is not in cartesian coordinates.

The only problem seems to be a possible headache caused by the identification of a pair of indexes (i,j) with a single index I. This problem is in fact is one of the primary source of errors in writing computer codes,  but it is not the worse problem, as we will see.

In particular let us note that in one-dimensional problems a typical dimension for the basis system was one hundred. In two dimensions this means 100 for each dimension, i.e. (10^4)  base functions and matrices of the order of (10^4)×(10^4). This can be  a serious problem, both for time and for storage requirements.

2 Approximate spectral method

The obvious generalization is the study of potentials of the general form V[x,y]. In general,  the integrals appearing in matrix element computations are not analytically known and one has to resort to numerical integration. For quite obvious time constraints we will limit ourselves to Gaussian integration techniques. In the cases shown below the execution time for solving a Schrödinger equation with this method will be two order of magnitude greater than the previous method, but reasonable fast anyway. The good point is that the execution time at fixed dimension of the basis is virtually of the same order of magnitude for every potential. The bad feature is that in numerical integration some precision will be lost, and this can be a limitation if subtle points have to be studied, as degeneracies for higher states, details of the eigenfunctions etc.

3 Pseudospectral method

Finally we will give a simple version of the pseudospectral method. This method is virtually equivalent to the previous one. It can have some advantages in some cases as fewer values of the potential have to be computed.

Functions

Hermite functions

hermiteMatrices[dim_,parity_Integer:-1]:=
Module[{x2,p2,x4,x1}, 
Block[{i,n,shift},
If[parity == 0 || parity == 1,n =2i-2+parity; shift = 1;,n=i-1;shift =2;];x2=SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,k_}/;k-i==shift->(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}];x2 =(x2+Transpose[x2]);
p2 =SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,k_}/;k-i==shift->-(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}];
p2 = (p2+Transpose[p2]);
x4=SparseArray[{{i_,i_}->(3/8)(2(n^2)+2n+1),{i_,k_}/;k-i==shift->((4n+6)/4)(Sqrt[(n+1)(n+2)]),
{i_,k_}/;k-i==2shift->(1/4)(Sqrt[(n+1)(n+2)(n+3)(n+4)])},{dim,dim}];
x4 =  (x4+Transpose[x4]);
If[parity == 0 || parity == 1,Return[{p2,x2,x4}];];
x1=SparseArray[{{i_,k_}/;k-i==1->(Sqrt[(n+1)/2])},{dim,dim}];
x1 = (x1+Transpose[x1]);
{p2,x2,x4,x1}]];

getscale[n_,g_,μ_]:= x/.FindRoot[(2n+1)(x-(μ/(x^3)))-((3g)/(x^5))(2(n^2)+2n+1),{x,(((3g))^(1/3))}]

recursiveHermitePolynomials= Compile[{{n,_Integer},{x,_Real,1}},
Module[{y1,y0,p1,ris,ris0},
y0 = 1.0 + 0.*x;
y1 = 2.x;
ris0={y0,y1};
ris=Table[
p1=2. x y1 - 2.  k y0;
y0=y1;
y1=p1;
p1,{k,1,n-2}];
ris=Join[ris0,ris];
ris]];

fBase[k_,x_]:= (1/(π^(1/4))) (1/(Sqrt[(2^(k-1))(k-1)!]))HermiteH[k-1,x]Exp[-x^2/2]

For fixed parity we can use the decomposition () to write

fBase[n_,x_,par_]:=(1/(π^(1/4)))(1/(Sqrt[(2^(2n-2+par))(2n-2+par)!]))HermiteH[2n-2+par,x]Exp[-x^2/2]

This function computes normalized Hermite polynomials, in the form needed for oscillator wave functions. We write two versions of it , one with arguments [dim, r] gives all polynomials up to degree dim-1.  The second with three arguments, [dim, r, par], gives the first dim polynomials of parity fixed by par, (0 for even polynomials, 1 for odd polynomials).

hermiteP[dim_,xIn_]:= Module[{cf,x,ris,coeff},
x = Flatten[{xIn}];
 coeff=Table[(1./((N[π])^(1/4))) (1/(Sqrt[(2.^(k-1))(k-1)!])),{k,1,dim}]; 
ris= recursiveHermitePolynomials[dim,x]*coeff;
If[Length[xIn]== 0,Flatten[ris],ris]];

hermiteP[dim_,xIn_,par_]:= Module[{cf,dimUp,index,x,ris},
If[par== 0 || par == 1,
x = Flatten[{xIn}];
dimUp=2 dim - 2 + par+1;
index=Range[par+1,2 dim-2+par+1,2];
ris=hermiteP[dimUp,x][[index]];
ris,
hermiteP[dim,xIn]]];

A more efficient version would use recursion relation only for even or odd polynomials: we leave it an exercise to write such a routine.

Finally getPsi gives the wave function for a coefficient list eigv at points r. Parity can be specified, if not all polynomials are considered. This is the default.

getPsi[eigv_,r_,par_: -1]:=
If[par== 1|| par==0, 
(eigv.hermiteP[Length[eigv],r,par])*Exp[-(r^2)/2],
(eigv.hermiteP[Length[eigv],r])*Exp[-(r^2)/2]];

eigfHO[k_,x_,ω_]:=(ω^(1/4)) (1/(π^(1/4))) (1/(Sqrt[(2^(k-1))(k-1)!]))HermiteH[k-1,(ω^(1/2))  x]Exp[-ω x^2/2]

hermiteRoots[n_]:= Module[{J,δ,γ,radEig,w},
  γ = Sqrt[Range[1,n-1]/2];J = Normal[SparseArray[{Band[{2,1}]-> γ,Band[{1,2}]-> γ}]];
radEig =Sort[ Eigenvalues[N@J]];
w=(((2^(n-1))(Sqrt[π]))/(n^2))((n!)/((HermiteH[n-1,radEig])^2));
{radEig,w}]

Functions for two dimensional systems

makeData[α_,x_,y_,v_,dim1_,dim2_,listpar_:{-1,-1}]:= Module[
{matf,v1,v2,values,pgrid,p0,arg,parx,pary},
matf = Partition[v,dim2];
{parx,pary} = listpar;
If[parx==0 || pary == 1,
v1=Transpose[hermiteP[dim1,α x,parx]]*Exp[- (α x)^2/2];,
v1=Transpose[hermiteP[dim1,α x]]*Exp[- (α x)^2/2];];
If[pary== 0|| pary== 1,
v2=Transpose[hermiteP[dim2,α y,pary]]*Exp[- (α y)^2/2];,
v2=Transpose[hermiteP[dim2,α y]]*Exp[- (α y)^2/2];];
(* a factor α is inserted to normalize eigenfunctions at this scale *)
values= α v1.matf.Transpose[v2];
arg=Arg[Select[Flatten[values],Abs[#1]>1./dim1&,1][[1]]];
values = Chop[values*Exp[-I arg]];
pgrid =Flatten[ Outer[List,x,y],1];
Partition[Flatten[Transpose[{pgrid,Flatten[values]}]],3]]

funzIndice[ω_]:=Module[{ind,dimIndR,fromInd},
ind=0;
fromInd = {};
Do[
Do[
If[i1+i2 ≤ ω ,
ind=ind+1;
indR[i1,i2]= ind;
fromInd=Append[fromInd,{i1,i2}];
];,
{i2,1,ω-i1}],{i1,1,ω}];
dimIndR = ind;
{dimIndR,indR,fromInd}];

getMatV[dim_,{parx_,pary_},{dimGauss1_,dimGauss2_},potV_,α_]:= Module[
{GaussGridx,weightsx,GaussGridy,weightsy,ff,xvalues,yvalues,i2,j2,
hPx,hPxW,hPy,hPyW,VontheGrid,matV,dimIndR,indR,fromInd,matV1,inda,indb},
(* compute grid and weights for Gaussian integration *)
{dimIndR,indR,fromInd}=funzIndice[dim];
{GaussGridx,weightsx} = hermiteRoots[dimGauss1];
{GaussGridy,weightsy} = hermiteRoots[dimGauss2];
(* normalized Hermite polynomials without and with weights *)
(* The format is  h[i,x] *)
hPx=hermiteP[dim,GaussGridx,parx];
hPxW =Transpose[Transpose[hPx]*weightsx];
hPy=hermiteP[dim,GaussGridy,pary];
hPyW =Transpose[Transpose[hPy]*weightsy];
(* compute potential on Gauss grid points, the format is V[x,y] *)
VontheGrid= Table[potV[GaussGridx[[i]]/α , GaussGridy/α],{i,1,Length[GaussGridx]}];
(* compute off diagonal matrix elements of the potential V[k,j] *)
xvalues = Flatten[Outer[Times,hPx,hPxW,1],1];
yvalues = Flatten[Outer[Times,hPy,hPyW,1],1];
matV1 = xvalues.VontheGrid.Transpose[yvalues];
matV = ConstantArray[0,{dimIndR,dimIndR}];
(* kx = fromInd[[ind1,1]],  jx= fromInd[[ind2,1]] etc *)
ff[{i_,j_}]:= matV1[[i,j]];
inda=Table[fromInd[[ind1]],{ind1,1,dimIndR}];
matV=Table[
{i2,j2} = fromInd[[ind2]];
indb=Transpose[{inda[[All,1]]+dim (i2-1),inda[[All,2]]+dim (j2-1)}];
ff/@indb,
{ind2,1,dimIndR}];
matV]

getp2[dim_,{parx_,pary_}]:= Module[{matp2,p2,x2,x4,x1,op,dimIndR,
indR,fromInd,unitM,p2x,p2y,index,matp20},
{dimIndR,indR,fromInd}=funzIndice[dim];
If[parx == 0|| parx== 1,
{p2x,x2,x4}=hermiteMatrices[dim,parx];,
{p2x,x2,x4,x1}=hermiteMatrices[dim];];
If[pary== 0|| pary== 1,
{p2y,x2,x4}=hermiteMatrices[dim,pary];,
{p2y,x2,x4,x1}=hermiteMatrices[dim];];
unitM = SparseArray[IdentityMatrix[dim]];
(* the role of px and py must be interchanged to use the index below *)
matp20 = KroneckerProduct[p2y,unitM]+KroneckerProduct[unitM,p2x];
(* kx = fromInd[[ind1,1]],  jx= fromInd[[ind2,1]] etc *)
index = Table[{i1,j1} = fromInd[[ind1]];i1+dim (j1-1),{ind1,1,dimIndR}];
matp2= matp20[[index,index]];
Print["Dimension of the basis: "<>ToString[dimIndR]];
matp2];

Functions for pseudospectral method

getDCheb[x0_,L_,dim_]:= Module[{D1,D2,grid,fdf},
grid= x0 +(1/2) L(1.0 -  Cos[π Range[0,dim+1]/(dim+1)]);
fdf = NDSolve`FiniteDifferenceDerivative[2, grid, DifferenceOrder->"Pseudospectral"];
     D2=fdf[DifferentiationMatrix[]];
{grid[[2;;dim+1]],D2[[2;;dim+1,2;;dim+1]]}]

getD2[dim_,rad_:0]:= Module[{matM2,matK,ξgrid,D2,D2X,xgridL1},
{ξgrid,D2}= getDCheb[-1,2,dim];
matM2 = DiagonalMatrix[(((1-((ξgrid)^2)))^3)];
matK = DiagonalMatrix[(3/4)(2- ((ξgrid)^2))(1-((ξgrid)^2))];
D2X =SparseArray[ matM2.D2 + matK];
xgridL1= (( ξgrid)/(Sqrt[1-((ξgrid)^2)]));
{ξgrid,xgridL1,D2X}];

Spectral method

General remarks

The tensor products in Mathematica

A first good news is that Mathematica has a built-in function which makes tensor products, it is named KroneckerProduct. It makes a product in all possible combinations of two lists, and arrange the results in a matrix form.

Here is how it works. Let us take two lists (vectors) a and b of length 2, their tensor product is

vecA =Array[(Subscript[a, ##])&,2] ;
vecB =Array[(Subscript[b, ##])&,2] ;

(coeffMatrixForm=KroneckerProduct[vecA,vecB])//MatrixForm

You can imagine that a and b are coefficient of a wave function expansion in some orthogonal basis. The final space has dimensions 2×2 = 4 and the coefficients in this bigger spaces can be represented in the above form or flattening the matrix. This is equivalent to express coefficients in the form (Subscript[c, ij])or (Subscript[c, I])as explained in the first section. Here I = 1,...4.

coeffVectorForm = Flatten[coeffMatrixForm]

Let us now consider operators, i.e. matrices. Let mA and mB two matrices acting on first and second space respectively, for instance they can be the results of x and y integrals respectively in the Schrödinger picture. Their tensor product mA⊗mB is simply

matA =Array[(Subscript[mA, ##])&,{2,2}] ;
matB =Array[(Subscript[mB, ##])&,{2,2}] ;

(Op2Dim = KroneckerProduct[matA,matB])//MatrixForm

It a 4×4 matrix, i.e. is written in the language of "big indexes"  (Subscript[M, IJ]), I, J = 1, ... 4.

To convince ourselves that everything is correct let us see the effect on a vector:

Simplify[Op2Dim.coeffVectorForm]

This is exactly what one expects  for factorized operators: matA acts only on vector a and matB only on vector b. As a check we can first act with matrices and after take the KroneckerProduct:

Flatten[KroneckerProduct[matA.vecA,matB.vecB]]

A simple example

Let us consider a two dimensional harmonic oscillator of the form

H = (1/2)(Subscript[P, x]^2) + (1/2)(x^2) + (1/2)(Subscript[P, y]^2) + ((ω^2)/2)(y^2)

We take ω ≠ 1 to distinguish clearly the two coordinates. The spectrum is given by all sums of the form (1/2 + n1) + (1/2+n2) ω. These can be generated by Table

ω = 0.7;

eneH0 = Sort[Flatten[Table[((1/2)+n1) + ((1/2)+n2)ω,{n1,0,10},{n2,0,10}]]][[1;;10]]

Or better with the function Outer (see Mathematica Help , KroneckerProduct is a particular case of Outer):

(Sort@Flatten@Outer[Plus,(1/2)+Range[0,10],ω((1/2)+Range[0,10])])[[1;;10]]

We know how to construct matrices for one dimensional oscillators.  We can use the function hermiteMatrices which essentially gives us matrix elements of an harmonic oscillator of unitary frequency. Let us take an approximation with a basis of 20 elements (for dimension)

dim = 20;

Here are the harmonic oscillator matrices

{p2,x2,x4,x1}=N/@hermiteMatrices[dim];

The Hamiltonian is trivially written. We need the identity matrix:

unitM = IdentityMatrix[dim];

To simplify the notations we introduce  (To obtain ⊗ from the keyboard you have to type  Esc-c*, or you can write explicitly KroneckerProduct in the formulas below) :

a_⊗b_ := KroneckerProduct[a,b]

In the usual variational calculation a scale parameter α is introduced, original coordinates X are expressed as x/α. We choose arbitrarily α = 0.8, to check normalization and consistency of the formulas. In the rescaled coordinates:

α = 0.8;

ham = ((α^2)/2)p2⊗unitM + ((α^2)/2)unitM⊗p2 +(1/(2(α^2)))x2⊗unitM+(ω^2)(1/(2(α^2)))unitM⊗x2;

Let us compute its eigenvalues and compare with the analytical result:

Sort[Eigenvalues[Normal[ham]]][[1;;10]]

These are in rough agreement with above results. Let us note that in this trivial example we have used 400×400 matrices.

The wave functions

Using Eigensystem[H] we will obtain eigenvalues and eigenvectors. Matrices tend to be quite large and it is more efficient to use the option Method → Arnoldi in Eigensystem (and also in Eigenvectors). This method of solution is of recursive nature and allows the computation of a subset of all eigenvectors. To improve its accuracy for high excited levels it is better to compute some more eigenvectors than those needed. An important point is the choice of a shift parameter: we will be interested on a set of levels starting from the ground states, in this case the shift can be settled to a number lower than energy of the ground state, for instance the minimum of the potential. Another possibility is to add a constant to the potential in such a way that all eigenvalues are positive. Eigenvalues are ordered in decreasing order of absolute value. With this choice specifying a parameter -n in Eigensystem we will have the first n eigenvalues. For harmonic oscillator there are nothing to do as the potential is positive definite. Here are the first eigenvectors and eigenvalues

dim = 20;
numEig = 30;

{eigH,eigV} = Eigensystem[ham,- Min[numEig+50,Floor[((Length[ham])^2)/2]],
Method-> Arnoldi];
order = Ordering[eigH];
eigH=eigH[[order]];
eigV = eigV[[order]];

We have computed numEig + 50 = 80 eigenvalues. A given eigenvector, let us say the first one to fix our ideas, is given by a list of coefficients (Subscript[c, I]) I= 1,... (dim^2), or in general I = 1, ... dim1×dim2 if two different dimensions have been chosen for x and y coordinates. To write a wave function we must do the inverse operation of Flatten, i.e. transform the list in a matrix: this is done using Partition.  For example

vecA = Array[(Subscript[a, ##])&,6]

Partition[vecA,3]

With Partition[eigV[[1]], dim2] we write a matrix and transform our coefficient (Subscript[c, I ])in the form (Subscript[c, ij]). With this reshape we can write

ψ[x,y] = ∑ (Subscript[f, i])[x] (Subscript[c, ij])(Subscript[f, j])[y]

This operation can be very slow, and essentially useless for doing plots. It is simpler to define a grid of points, xg ⊗ yg and compute the function at once on the whole grid, using listability properties of built in Mathematica functions. As for the base functions we use the version introduced in notebook [], i.e. we compute recursively Hermite polynomials. The reader can use built-in function HermiteH, and this function must be used if the results to arbitrary precision are requested. With the function hermiteP we can compute the whole set of Hermite polynomials at once on the grid, substitution in () gives wave function. All this work is done in the function makeData, which takes into account some other points.

The function must be called as

data = makeData[α, xg, yg, vec, dim1, dim2, {parx, pary}]

The output is a list of triplet in the form {x,y,f}, adapted to be used in ListPlot3D. The last argument, the list of parities, is optional.

Here is an example. We used α = 0.8 and all states in the above computation, then:

xg = Range[-5,5,0.2]; 
yg = Range[-5,5,0.2];

data = makeData[α,xg,yg,eigV[[1]],dim,dim,{-1,-1}];

The utilities absData and probData transform the list in absolute values and square values for ψ.

absData[data_]:= (#*{1,1,Sign[Last[#]]})&/@data;
probData[data_]:= (#*{1,1,Last[#]})&/@data;

Here is a plot of the wave function

ListPlot3D[data,PlotRange->All,Mesh->None,Axes->False,Boxed->False]

and here is a bunch of states

data = makeData[1,xg,yg,eigV[[#1]],dim,dim,{-1,-1}]&/@{3,6,12};
GraphicsRow[ListPlot3D[absData[#1],PlotRange->All,Mesh->None,AxesLabel->{"x","y",""},Boxed->False]&/@data,ImageSize->8 72]

We note that as we have fixed ω = 0.7 for the y-oscillator and ⟨(y^2)⟩∼ 1/ω, the bell-like curves are wider in y direction.

Check on the absolute normalization

To check the normalization of our results let us consider the exact eigenstates of harmonic oscillator

eigfHO[k_,x_,ω_]:=(ω^(1/4)) (1/(π^(1/4))) (1./(Sqrt[(2^(k-1))(k-1)!]))HermiteH[k-1,(ω^(1/2))  x]Exp[-ω x^2/2]

Let us consider the ground state of the system

data = makeData[α,xg,yg,eigV[[1]],dim,dim,{-1,-1}];

and take a sample of values from our vector data

offset = Length[xg]*(Floor[Length[xg]/2]-1);
sample=data[[offset+1;;offset+10]]

In the sample the x coordinate is fixed to - 2 while y range from - 5 to - 3.2 . The theoretical expected values are

expectedValues=Table[eigfHO[1,-0.2,1]eigfHO[1,y,ω],{y,Range[-5,-3.2,0.2]}]

to be compared to

sample[[All,3]]

Max[Abs[sample[[All,3]]-expectedValues]]

We stress that the good agreement obtained with this basis depends on the value of α, here the value 0.8 has been chosen after a look at the stability region for eigenvalues, as usual in variational calculations. This point has been omitted as it is well known to the reader for a harmonic oscillator. The correlation between eigenvalues and scale factor will be explicitly shown in the more complicated model below.

Study of a model case

The model

The model we want to take for example is described by the potential

V[x,y]= α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1(x^3-3 x (y^2))= α1 (r^4) + β1 (r^2) + γ1 (r^3) Cos[3θ]

r, θ are polar coordinates in the plane. For α1 = 0 this is known as Hénon-Heiles potential. It is an example where a chaotic motion occurs in Classical Mechanics. Its quantum regularized version has been studied by several groups, in particular we mention the work of M.Trott [1], where a beautiful presentation of the model is given using Mathematica . Several references and a brief presentation of some classical aspects of the model are given there. The study is performed in radial coordinates. For reasons that will be explained at the end of this notebook the radial coordinates are in general more difficult to handle in two dimensional problems, even if in this case matrix elements can be computed analytically.

We are interested in this model for its symmetry. As it is apparent from () the potential is invariant under rotations ℛ of 2π/3 in the plane, and these form an abelian group with elements (𝟙, ℛ, (ℛ^2)) , and it is also invariant under the reflection 𝒫: y→ - y. The six elements (𝟙, ℛ, (ℛ^2), 𝒫, 𝒫 ℛ, 𝒫 (ℛ^2)) form a non abelian group, as σ and ℛ do not commute. An explicit representation of the operators can be given in the two dimensional space where (x,y) are the coordinates of a point. We have

ℛ = (({{Cos[2π/3],  - Sin[2π/3]}, {Sin[2π/3], Cos[2π/3]}}))= (({{-(1/2), -((Sqrt[3])/2)}, {((Sqrt[3])/2), -(1/2)}})) =- (1/2)(Subscript[σ, z]) + ((Sqrt[3])/2) (Subscript[σ, y]);
𝒫 = (({{1, 0}, {0, -1}})) = (Subscript[σ, z]);     [ ℛ, 𝒫] ≠ 0.

The presence of two noncommuting symmetries implies a degeneracy in the Hamiltonian. A complete set of observables can be taken as the Hamiltonian H and the "parity" 𝒫.  A mode detailed study of this model has been presented in notebook [2.1]. Here we simply write a prototype Hamiltonian with this kind of invariance

H = (({{1, -δ, -δ}, {-δ, 1, -δ}, {-δ, -δ, 1}}))

which describe three almost degenerate states, δ being a tunneling amplitude. That this is the case is shown by plotting the potential ()

α1=1;β1=- 1/2;γ1 = -3;

minV = First[FindMinimum[α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1(x^3-3 x (y^2)),{x,y}]];

Plot3D[α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1(x^3-3 x (y^2))-minV,{x,-3,3},{y,-3,3},PlotRange-> {{-3,3},{-3,3},{0,15}}]

We see three separate minima, each one is semiclassically a stable configuration. We expect an almost threefold degeneracy broken by tunneling effects. The degeneracy is not completely removed, as it happens in the analogous one dimensional case: states non invariant under ℛ are degenerate. We can classify the eigenstates with |E , 𝒫⟩.

Preliminary study : scale dependence

Let us take the parameters of the potential fixed at above given values

α1=1;β1=- 1/2;γ1 = -3;
minV = First[FindMinimum[α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1(x^3-3 x (y^2)),{x,y}]];

First we build one dimensional matrices. We need also the operator x^3 which is not directly given by hermiteMatrices. As x - operators connects the states which differ at most by 1 in principal quantum number, the easiest way to write x^3 is x^2 . x where matrices are written in dimensions dim+3 and finally cut out as dim×dim matrices: this is an exact procedure for polynomial interactions in Hermite harmonic oscillator basis.

dim = 20;
{p2,x2,x4,x1}=hermiteMatrices[dim+4];
x3 = x2.x1;
{p2,x2,x4,x1,x3}=Take[#1,{1,dim},{1,dim}]&/@N/@{p2,x2,x4,x1,x3};

Here we write the Hamiltonian using the function KroneckerProduct. First we write single operators. We shift the potential by its minimum value:

𝒰1 =SparseArray@IdentityMatrix[dim];
kin = (1/2) p2⊗𝒰1 +(1/2) 𝒰1⊗p2;
V0 = -minV 𝒰1⊗𝒰1;  
V2 = β1 (x2⊗𝒰1 + 𝒰1⊗x2);
V4 = α1 (x4⊗𝒰1 + 2 x2⊗x2+𝒰1⊗x4);
V3 = γ1 (x3⊗𝒰1 - 3 x1⊗ x2 );

We use a scale parameter α for optimization :

hamHH[α_]:= (α^2) kin + V0 + (1/(α^2))V2 + (1/(α^3))V3 + (1/(α^4))V4;

Here is a simple calculation of eigenvalues:

Re[ Sort[Eigenvalues[hamHH[1.5],-dim,Method-> Arnoldi]]][[1;;Min[25,dim]]]

We see that in fact states appear in "triplets" almost degenerate. A more precise computation will show that a two of the three states are indeed degenerate.

As usual before starting a computation is a good idea to look at scale stability:

alphaRange=Range[0.5,5,0.1];

eneA = Table[
Re[ Sort[Eigenvalues[hamHH[α],-dim,Method-> Arnoldi]]][[1;;Min[6,dim]]],{α,alphaRange}];

ListPlot[Table[Transpose[{alphaRange,eneA[[All,j]]}],{j,1,6}],
Joined-> True]

We see a good stability region for α around 1.5, and this will be our preferred value. The region will grow augmenting the basis's dimension of course.

Eigenvalues

We will need matrices at various dimensions. It is simpler to write a function which collect all the above tensor products. As the reader can imagine we want to have a freedom to fix parity with respect inversion of y coordinates, in this case y-matrices are computed with only even or odd states. The following function collects tensor products and allows a choice of parity (this will be used in the next subsection):

getMatricesHamiltonian[dim1_,dim2_,minV_,par2_:-1]:= Module[
{p2x,x2,x4,x1,x3,kin,V2,V4,V3,V0,
p2y,y2,y4,y1},
(* compute x matrices *)
{p2x,x2,x4,x1}=hermiteMatrices[dim1+4];
x3 = x2.x1;
{p2x,x2,x4,x1,x3}=Take[#1,{1,dim1},{1,dim1}]&/@N/@{p2x,x2,x4,x1,x3};
𝒰1x =SparseArray@IdentityMatrix[dim1];
(* compute y matrices *)
If[par2 == 0 || par2 == 1,
{p2y,y2,y4}=N/@hermiteMatrices[dim2,par2];,
{p2y,y2,y4,y1}=N/@hermiteMatrices[dim2];];
𝒰1y = SparseArray@IdentityMatrix[dim2];
kin = (1/2) p2x⊗𝒰1y +(1/2) 𝒰1x⊗p2y;
     V2 = β1 (x2⊗𝒰1y + 𝒰1x⊗y2);
     V4 =  α1 (x4⊗𝒰1y + 2 x2⊗y2+𝒰1x⊗y4);
V3 =  γ1 (x3⊗𝒰1y - 3 x1⊗ y2 );
V0 =  -minV 𝒰1x⊗𝒰1y;
{kin,V0,V2,V3,V4}]

We will use the same symbols V2, V3 etc. throughout then the definition of hamHH[α] will be at work.

Let us do a computation with a bigger basis, dim = 50, i.e. 2500 basis functions :

dim=50; α = 1.5;

{kin,V0,V2,V3,V4}= getMatricesHamiltonian[dim,dim,minV];

Dimensions[hamHH[1.5]]

The Hamiltonian is quite huge, but it is very sparse :

ham =hamHH[1.5];
ArrayPlot[ham,ColorFunction->(If[#==0,White,Black]&),ColorFunctionScaling->False,
FrameTicks->{{Range[0,Length[ham],Floor[Length[ham]/5]],None},{Range[0,Length[ham],Floor[Length[ham]/5]],None}}]

Here are the first eigenvalues

ene1= Re[Sort[Eigenvalues[hamHH[α],-100,Method-> Arnoldi]]];//Timing

ene1[[1;;27]]

We see the triplet structure and the predicted degeneracy. Let us check that the situation is stable against dimensions

dim = 70;
{kin,V0,V2,V3,V4}= getMatricesHamiltonian[dim,dim,minV];
Dimensions[hamHH[1.5]]

ene1dim70= Re[Sort[Eigenvalues[hamHH[α],-100,Method-> Arnoldi]]];//Timing

The two values are in reasonable agreement :

GraphicsRow[ListPlot[Abs[ene1dim70-ene1][[1;;#1]]]&/@{100,30},ImageSize-> 8 72]

Things go better clearly for higher dimensions. With dim = 100, i.e. a basis of (10^4)elements,

dim = 100;
{kin,V0,V2,V3,V4}= getMatricesHamiltonian[dim,dim,minV];
Dimensions[hamHH[1.5]]

ene1dim100= Re[Sort[Eigenvalues[hamHH[α],-100,Method-> Arnoldi]]];//Timing

GraphicsRow[ListPlot[Abs[ene1dim70-ene1dim100][[1;;#1]]]&/@{100,30},ImageSize-> 8 72]

For our pedagogical needs the use of dim = 50 will be in general sufficient.

Eigenstates and parity

As stated above the degeneracy has not been removed, so our eigenstates are not completely characterized, except for non degenerate levels. It is not difficult to construct a parity operator which implement the inversion of y coordinate.  This can be used to select parity eigenstates from our basis.

Each Hermite polynomial has parity (((-1))^deg)where deg is its degree. Our counting of states start from 1, i.e. degree 0 corresponds to the first state. Then y inversion is trivially written as (in y coordinates):

dim = 50;
inversionY = SparseArray@DiagonalMatrix[Table[(((-1))^(i+1)),{i,1,dim}]];

As an operator in tensor space :

𝒫y = (SparseArray@IdentityMatrix[dim])⊗ inversionY;

Let us first compute eigenvectors.

{kin,V0,V2,V3,V4}= getMatricesHamiltonian[dim,dim,minV];
Dimensions[hamHH[1.5]]

{eigH,eigV} = Eigensystem[hamHH[α],-100,Method-> Arnoldi];
order = Ordering[eigH];
eigH= eigH[[order]];
eigV = eigV[[order]];

Parity eigenstates must satisfy 𝒫y.v = v or 𝒫y.v = -v. The following function tests the two possibilities computing the norm of the difference between left hand side and right hand side, printing the most favorable, i.e. the lowest norm.

getParity[v_]:= Module[{ris},ris=Norm[#1]&/@{𝒫y.v-v,𝒫y.v+v};
If[ris[[1]]< ris[[2]],Row[{1,"; prec=",ris[[1]]}],Row[{-1,"; prec=",ris[[2]]}]]]

For the first 15 states, for instance

Partition[getParity[eigV[[#1]]]&/@Range[1,18],3]

Comparing with eigenvalues we see that in fact the isolated states are, numerically,  eigenstates of the parity operator, while "doublets" suffer a certain mixing. The reader must note the states number 12 and 18: these are negative parity states, and their energy is higher than energy of the states of the corresponding doublet, as it can be checked from the table below:

Partition[eigH[[1;;18]],3]//TableForm

In the model Hamiltonian () this corresponds to δ < 0.

Let us construct explicitly even and odd states for the first excited level. We must have, apart from an overall normalization:

vecEven = c2 (Subscript[V, 2]) + (Subscript[V, 3]) ;  𝒫y (c2 (Subscript[V, 2]) + (Subscript[V, 3])) = c2 (Subscript[V, 2]) + (Subscript[V, 3])

vecOdd = (Subscript[V, 2]) + c3 (Subscript[V, 3]);     𝒫y ((Subscript[V, 2]) + c3 (Subscript[V, 3]))= - ((Subscript[V, 2]) + c3 (Subscript[V, 3])).

The two equations are highly nontrivial, as they are, in our case, 2500 equations for a single variable! They will be only approximately true due to numerical approximations.

We can find the coefficients by minimizing the norms :

minc2=FindMinimum[Norm[𝒫y.(cvar eigV[[2]]+eigV[[3]])- (cvar eigV[[2]]+eigV[[3]])],{cvar,0,1}]

minc3=FindMinimum[Norm[𝒫y.(cvar eigV[[3]]+eigV[[2]])+ (cvar eigV[[3]]+eigV[[2]])],{cvar,0,1}]

c2=cvar/.Last[minc2];
c3 = cvar/.Last[minc3];

vecEven = c2*eigV[[2]]+ eigV[[3]]; vecEven = vecEven/Norm[vecEven];
vecOdd = eigV[[2]]+ c3 eigV[[3]]; vecOdd = vecOdd/Norm[vecOdd];

We can verify that the linear combinations continue to be eigenstates of the Hamiltonian, within errors of course.

Norm[hamHH[α].#1 - eigH[[3]]#1]&/@{vecEven,vecOdd}

Eigenfunctions plots

We can repeat the procedure used for the harmonic oscillator. First let us look at the ground state:

xg = Range[-5,5,0.2]; 
yg = Range[-5,5,0.2];

data = makeData[α,xg,yg,eigV[[1]],dim,dim,{-1,-1}];

ListPlot3D[data,PlotRange->All,Mesh->None,AxesLabel->{"x","y",""},Boxed->False]

We see the triangular symmetry and the even shape under y → - y. Even more interesting are graphs for even and odd states

data = makeData[α,xg,yg,vecEven,dim,dim,{-1,-1}];
ListPlot3D[data,PlotRange->All,Mesh->None,AxesLabel->{"x","y",""},Boxed->False]

data = makeData[α,xg,yg,vecOdd,dim,dim,{-1,-1}];
ListPlot3D[data,PlotRange->All,Mesh->None,AxesLabel->{"x","y",""},Boxed->False]

A more detailed look at  higher states will be given in the next section.

Fixed parity computations

The reader has no doubt noticed that  it is quite ineffective to compute all the spectrum as degenerate states are not separated. It is clear that we can separate a priori the states of different parity simply changing the set of base functions: it is enough to use even hermite polynomials in y for even states and odd polynomials for odd states. This will also improve the precision, as at fixed total number of states a larger portion of Hilbert space is covered.

All the computations are similar to the previous ones, we have simply to specify parity. Let us choose for example even states.

α1=1;β1=- 1/2;γ1 = -3;
minV = First[FindMinimum[α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1(x^3-3 x (y^2)),{x,y}]];
dim1 = 50; dim2=50;
parity = 0;
{kin,V0,V2,V3,V4}=getMatricesHamiltonian[dim1,dim2,minV,parity];

eneEven = Re[Sort[Eigenvalues[hamHH[α],-100,Method-> Arnoldi]]];
eneEven[[1;;27]]

We see the absence of degeneracy. The first levels of both parity are:

ene1[[1;;20]]

and we have a good agreement between the energy levels.

Wave functions

Let us compute the eigenvectors and show some different ways for their visualization.

dim1 = 50; dim2=50;
parity = 0;
{kin,V0,V2,V3,V4}=getMatricesHamiltonian[dim1,dim2,minV,parity];
{eigH,eigV} = Eigensystem[hamHH[α],-150,Method-> Arnoldi];
order = Ordering[eigH];
eigH= eigH[[order]];
eigV = eigV[[order]];

One can use three dimensional plots, remembering to specify parity. We do plots at higher resolution than before, so the process will take some more time

xg = Range[-5,5,0.1]; 
yg = Range[-5,5,0.1];

We plot the absolute value of ψ .

data = makeData[α,xg,yg,eigV[[55]],dim1,dim2,{-1,parity}];
ListPlot3D[absData[data], PlotRange -> All, Mesh -> None, AxesLabel -> {"x", "y", ""}, Boxed -> False]

A second possibility is to use ListContourPlot. Moving with the mouse on the graphs you will see the various level of  |ψ|.

data = makeData[α,xg,yg,eigV[[55]],dim1,dim2,{-1,parity}];
ListContourPlot[absData[data],
ColorFunctionScaling->True,
ContourStyle->{{GrayLevel[0],Thickness[0.002]}}]

A third, faster possibility, is to use directly ArrayPlot or ReliefPlot. We have to extract from makeData the values of the function, and this can be done taking the component 3, 6, 9 of the flattened matrix, as makeData produce {x,y,f} lists. Thne we use partition to divide data following the number of grid points;  finally we plot the transpose of the result, as ArrayPlot and ReliefPlot are row graphs and interpret column number as horizontal coordinate, while in our scheme column number denote y coordinate. We plot probabilities as the graph is more significant

data = makeData[α,xg,yg,eigV[[55]],dim1,dim2,{-1,parity}];
data1 =Partition[Flatten[ data][[3;; ;;3]],Length[yg]];

ReliefPlot[Transpose[data1^2]]

We can encode all the previous stuff in

onlyData[α_,xg_,yg_,v_,dim1_,dim2_,{parx_,pary_}]:= Module[{data},
data =  makeData[α,xg,yg,v,dim1,dim2,{parx,pary}];
Transpose[Partition[Flatten[ data][[3;; ;;3]],Length[yg]]]]

Here is an example of quite complex structures

ArrayPlot[onlyData[α,xg,yg,eigV[[#1]],dim1,dim2,{-1,parity}]^2]&/@(70+Range[0,5])

The reader can play with parameters and in particular study the sensitivity of the wave function to small changes in the parameters α1, β1, γ1.

Clear[data,data1,"Global`eigH*","Global`eigV*",kin,V0,V2,V3,V4];

Approximate spectral method

Introduction

In the general case analytical matrix elements of the potential are not available and we have to resort to numerical evaluation of integrals. We will use the Gaussian integration formulas based on zeros of Hermite polynomials. A generic matrix element has the form

V[[((Subscript[k, 1]),(Subscript[k, 2])),((Subscript[j, 1]),(Subscript[j, 2]))]]=(∫_(-∞)^∞)𝕕x 𝕕y f[(Subscript[k, 1]);x]f[(Subscript[k, 2]);y]V[x,y] f[(Subscript[j, 1]);x]f[(Subscript[j, 2]);y]

f' s are the basis functions, in our case Exp[-x^2/2] h[k;x]  where h[k;x] is the normalized form of Hermite polynomials. The application of Gauss formulas to () is straightforward. Let ((Subscript[x, a]),(Subscript[w, a])) the set of roots and weights for x coordinates, and ((Subscript[y, b]),(Subscript[w, b]))a similar set for y coordinates, then () is approximated by

V[[((Subscript[k, 1]),(Subscript[k, 2])),((Subscript[j, 1]),(Subscript[j, 2]))]]= (∑_(a,b))(Subscript[w, a])(Subscript[w, b]) h[(Subscript[x, a]),k1]h[(Subscript[y, b]),k2]V[(Subscript[x, a]),(Subscript[y, b])] h[(Subscript[x, a]),j1]h[(Subscript[y, b]),j2]

To compute matrix elements and to write Hamiltonian we have to solve a couple of practical problems:

We have to fix a correspondence  (m,n) → I between pairs of indices and one index I, in such a way that matrices can be written in usual form as H[[I,J]] .

We have to write expression () in a vectorized form, in such a way that listability properties of Mathematica functions can be used, otherwise the execution's times will be not realistic.

In the two following subsections we give some details on the functions getp2 and getMatV which compute matrix elements for kinetic term and potential respectively, the reader not interested in these details can skip the paragraphs and go to section Simple applications, where some examples are given.

Indexes

There are several ways to trade a pair (or more) of indices for one index. If i , j  run from 1 to dim1, dim2 respectively one possibility is to write

I = i + dim1 (j-1)

This is equivalent to writing the sequence Flatten - Partition to the original two index structure. This possibility will be used for the Hylleraas variables in a variational calculation of the helium spectrum. A second possibility is to use a "circular" counting, i.e. fix a priori a maximum index, ω, and consider all pairs such that  i + j ≤ ω. This possibility will be used with the Pekeris variables for the helium problem and is the one adopted in this notebook. For the harmonic oscillator type of basis this procedure amounts to put a cutoff on the total "energy" of the basis states, as the energy is proportional to the index. The first possibility instead puts a cutoff on the individual energies.

The practical implementation we will use is quite trivial: make a double loop on i and j;  at each stage increment the index I only if the constraint is satisfied, that's all. We will also need the converse: given an index I, write the pair (i,j) which originates it. This is done by a vector, in our routine named fromInd which contains for each successful test the pair {i,j} which has passed the test. There exist certainly more sophisticated and efficient ways to produce the same result but we prefer to keep  things simple at this stage.

The function contained in the section functions,  needs as input the maximum index we want and produces a triplet: {dimIndR, indR, fromInd}, the first argument being the dimension of the base, the second an array which contains indices in the form ind[i,j] = =I.  We will not use this array, and the third is a vector formInd, such that fromInd[[I]] = {i,j}. Clearly dimIndR = Length[fromInd].

A couple of examples can help. Put a cutoff of 3, i.e. we accept only indices such that their sum is less than 3.  We have only three possible pairs, (1,1), (1,2) and (2,1). Here is how it works:

{dimIndR,indR,fromInd}= funzIndice[3];

fromInd

fromInd[[2]]

dimIndR

Matrix elements

Let us start from the kinetic term. We know exactly the matrix elements of p^2 for each dimension. We can write tensor products of each of them,

p2Tensor= p2x⊗𝟙 + 𝟙⊗p2y

and all our elements are written inside this matrix. We have only to find a quick way to extract them. This is done by the following trick. Construct an index in the form

index = Table[{i1,j1} = fromInd[[ind1]];i1+dim (j1-1),{ind1,1,dimIndR}];

This essentially transforms the circular ordering in the ordering adapted to Kronecker products. The matrix we look for is selected by taking rows and columns enumerated by index:   p2Tensor[[index, index]].  We do not expect that the reader accepts this on faith, or that it is simple to visualize. Below some checks of the procedure will be provided.

Let us now consider the potential terms.  We report here the expression to be computed, which we rewrite moving  x - variables to the same side of the sum

V[[((Subscript[k, 1]),(Subscript[k, 2])),((Subscript[j, 1]),(Subscript[j, 2]))]]= (∑_(a,b))(Subscript[w, a]) h[(Subscript[x, a]),k1]h[(Subscript[x, a]),j1]V[(Subscript[x, a]),(Subscript[y, b])] (Subscript[w, b]) h[(Subscript[y, b]),k2]h[(Subscript[y, b]),j2]

The problem is clear : we want to sum using matrix multiplication on the Gauss grid. We can write easily a matrix V[(Subscript[x, a]),(Subscript[y, b])]by something like

VontheGrid = Table[V[GaussGridx[[i]] , GaussGridy], {i, 1, Length[GaussGridx]}];

In this way we need two vectors respectively on x - grid and y - grid

v1[(Subscript[x, a])]= (Subscript[w, a]) h[(Subscript[x, a]),k1]h[(Subscript[x, a]),j1];     v2[(Subscript[y, b])] = (Subscript[w, b]) h[(Subscript[y, b]),k2]h[(Subscript[y, b]),j2].

This for each pair, to be obtained with the function fromInd. The problem is that initial and final indexes are entangled, as the x grid vector contains both j and k indexes, the same for y grid vector. We can compute each vector inside a double loop, on big indices K, J for the initial and final state, extract for each index the pairs ((Subscript[k, 1]),(Subscript[k, 2]))= fromInd[[K]], and ((Subscript[j, 1]),(Subscript[j, 2]))= fromInd[[J]] and make the final dot product

(∑_(a,b)) v1[(Subscript[x, a])]V[(Subscript[x, a]),(Subscript[y, b])]v2[(Subscript[y, b])]

This can be done but turns out to be quite slow, the reader can try himself. We adopted another strategy, not very efficient, maybe,  but faster. We compute first all polynomials on the grid. This is easy as our function hermiteP does exactly this. On the left hand side of () we have two such matrices: hPx= h[(Subscript[x, a]),k1] and  hPxW = (Subscript[w, a]) h[(Subscript[x, a]),j1],  the second one is the first equipped with weights. Then we take the tensor product of the two with

xvalues = Flatten[Outer[Times, hPx, hPxW, 1], 1];

We have a big matrix which contains grid values for all combinations of indexes (Subscript[k, 1]) and (Subscript[j, 1]). This is more than we need as (Subscript[k, 1])and (Subscript[j, 1])appear in less combinations, due to our circular ordering. The same is done on the right hand side, obtaining a corresponding matrix yvalues. Finally a dot product is made:

matV1 = xvalues.VontheGrid.Transpose[yvalues];

This explicitly performs the sums over the Gauss grids.  matV1 contains all our matrix elements but in the form of tensor product indices.  We have to extract them in terms of our big indexes K, J. This is particularly tricky as the tensor product on the left hand side is the tensor product of ((Subscript[k, 1])(Subscript[j, 1]))and similarly the right hand side, i.e. indices are "entangled".  We have separated the indices in the following way. We construct a vector of couples of indices

inda = Table[fromInd[[ind1]], {ind1, 1, dimIndR}];

We do a loop on the second (big) index ind2, and inside this we construct a second vector, indb which disentangle components by reporting them to a tensor like structure:

Table[
  {i2, j2} = fromInd[[ind2]];
  indb = Transpose[{inda[[All, 1]] + dim (i2 - 1), inda[[All, 2]] + dim (j2 - 1)}];
  ff /@ indb,
  {ind2, 1, dimIndR}];

ff is an easy way to write matV1[] on the indices generated with indb (see the code of the function).  dim is the maximum dimension chosen with funzIndice.

This tour de force on the indices exceeds by far the author's intuitive grasp of the subject and the procedure outlined above is the result of successive improvements on a basic "loop- and - go" strategy.

Two more remarks:

If variational-type calculations a scale factor α is needed.  This is simply included by computing the potential V[x,y] on a rescaled grid, V[x/α , y/α].

It is useful to have the possibility to switch to fixed parity functions in one or both dimensions. This affects only the choice of base functions, nothing else is changed, and the choice is included in the function getMatV.

Checks

loopgetMatV[dim_,{dimGauss1_,dimGauss2_},potV_,α_]:= Module[
{GaussGridx,weightsx,GaussGridy,weightsy,
hPx,hPxW,hPy,hPyW,VontheGrid,matV,dimIndR,indR,fromInd},
(* compute grid and weights for Gaussian integration *)
{dimIndR,indR,fromInd}=funzIndice[dim];
{GaussGridx,weightsx} = hermiteRoots[dimGauss1];
{GaussGridy,weightsy} = hermiteRoots[dimGauss2];
(* normalized Hermite polynomials without and with weights *)
(* The format is  h[i,x] *)
hPx=hermiteP[dim,GaussGridx];
hPxW =Transpose[Transpose[hPx]*weightsx];
hPy=hermiteP[dim,GaussGridy];
hPyW =Transpose[Transpose[hPy]*weightsy];
(* compute potential on Gauss grid points, the format is V[x,y] *)
VontheGrid= Table[potV[GaussGridx[[i]]/α , GaussGridy/α],{i,1,Length[GaussGridx]}];
(* compute off diagonal matrix elements of the potential V[k,j] *)
matV = ConstantArray[0,{dimIndR,dimIndR}];
(* kx = fromInd[[ind1,1]],  jx= fromInd[[ind2,1]] etc *)
Table[matV[[ind1,ind2]]=
(hPx[[fromInd[[ind1,1]]]]*hPxW[[fromInd[[ind2,1]]]]).VontheGrid.
(hPy[[fromInd[[ind1,2]]]]*hPyW[[fromInd[[ind2,2]]]]),{ind1,1,dimIndR},
{ind2,ind1+1,dimIndR}];
matV = matV + Transpose[matV];
Table[matV[[ind1,ind1]]=
(hPx[[fromInd[[ind1,1]]]]*hPxW[[fromInd[[ind1,1]]]]).VontheGrid.
(hPy[[fromInd[[ind1,2]]]]*hPyW[[fromInd[[ind1,2]]]]),{ind1,1,dimIndR}];
matV]

checkTensor[A_,B_]:=  Module[{matp2,dimIndR,dim,
indR,fromInd,unitM},
dim = Length[A];
{dimIndR,indR,fromInd}=funzIndice[dim];
(* kx = fromInd[[ind1,1]],  jx= fromInd[[ind2,1]] etc *)
matp2 = ConstantArray[0,{dimIndR,dimIndR}];
Table[matp2[[ind1,ind2]]=
A[[fromInd[[ind1,1]],fromInd[[ind2,1]]]]*B[[fromInd[[ind1,2]],fromInd[[ind2,2]]]];
,{ind1,1,dimIndR},
{ind2,ind1+1,dimIndR}];
matp2 = matp2 + Transpose[matp2];
Table[matp2[[ind1,ind1]]=
A[[fromInd[[ind1,1]],fromInd[[ind1,1]]]]*B[[fromInd[[ind1,2]],fromInd[[ind1,2]]]];
,{ind1,1,dimIndR}];
matp2];

dim=4; dimGauss1=50; dimGauss2= 48;
α = 1.;

potV[x_,y_]:= (1/2)x^2 + (1/4)y^2 + x y^2;

{p2,x2,x4,x1}=N/@ hermiteMatrices[dim];
𝒰1 = IdentityMatrix[dim];

Checks for kinetic term

checkp2=(checkTensor[p2,𝒰1]+checkTensor[𝒰1,p2]);
matp2 =N@ Normal[getp2[dim,{-1,-1}]];

compare[a_,b_]:= Union[Thread[Chop[Flatten[a-b]]== 0]]

compare[checkp2,matp2]

Checks for the potential term

check1potV = (1/(2(α^2)))checkTensor[x2,𝒰1] + (1/(4 (α^2)))checkTensor[𝒰1,x2] + (1/(α^3))checkTensor[x1,x2];

check2potV = loopgetMatV[dim,{dimGauss1,dimGauss2},potV,α];

matV = getMatV[dim,{-1,-1},{dimGauss1,dimGauss2},potV,α];

MapThread[compare,{{check1potV,check2potV},{matV,matV}}]

Times

loopgetMatV[30,{dimGauss1,dimGauss2},potV,α];//Timing

{Timing[getMatV[20,{-1,-1},{dimGauss1,dimGauss2},potV,α];],
Timing[loopgetMatV[20,{dimGauss1,dimGauss2},potV,α];]}

A simple application

As a first application let us compare the exact spectral method of previous sections to this approximate one in the case of Henon-Heiles -like potential.

Exact spectral method

α1=1;β1=- 1/2;γ1 = -3;
minV = First[FindMinimum[α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1(x^3-3 x (y^2)),{x,y}]];

α = 1.5;
dim = 50;

({kin,V0,V2,V3,V4}= getMatricesHamiltonian[dim,dim,minV];
eneExactSM=Sort[Eigenvalues[hamHH[α],-100,Method-> Arnoldi]];)//Timing

Dimensions[V0]

Here is the computation with approximate evaluation of matrix elements. We choose 100 points for Gaussian integration.

dimGauss1=100; dimGauss2=100;

dim=50;

parity=-1; (* all states *)

potV[x_,y_]:= α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1(x^3-3 x (y^2)) - minV;

(p2 = getp2[dim,{parity,parity}];
matV=getMatV[dim,{parity,parity},{dimGauss1,dimGauss2},potV,α];
ham = ((α^2)/2)p2 + matV;
eneApproxSM=Sort[Eigenvalues[ham,-100,Method-> Arnoldi]];)//Timing

We see that circular ordering has limited the number of states to 1225, compared to 2500 used with dim = 50 in previous calculation. A comparison of the two methods follows:

NumberForm[eneExactSM[[1;;21]],10]

NumberForm[eneApproxSM[[1;;21]],10]

GraphicsRow[ {ListPlot[(eneApproxSM-eneExactSM)[[1;;30]] ],
ListPlot[eneApproxSM-eneExactSM ]},ImageSize-> 8 72]

We see a very good agreement for low levels and a somewhat larger discrepancy for higher states. In this model this is due to the difference in dimensions, in fact for dim =70 we have a comparable number of base function and an excellent agreement on the eigenvalues:

dim=70;

(p2 = getp2[dim,{parity,parity}];
matV=getMatV[dim,{parity,parity},{dimGauss1,dimGauss2},potV,α];
ham = ((α^2)/2)p2 + matV;
eneApproxSM=Sort[Eigenvalues[ham,-100,Method-> Arnoldi]];)//Timing

GraphicsRow[ {ListPlot[Abs[eneApproxSM-eneExactSM][[1;;30]] ],
ListPlot[Abs[eneApproxSM-eneExactSM ]]},ImageSize-> 8 72]

Fixed parity

dim1 = 50; dim2=50;
parity = 0;
{kin,V0,V2,V3,V4}=getMatricesHamiltonian[dim1,dim2,minV,parity];

eneEvenESM = Re[Sort[Eigenvalues[hamHH[α],-100,Method-> Arnoldi]]];

dim=50;

(p2 = getp2[dim,{-1,parity}];
matV=getMatV[dim,{-1,parity},{dimGauss1,dimGauss2},potV,α];
ham = ((α^2)/2)p2 + matV;
eneEvenASM=Sort[Eigenvalues[ham,-100,Method-> Arnoldi]];)//Timing

GraphicsRow[ {ListPlot[Abs[eneEvenASM-eneEvenESM][[1;;30]] ],
ListPlot[Abs[eneEvenASM-eneEvenESM ]]},ImageSize-> 8 72]

A second example

The approximate method we are using is more flexible, and it allows to consider any kind of potential. The only warning is to use even dimensions for a potential singular at the origin, otherwise the point (0,0) would be part of the Gauss grid and indefinite expressions like 0/0 would appear in the Gauss sum.

For example let us consider a generalization of previous potential, having a 5-fold symmetry. In radial coordinates it is written as

V[r,θ] = (Subscript[α, 1]) (r^4) + (Subscript[β, 1])(r^2) + (Subscript[γ, 1]) Cos[5θ] (r^3)

To write the potential in cartesian coordinates we can use the Exp[i θ] = (x + i y)/r and with

Simplify[(1/(2((x^2)+(y^2))))Expand[(((x + I y))^5)+ (((x- I y))^5)]]

we have

V[x,y] = α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1((((x^5)-10 (x^3) (y^2)+5 x (y^4)))/((x^2)+(y^2)))

The computation is exactly the same as before. We expect at low energies for the usual values of the parameters, five equivalent minima.

α1=1;β1=- 1/2;γ1 = -3;
minV = First[FindMinimum[α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1((x^5)-10 (x^3) (y^2)+5 x (y^4))/(x^2+y^2),{x,y}]]

potV[x_,y_]:= α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) +  γ1((x^5)-10 (x^3) (y^2)+5 x (y^4))/(x^2+y^2)-minV

Plot3D[potV[x,y],{x,-5,5},{y,-5,5},PlotRange-> {Automatic,Automatic,{0,25}}]

The model Hamiltonian describing the system is

H = (({{(Subscript[ϵ, 0]), δ, 0, 0, δ}, {δ, (Subscript[ϵ, 0]), δ, 0, 0}, {0, δ, (Subscript[ϵ, 0]), δ, 0}, {0, 0, δ, (Subscript[ϵ, 0]), δ}, {δ, 0, 0, δ, (Subscript[ϵ, 0])}}))

where  δ stands for the tunneling amplitude between adjacent minima. This Hamiltonian has a spectrum

E = (Subscript[ϵ, 0]) + 2 Cos[2π k (δ/5)];  k = 0,1,2,3,4

two pairs of degenerate doublets. Let us verify this on the numerical computation

dim=50;
parity=-1;

p2 = getp2[dim,{parity,parity}];
matV=getMatV[dim,{parity,parity},{dimGauss1,dimGauss2},potV,α];

ham = ((α^2)/2)p2 + matV;

ene5 =Re[ Sort[Eigenvalues[ham,-dim,Method-> Arnoldi]]];

Here we print eigenvalues, 5 in each row to emphasize degenerate doublets in a complex of five states, as expected.

Partition[ene5,5]//TableForm

A bit different is the case of cubic term with an even number of minima, as Cos[4θ] (r^3). In this case the spectrum of the model Hamiltonian is

E = (Subscript[ϵ, 0]) + 2 Cos[2π k (δ/4)];  k = 0,1,2,3

and we have only one doublet and two non degenerate levels. Again this is easily checked in the model.

First we translate the potential in cartesian coordinates using

Simplify[(1/(2(Sqrt[((x^2)+(y^2))])))Expand[(((x + I y))^4)+ (((x- I y))^4)]]

α1=1;β1=- 1/2;γ1 = -3;
minV = First[FindMinimum[α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1(((x^4)-6 (x^2) (y^2)+(y^4))/(Sqrt[(x^2)+(y^2)])),{{x,10},{y,0}}]]

This is the potential and its graph

potV[x_,y_]:= α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1(((x^4)-6 (x^2) (y^2)+(y^4))/(Sqrt[(x^2)+(y^2)]))-minV

Plot3D[potV[x,y],{x,-5,5},{y,-5,5},PlotRange-> {Automatic,Automatic,{0,25}}]

The computation is always the same

dim=50;
parity=-1;

p2 = getp2[dim,{parity,parity}];
matV=getMatV[dim,{parity,parity},{dimGauss1,dimGauss2},potV,α];

ham = ((α^2)/2)p2 + matV;

ene4 =Re[ Sort[Eigenvalues[ham,-dim,Method-> Arnoldi]]];

Here we print the quartet of levels and a single degenerate doublet is clearly visible

Partition[ene4,4]//TableForm

Clear[data,data1,"Global`eigH*","Global`eigV*",kin,V0,V2,V3,V4,ham,p2,matV];
Clear[x1,x2,x3,x4];

Pseudospectral Method

General remarks

As a last method we give a brief account of the collocation or pseudospectral method. As explained in notebooks [NB-10.4-PSmethods.nb],  this method approximates the equation directly on a space grid, built with the zeros of the Chebyshev polynomials, i.e. the Gauss grid is at the same time the grid of points. The Schrödinger equation takes the form

-(1/2)(Subscript[D, IJ]) (Subscript[ψ, J]) + V[(Subscript[P, I])](Subscript[ψ, I]) = E (Subscript[ψ, I])

where D is the spectral derivative. In one dimension we have used the functions provided by Mathematica to construct the derivative operator, in two dimensions it is sufficient to take the tensor product of the one dimensional result. The function getD2 used in notebook [13] provides the one dimensional derivative as well as ξ grid, in Chebyshev variable running from -1 to +1 and the x grid, without rescaling, ranging from -∞ to +∞.

The only real problem, as usual, is to connect the double index (i,j) of the grid to a vector-like index I, to use in (). Here we use a plain index procedure, without circular cut, in the form

I = i + dim (j-1)

From index I position i and j are easily extracted as

i = Mod[I-1,dim]+1;    j = Quotient[I-1,dim]+1.

The reader can convince herself (or himself) of the correctness of this relation.  It is easier to consider indices starting form 0.

The tensor product D ⊗ 𝟙 + 𝟙 ⊗ D gives directly the matrix order with indexes I,J. For the potential term one has to construct a vector (Subscript[V, I]). Using a vectorized instruction (i.e. Table) one has to write

Table[V[xgrid, xgrid[[i]] ],{i,1,dim}]

One has to fill the vector with all x coordinates at each fixed y, which indeed is what is written in ().

Checks

Before doing any calculation let us check the index ordering on low dimensional matrices.

dim=5;

Here we get one dimensional spectral derivative and points grid.

{ξgrid,xgridL1,D2}= getD2[dim];

This is explicitly the vector which contains the conversion from (i, j) to I indexes

fIn=Table[{Mod[k,dim]+1,Quotient[k,dim]+1},{k,0,dim^2-1}];

As a reference potential we use

potV[x_,y_]:= x y^2;

First we make the tensor product of one dimensional derivatives

𝒰1 =SparseArray[ IdentityMatrix[dim]];
matD2=Normal[(D2⊗𝒰1+ 𝒰1⊗D2)];

and we fill the matrix checkD2 with explicit values of the operator, which in (i,j) coordinates reads:

⟨(Subscript[i, 1]),(Subscript[i, 2])|D2⊗𝒰1+ 𝒰1⊗D2|(Subscript[j, 1]),(Subscript[j, 2])⟩ = (Subscript[D2, (Subscript[i, 1])(Subscript[j, 1])])(Subscript[δ, (Subscript[i, 2])(Subscript[j, 2])])+ (Subscript[δ, (Subscript[i, 1])(Subscript[j, 1])])(Subscript[D2, (Subscript[i, 2])(Subscript[j, 2])])

checkD2 = ConstantArray[0,{dim^2,dim^2}];
Table[checkD2[[i,j]]= D2[[fIn[[i]][[1]],fIn[[j]][[1]]]]KroneckerDelta[fIn[[i]][[2]],fIn[[j]][[2]]]+
+D2[[fIn[[i]][[2]],fIn[[j]][[2]]]]KroneckerDelta[fIn[[i]][[1]],fIn[[j]][[1]]];,
{i,1,dim^2},{j,1,dim^2}];

For the potential we first construct the vector along the lines written above :

vPot = Flatten[Table[potV[xgridL1, xgridL1[[i]] ],{i,1,dim}]];

then explicitly fill a check vector with the values V[(Subscript[x, (Subscript[i, 1])]),(Subscript[y, (Subscript[i, 2])])]

checkPot = Table[{s1,s2} = fIn[[i]];
potV[xgridL1[[s1]], xgridL1[[s2]]],{i,1,dim^2}];

This is the result of the comparison

MapThread[compare,{{checkPot,checkD2},{vPot,matD2}}]

Applications

Harmonic oscillator

Let us start with a simple non isotropic oscillator, just to test the procedure

ω = 0.7;
potV[x_,y_]:= (1/2)x^2 + (ω^2) (1/2)y^2;

α=1.1;
dim = 20;

The derivative matrix

{ξgrid,xgridL1,D2}= getD2[dim];
𝒰1 =SparseArray[ IdentityMatrix[dim]];
matD2=(D2⊗𝒰1+ 𝒰1⊗D2);

and the diagonal matrix constructed with the potential:

matV = SparseArray@DiagonalMatrix@ Flatten[Table[potV[xgridL1/α, xgridL1[[i]]/α ],{i,1,dim}]];

The eigenvalues are

Sort[Eigenvalues[ -((α^2)/2)matD2 + matV,-20,Method-> Arnoldi]][[1;;10]]

to be compared to

eneH0 = Sort[Flatten[Table[((1/2)+n1) + ((1/2)+n2)ω,{n1,0,10},{n2,0,10}]]][[1;;10]]

We see that in this case probably quite large basis are needed to achieve a reasonable precision.

Quartic potential

Let us come back to our standard example

α1=1;β1=- 1/2;γ1 = -3;
minV = First[FindMinimum[α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1(x^3-3 x (y^2)),{x,y}]];

potV[x_,y_]:= α1 ((((x^2)+(y^2)))^2)+β1((x^2)+(y^2)) + γ1(x^3-3 x (y^2))-minV;

We choose a grid of 40 points per dimension, i.e.  1600 grid points and give a look at scale dependence, to choose a reasonable value for the scale parameter α.

dim = 40;

{ξgrid,xgridL1,D2}= getD2[dim];
𝒰1 =SparseArray[ IdentityMatrix[dim]];
matD2=(D2⊗𝒰1+ 𝒰1⊗D2);

alphaRange=Range[0.1,2,0.1];

eneA=Table[matV = SparseArray@DiagonalMatrix@ Flatten[Table[
potV[xgridL1/α, xgridL1[[i]]/α ],{i,1,dim}]];
Sort[Eigenvalues[ -((α^2)/2)matD2 + matV,-20,Method-> Arnoldi]][[1;;10]],{α,alphaRange}];

ListPlot[Table[Transpose[{alphaRange,eneA[[All,j]]}],{j,1,6}],
Joined-> True]

A stable plateau is around α = 0.5, which we adopt as the scale parameter. The reader has noted the quite big instability after α = 1, and the relative long execution time. This is due to the fact that derivative matrix is a full matrix and even the iterative Arnoldi method require plenty of time.

To do a high precision computation we take dim = 70, i.e. 4900 points on the grid. This computation requires several seconds on a 3 GHz computer, if the reader has a slow machine he can lower the parameter dim, or wait some minutes for the result

α = 0.5;

dim = 70;

{ξgrid,xgridL1,D2}= getD2[dim];
𝒰1 =SparseArray[ IdentityMatrix[dim]];
matD2=(D2⊗𝒰1+ 𝒰1⊗D2);

matV = SparseArray@DiagonalMatrix@ Flatten[Table[
potV[xgridL1/α, xgridL1[[i]]/α ],{i,1,dim}]];
(enePSM=Sort[Eigenvalues[ -((α^2)/2)matD2 + matV,-100,Method-> Arnoldi]];)//Timing

Here is the comparison with the spectral method results. The agreement is of the order of 1 part per million, which is not a bad result.

NumberForm[enePSM[[1;;21]],10]

NumberForm[eneExactSM[[1;;21]],10]

For the wave function one has to multiply by the relevant normalization factor. The function makeData can then be used to produce the plots. We leave this task to the interested reader.

This method is not the fastest among the ones we have presented but it is the simpler to encode.  We have used essentially only two lines: the definition of D2 and the definition of the vector for potential's values. The method can be very useful at least for a first analysis of the problems.

References

1) Trott, M. "The Mathematica Guidebooks Additional Material: Hénon-Heiles Eigenfunctions." ,

http://www.mathematicaguidebooks.org/additions.shtml#S_2_01.