Lo Surdo - Stark effect for Hydrogen atom

Quantum Mechanics
                      K. Konishi, G. Paffuti

Use parabolic coordinates to compute to arbitrary order in perturbation theory the energy shift due to an electric field on an arbitrary bound level of the Hydrogen atom. Do the computation in numerical and analytical form.

Introduction

In this notebook  atomic units will be used : |e| = ℏ = m = 1.

In problem [*] we have seen that the computation of perturbative corrections to atomic levels of Hydrogen in presence of a constant electric field can be carried out analytically using parabolic coordinates:

x = (Sqrt[ξ η]) Cos[φ]; y = (Sqrt[ξ η]) Sin[φ]; z = (1/2)(ξ-η);
ξ = r+z; η = r-z; φ=ArcTan[(y/x)];  r = (1/2)(ξ+η); 
Δ = (4/(ξ+η))(∂/(∂ξ))(ξ(∂/(∂ξ))) + (4/(ξ+η))(∂/(∂η))(η(∂/(∂η)))+ (1/(ξ η))((∂^2)/(∂(φ^2))).

The range for these coordinates is

0≤ (ξ, η ) ≤ ∞ ;   0≤ φ ≤ 2π

and volume element is

(𝕕^3)x = ((ξ+η)/4) 𝕕ξ 𝕕η 𝕕φ.

The electric interaction - e ℰ = |e| ℰ will be written as

V = (ℰ/2)(ξ - η).

In the rescaled variables

ϵ= (Sqrt[-2 E]);  x = ϵ ξ;  y = ϵ η;  (Subscript[β, i]) =((Subscript[Z, i])/ϵ);  ℱ = (ℰ/(ϵ^3)).

The eigenvalues problem was brought to the solution of a system of two differential equations, the energy eigenvalues of he original problem is found by the last constraint in the following equations:

(d/(d x))(x((Subscript[df, 1])/(d x)))+(-(1/4)x-((m^2)/(4x)) - (ℱ/4)(x^2)+ (Subscript[β, 1]))(Subscript[f, 1])= 0;
(d/(d y))(y((Subscript[df, 2])/(d y)))+(-(1/4)y-((m^2)/(4y)) + (ℱ/4)(y^2)+ (Subscript[β, 2]))(Subscript[f, 2])= 0;
(Subscript[β, 1])+ (Subscript[β, 2])= (Z/ϵ).

The leading order result (ℰ = 0) is

(Subscript[β, 1])= (Subscript[n, 1])+ ((|m|+1)/2);  (Subscript[β, 2])= (Subscript[n, 2])+ ((|m|+1)/2);       ⇒          (Subscript[β, 1])+(Subscript[β, 2]) = (Subscript[n, 1])+(Subscript[n, 2])+ (|m|+1) ≡ n;

then

(Subscript[β, 1])+(Subscript[β, 2]) = (Z/ϵ)       ⇒           ϵ = (Z/n)        ⇒        E = - (1/(2 (n^2)))

Normalized eigenfunctions for equations () are:

(Subscript[f, (Subscript[n, i]),m])[x] =(Sqrt[((n!)/((n+|m|)!))]) (x^(|m|/2)) Exp[-(x/2)](Subscript[L, n]^(|m|))(x)≡(Subscript[c, (Subscript[n, i]),m])  (x^(|m|/2)) Exp[-(x/2)](Subscript[L, n]^(|m|))(x)
(∫_0^∞)(Subscript[f, (Subscript[n, 1]),m])(Subscript[f, (Subscript[n, 2]),m])dx = (Subscript[δ, (Subscript[n, 1])(Subscript[n, 2])]).

Method of computation

The key point to do analytical computations lie in the fact that matrix elements of perturbation ((x^2))connect only neighboring states, this means that at any finite order a finite matrix computation will give exact results. This situation has already been exploited in notebook [NB-9.1] for polynomial perturbations to an harmonic oscillator, the strategy we will follow in this notebook is similar:

We will compute perturbation theory for equations () getting approximate values for eigenvalues (Subscript[β, 1])and (Subscript[β, 2]).

We will insert the obtained series in the constraint, determining energies.

As in notebook [NB-9.1] we will present two different approach: in the first we develop a numerical matrix calculation, in the second we build an operator approach, best suited to obtain analytical results.

Matrix approach

N.B. Matrices have indices which run from 1, then in this section the states follow the same convention, state number 1 will be the ground state etc. in the following section we come back to usual convention of naming quantum numbers starting from 0, in this case.

For simplicity in this section we will consider only the m = 0 case. The reader can extend results to other cases,the only thing to change are basis functions and their normalizations.

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

m=0;

Let us call ψB[i] the basis functions, first we have to construct the matrices which appear in he Hamiltonian. Let us give a name to quadratic combinations which enter in the Hamiltonian

D2ij[i_,j_]:= D[ψB[i],x] x D[ψB[j],x];
Vm1ij[i_,j_]:= ψB[i] (1/x)  ψB[j];
V1ij[i_,j_]:= ψB[i] x  ψB[j];
V2ij[i_,j_]:= ψB[i] x^2  ψB[j];

The basis functions are given by

ψB[i_]:= LaguerreL[i-1,x] Exp[-x/2];

Using recursion equations for Laguerre polynomials it is easy to show that only matrix elements (Subscript[M, ij]) with |i-j| ≤ 2 are different from zero, and all the integrals could be computed using these recursion relations, but we follow a more elementary approach, just computing integrals. It would be quite inefficient to use the general Integrate routine here: all functions are of the form (x^a)Exp[-x] and integrals are computed just by inspection using this function:

myInt[f_]:=  Together[Expand[f]//.{x^a_. Exp[-x]-> a!, Exp[-x]-> 1}];

Now we need a function which do integrations for previously defined combinations :

getMatrices[X_]:= Module[{mat},
mat = 0*IdentityMatrix[dim];
Table[mat[[i,j]]=myInt[X[i,j]];,{i,1,dim},{j,i,dim}];
mat + Transpose[mat] - DiagonalMatrix[Tr[mat,List]]];

Once defined the matrices the only thing we need to construct a recursive formulation of perturbation theory is the operator

Λ = Q (1/((Subscript[E, 0])- (Subscript[E, n])))Q

Where Q projects out of the original state, see notebook [*]. In "energy" eigenstates (here the role of energy is played by β parameters):

Gmatrix[k_,dim_]:= Module[{mat,axx},
Block[{nn},
nn = k - 1/2+m/2;
mat= DiagonalMatrix[1/(nn-(2*Range[dim]-1 +m)/2+axx)];
mat[[k,k]]=0; mat = mat/.axx-> 0]]

The value used, nn = k-1+m/2, is the value given above for unperturbed values of β, shifted by -1 because our indices start from 1.

The computation

As an example we compute energy corrections to ground state energy. Let us choose a relatively big dimension for our matrices:

dim=20;

and compute the necessary quadratic forms

matD2 = getMatrices[D2ij];
matV2 = getMatrices[V2ij];
matV1 = getMatrices[V1ij];

Matrix elements (quadratic forms) of equations () are:

∫ (D (Subscript[ψ, i]))x (D (Subscript[ψ, j])) + ∫(Subscript[ψ, i])(+(1/4)x+ ((m^2)/(4 x)) + (ℱ/4)(x^2) )(Subscript[ψ, j])=  β1 ∫(Subscript[ψ, i])(Subscript[ψ, j]);
 ∫ (D (Subscript[ψ, i]))y (D (Subscript[ψ, j])) + ∫(Subscript[ψ, i])(+(1/4)y+ ((m^2)/(4 y)) - (ℱ/4)(y^2) )(Subscript[ψ, j])=  β2 ∫(Subscript[ψ, i])(Subscript[ψ, j]);
β1 + β2 = 1/ϵ;

β1 e β2 play the role of eigenvalues. Let has define the two Hamiltonians (for m=0):

Ham1 = matD2 + (1/4)matV1+(ℱ/4)matV2;
Ham2 = matD2 + (1/4)matV1-(ℱ/4)matV2;

Before we start let us check that states (1,0,0...) , or states with 1 at k-th place are in effect eigenstates of H for zero external field:

test[k_]:= (temp = IdentityMatrix[dim][[k]];
Expand[(Ham1/.ℱ->0).temp] ==Expand[( k - 1/2 + m/2)temp])

Table[test[k],{k,1,dim/2}]

Now let choose our starting state (here the ground state)

state = 1;
ψ0 = Table[0,{dim}]; ψ0[[state]] = 1;

We compute the relative Gmatrix and the perturbation matrix

Λ = Gmatrix[1,dim];
v = (ℱ/4)matV2;

We can now write the starting point of our recursion relations, i.e. first order in perturbation theory

ψ[1] = Λ.v.ψ0;
λ[1] = ψ0.v.ψ0;

The recursion relations are :

ψ[n_]:=ψ[n]=Λ.v.ψ[n-1] - Sum[λ[n-k]Λ.ψ[k],{k,1,n-1}];
λ[n_]:=λ[n] = ψ0.v.ψ[n-1] ;

Let us fix the order we want to compute

pertOrder = 8;

and the unperturbed β value:

val = ((Ham1/.ℱ->0).ψ0)[[state]];

We ask for a list of eigenvalues and prepend the 0 order value to the list

β1=Prepend[Table[λ[k],{k,1,pertOrder}],val]

(Subscript[β, 2])is obtained simply by changing the sign in external field:

β2=β1/.ℱ-> -ℱ

To solve iteratively the constraint we write a series expansion for ϵ:

coeffeps =  Array[a,pertOrder];
eps = 1 + coeffeps.(ℰ^( Range[pertOrder]));

Now we put the constraint (Subscript[β, 1 ])+ (Subscript[β, 2])= 1/ ϵ

espr = Normal[Series[eps*(Tr[β1+β2]/.{ℱ-> ℰ/eps^3}),{ℰ,0,pertOrder}]]-1;

and compute the coefficients that satisfy this constraint:

coeffeps=coeffeps/.Solve[CoefficientList[espr,ℰ]== 0,coeffeps][[1]]

Table[a[i]=coeffeps[[i]];,{i,1,pertOrder}];

We can now compute the energy, E = - (ϵ^2)/2

eneMat =Expand[ -(1/2)Normal[Series[eps^2,{ℰ,0,pertOrder}]]]

The first two corrections coincide with those computed with Dalgarno - Lewis method, see problems  and notebook [NB-9.2].

Operator approach

Let us now consider the operator approach, which is even simpler than matrix approach. First of all let us switch to usual convention n =quantum numbers, they start from 0. Here we report for convenience of the reader eigenfunctions, eigenvalues and normalizations:

f[n,m,x]= c[n,m] (x^(|m|))Exp[-x/2]LaguerreL[n,|m|,x];
β = n + (1/2)+(m/2); n = 0,1,2,...;  m=0,1,...;
c[n,m] = (Sqrt[((n!)/((n+|m|)!))]);

To write a perturbation expansion we need to know how the operator (x^2)acts on the basis, or the operator x. This can be deduced from the recurrence relation

x L[n,α] = (n+α)(L[n,α]-L[n-1,α]) - (n+1)(L[n+1,α]-L[n,α]) =
- (n+1)L[n+1,α] - (n+α)L[n-1,α] + L[n,α] (2 n + α + 1)

Let us check this relation in a few cases:

Clear["n*","m*","λ*"];

espr1 = x LaguerreL[n,m,x] - 
( - (n+1)LaguerreL[n+1,m,x] - (n+m)LaguerreL[n-1,m,x] + 
LaguerreL[n,m,x] (2 n + m + 1) );
Simplify[{espr1/.{n->5,m->2},espr1/.{n->4,m->0},espr1/.{n->7,m->3}}]

Taking into account the normalizations coefficients c[n,m] it is quite simple to write the action of position operator 𝒬 in this representation:

Clear[n,mϕ,n1,n2]

𝓆[x_?NumberQ] = 0;
𝓆[x_.ket[n_]]:=-x  (n+1)(Sqrt[((n+mϕ+1)/(n+1))])ket[n+1] -x  (n+mϕ)(Sqrt[(n/(n+mϕ))])ket[n-1] + x(2n+mϕ+1)ket[n];
𝒬[x_]:= Distribute[𝓆[x]];

Unperturbed eigenvalues are β = n1 + (1/2)+(m/2); n1 = 0,1,2,...;  m=0,1,...; then the operator 1/(E- H0), projected on the space orthogonal to the level to be studied, is simply:

ℊ[k_][x_?NumberQ] = 0;
ℊ[k_][x_.ket[j_]]:=If[j===k,0,(x/((k-j)))   ket[j]];
𝒢[k_][x_]:= Expand[Distribute[ℊ[k][(x)]]];

We add an utility function, proj, which just project the result on the level to be studied:

proj[n_,espr_]:=Coefficient[espr,ket[n]];

Perturbation Hamiltonian is:

HI[x_]:= Collect[ (ℱ/4)*𝒬@(𝒬[x]),ket[_]]

For (Subscript[β, 2]) it is sufficient to change ℱ sign.  Function Collect is used to collect all expressions with  ket[k]. 
It is also useful to adopt some standard assumptions for Simplify, at the end of the section we resume standard assumptions to simplify factors in square roots:

old$Assumptions = $Assumptions;
$Assumptions = {{mϕ,n}∈Integers &&  n > 100,mϕ >100};

Example: first  order

This is first order in ℱ for (Subscript[β, 1]):

Simplify[proj[n,HI@ket[n]]]

Recursive definitions

Let us now write recursive relations of perturbation theory as operators equations:

psi[1] = 𝒢[n] @ HI@ket[n];

λ[1]=Simplify[proj[n,HI@ket[n]]];

psi[k_]:=psi[k]= Expand[𝒢[n] @ HI@psi[k-1] - Sum[λ[k-s]𝒢[n]@psi[s],{s,1,k-1}]]

λ[k_]:= λ[k]=Simplify[Expand[proj[n,HI@psi[k-1]]]]

λ[1]

Computation at  order 4

As an example let us compute fourth order correction to energy levels.

pertOrder = 4;

We ask for a list of eigenvalues and prepend the unperturbed value (the reader can also define λ[0] and start lists from 0 ).

The computation can take a few seconds

β1=Prepend[Table[λ[k],{k,1,pertOrder}],n+(1/2)+((mϕ)/2)]/.{mϕ-> m, n-> n1}

For (Subscript[β, 2])

β2 = β1/.{ℱ-> -ℱ, n1-> n2};

It is convenient to define principal quantum number n = n1+n2+1+m; and electric quantum number nE = n1-n2:

sost=Solve[{n1+n2+1+m == n,n1-n2 == nE},{n1,n2}][[1]]

Now we proceed as in the matrix case. At zero order  1/ϵ = n = n1+n2+1+m, then we can define

Clear["coeffeps*","a*"]
coeffeps =  Array[a,pertOrder];
eps = 1/n + coeffeps.(ℰ^( Range[pertOrder]));

Here is the constraint

espr = Normal[Series[eps*(Tr[β1+β2]/.{ℱ-> ℰ/eps^3}),{ℰ,0,pertOrder}]]-((1+m+n1+n2)/n);

We compute coefficients needed to satisfy constraint:

coeffeps=(coeffeps/.Solve[CoefficientList[espr,ℰ]== 0,coeffeps][[1]])/.sost;

Table[a[i]=coeffeps[[i]];,{i,1,pertOrder}];

and finally for the energy:

ene =Collect[Expand[ -(1/2)Normal[Series[eps^2,{ℰ,0,pertOrder}]]],ℰ,Simplify]

For ground state:

ene/.{n-> 1,nE-> 0,m-> 0}

to compare with previous result:

Normal[Series[eneMat,{ℰ,0,4}]]

$Assumptions = old$Assumptions;