Hydrogen atom in parabolic coordinates

Quantum Mechanics
                       K. Konishi, G. Paffuti

Check in Mathematica the solution in parabolic coordinates for Hydrogen atom obtained in problem 6.4.

Introduction

Parabolic coordinates (ξ, η, φ) are defined by

x = (Sqrt[ξ η]) Cos[φ]; y = (Sqrt[ξ η]) Sin[φ]; z = (1/2)(ξ-η);

ξ = r+z; η = r-z; φ=ArcTan[(y/x)];  r = (1/2)(ξ+η);

In this coordinates the eigenvalue problem takes the form, see problem 6.4

(d/(d ξ))(ξ((Subscript[df, 1])/(d ξ)))+((E/2)ξ-((m^2)/(4ξ)) + (Subscript[Z, 1]))(Subscript[f, 1])= 0;
(d/(d η))(η((Subscript[df, 2])/(d η)))+((E/2)η-((m^2)/(4η)) + (Subscript[Z, 2]))(Subscript[f, 2])= 0;
(Subscript[Z, 1])+ (Subscript[Z, 2])= Z.

The normalized eigenstates are, with  ϵ=(Sqrt[-2E]):

(Subscript[ψ, (Subscript[n, 1])(Subscript[n, 2])m])[ξ,η,φ]:= (Subscript[f, (Subscript[n, 1])m])[ξ] (Subscript[f, (Subscript[n, 2])m])[η]((e^(i m φ))/(Sqrt[2π]));
(Subscript[f, (Subscript[n, 1])m])[ξ]= (2^(1/4))ε (Sqrt[(( n!)/((n+|m|)!))]) (((ε ξ))^(|m|/2))(e^(- ε ξ/2))(Subscript[L, n]^(|m|))(ε ξ)

while constants (Subscript[Z, i])are related to quantum numbers by

((Subscript[Z, 1])/ϵ)- ((|m|+1)/2)=(Subscript[n, 1]);    0≤ (Subscript[n, 1])∈ ℕ.

((Subscript[Z, 2])/ϵ)- ((|m|+1)/2)=(Subscript[n, 2]);    0≤ (Subscript[n, 2])∈ ℕ.

(Z/ϵ)= |m|+1 +(Subscript[n, 1])+ (Subscript[n, 2]) ≡ n ≥ 1;     ⇒  - ((Z^2)/(2E))= (n^2).

Numerical checks

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

Wave equation

Let us first write one eigenfunction

fb[ξ_,n_,m_,ε_]:=(2^(1/4)) ε(Sqrt[(( n!)/((n+m)!))]) (((ε ξ))^(m/2))Exp[- ε  ξ/2]LaguerreL[n,m,ε  ξ]

The energy eignevalues are E = -1/2(n^2) and this gives for ϵ the value ϵ = 1/n. The separation constants are then given by:

Z1[n1_,m_,n_]:= (n1 + ((m+1)/2))/n

Now we can verify that functions satisfy differential equation ()

FullSimplify[(∂_ξ)(ξ (∂_ξ)fb[ξ,n1,m,1/n]) + (- (1/(4 (n^2)))ξ - ((m^2)/(4 ξ))+ Z1[n1,m,n])fb[ξ,n1,m,1/n],
Assumptions-> {n1∈ Integers,n1>0,m∈Integers, m> 0,n∈Integers, n> 0}]

Orthogonality

We can check directly orthogonality computing integrals. The following lines check orthogonality for a given value of m and for a given maximum value of n, nMax

nMax=3; m=1;

scalar[i1_,j1_,i2_,j2_,m_]:= Module[{ε1,ε2,ρ},
ε1 = 1/(i1+j1+m+1);
ε2=1/(i2+j2+m+1);
ρ = 2/(ε1+ε2);
a1 = Integrate[fb[ x,i1,m,ε1]*x*fb[ x,i2,m,ε2],{x,0,∞},GenerateConditions-> False];
a2 = Integrate[fb[ x,j1,m,ε1]*fb[ x,j2,m,ε2],{x,0,∞},GenerateConditions-> False];
b1 =Integrate[fb[ x,i1,m,ε1]*fb[ x,i2,m,ε2],{x,0,∞},GenerateConditions-> False];
b2 = Integrate[x*fb[ x,j1,m,ε1]*fb[ x,j2,m,ε2],{x,0,∞},GenerateConditions-> False];
(a1*a2+b1*b2)/4]

((tab=Table[scalar[n1,3-n1,n2,3-n2,m],{n1,0,nMax},{n2,0,nMax}])//TableForm)