Polarization of hydrogen atom by Dalgarno-Lewis method

Quantum Mechanics
                       K. Konishi, G. Paffuti

Use the Dalgarno-Lewis method to compute up to order (ℰ^4) the polarization of the ground state Hydrogen atom.

Utilities

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

Introduction

In problems of this chapter it has been shown that the computation of higher order corrections due to an electric field to ground state wave function in a central potential can be computed using the expansion

ψ = (Subscript[ψ, 0])+ Cos[θ] (Subscript[F, 1])[r](Subscript[ψ, 0])[r] + ((Subscript[P, 2])[Cos[θ]](Subscript[F, 22])[r]+ (Subscript[F, 20])[r])(Subscript[ψ, 0])+ …

Using the fact that first order correction is zero for parity reasons Schrödinger equation for ψ can be transformed in a series of inhomogeneous differential equations

( (1/2)(Subscript[F, 1]^(''))+(1/r)(Subscript[F, 1]^('))- (2/(r^2))(Subscript[F, 1])+((Subscript[ψ, 0]^('))/(Subscript[ψ, 0]))(Subscript[F, 1]^(')))+ ((e ℰ m)/(ℏ^2))r = 0;

( (1/2)(Subscript[F, 22]^(''))+(1/r)(Subscript[F, 22]^('))- (6/(r^2))(Subscript[F, 22])+((Subscript[ψ, 0]^('))/(Subscript[ψ, 0]))(Subscript[F, 22]^(')))+ (2/3)((e ℰ m)/(ℏ^2))r (Subscript[F, 1])= 0;

( (1/2)(Subscript[F, 20]^(''))+(1/r)(Subscript[F, 20]^('))+((Subscript[ψ, 0]^('))/(Subscript[ψ, 0]))(Subscript[F, 20]^(')))+ (1/3)((e ℰ m)/(ℏ^2))r (Subscript[F, 1]) + (Subscript[ϵ, 2])(m/(ℏ^2))=0

In the following we use atomic units, with e = ℏ=m=1. To come back to usual unities it is sufficient to note that the dimensionless expansion parameter is (a is the Bohr radius).

ℱ = e ℰ a/((e^2)/a).

We want to compute corrections to the ground state of hydrogen atom, for which  (Subscript[ψ, 0])∝ Exp[-r/a]. In this case previous equations for hydrogen atom take the simpler form

( (1/2)(Subscript[F, 1]^(''))+(1/r)(Subscript[F, 1]^('))- (1/(r^2))(Subscript[F, 1])-(Subscript[F, 1]^(')))+ ℰ r = 0;

( (1/2)(Subscript[F, 22]^(''))+(1/r)(Subscript[F, 22]^('))- (3/(r^2))(Subscript[F, 22])-(Subscript[F, 22]^(')))+ (2/3)ℰ r (Subscript[F, 1])= 0;

( (1/2)(Subscript[F, 20]^(''))+(1/r)(Subscript[F, 20]^('))-(Subscript[F, 20]^(')))+ (1/3)ℰ r (Subscript[F, 1]) + (Subscript[ϵ, 2])=0

Perturbation computation

First order in ψ

The equation to be solved is

( (1/2)(Subscript[F, 1]^(''))+(1/r)(Subscript[F, 1]^('))- (1/(r^2))(Subscript[F, 1])-(Subscript[F, 1]^(')))+ ℰ r = 0;

We write the wave function without the inessential factor 1/(Sqrt[2π]) which just cancel the trivial azimuthal angle integral. The overall 1/(Sqrt[4π])normalization factor do to angular variables reduces to 1/(Sqrt[2]).

ψ0[r_]:= 2 Exp[-r]/(Sqrt[2]);

We look for a simple polynomial solution

maxedg=2;

fp[r_]:= Sum[a[k](r^k),{k,0,maxedg}];
varA = Table[a[i],{i,0,maxedg}];

The equation is

eq1=CoefficientList[Expand[(r^2)((1/2)fp''[r]+(1/r)fp'[r]-(1/(r^2))fp[r]- fp'[r]+ℰ r)],r];

We use Solve to determine the coefficients :

a1=varA/.First[Solve[eq1== 0,varA]]

The radial function F1[r] is then

F1[r_]:=  Sum[a1[[k+1]](r^k),{k,0,2}]

For simplicity we put Cos[θ] → t, then the first order correction to wave function is

ψ1[r_,t_]:= t*F1[r]*ψ0[r]

Second order correction to energy level is given by

(Subscript[ϵ, 2]) = ⟨(Subscript[ψ, 0])|- ℰ r t |(Subscript[ψ, 1])⟩

In performing angular integration we have -1 ≤ t = Cos[θ] ≤ 1.

ϵ2=-ℰ Integrate[(r^2)*ψ1[r,t]*ψ0[r] r t,{t,-1,1},{r,0,∞}]

Second order in ψ

The equations to be solved are :

( (1/2)(Subscript[F, 22]^(''))+(1/r)(Subscript[F, 22]^('))- (3/(r^2))(Subscript[F, 22])-(Subscript[F, 22]^(')))+ (2/3)ℰ r (Subscript[F, 1])= 0;

( (1/2)(Subscript[F, 20]^(''))+(1/r)(Subscript[F, 20]^('))-(Subscript[F, 20]^(')))+ (1/3)ℰ r (Subscript[F, 1]) + (Subscript[ϵ, 2])=0

We proceed as before

maxedg=4;

fp[r_]:= Sum[a[k](r^k),{k,0,maxedg}];
varA = Table[a[i],{i,0,maxedg}];
eq1=CoefficientList[Expand[(r^2)((1/2)fp''[r]+(1/r)fp'[r]-(3/(r^2))fp[r]- fp'[r]+(2/3)ℰ r F1[r])],r]

a22=varA/.First[Solve[eq1== 0,varA]]

This gives the function (Subscript[F, 22])[r]

F22[r_]:=  Sum[a22[[k+1]](r^k),{k,0,maxedg}]

Let us now consider (Subscript[F, 20])[r]

fp[r_]:= Sum[a[k](r^k),{k,0,maxedg}];
varA = Table[a[i],{i,0,maxedg}];
eq1=CoefficientList[Expand[(r^2)((1/2)fp''[r]+(1/r)fp'[r]- fp'[r]+(1/3)ℰ r F1[r]+ϵ2)],r];

a20=varA/.First[Solve[eq1== 0,varA]]

The coefficient a[0] has not been fixed, because this part of the solution is not automatically orthogonal to (Subscript[ψ, 0]). We have to fix a[0] by imposing the orthogonality:

F20[r_]:=  Sum[a20[[k+1]](r^k),{k,0,maxedg}]

eqfora0=First[Solve[Integrate[F20[r]*((ψ0[r])^2)r^2,{r,0,∞}]== 0,a[0]]]

a20[[1]]= a[0]/.eqfora0;

The function (Subscript[F, 20])[r]in now completely fixed:

F20[r_]:=  Sum[a20[[k+1]](r^k),{k,0,maxedg}]

The full second order correction is :

ψ2[r_,t_]:=( F20[r] + LegendreP[2,t]F22[r])ψ0[r]

From (Subscript[ψ, 2])we can compute fourth order Stark effect using perturbation theory recurrence relation

(Subscript[ϵ, 4])= ⟨(Subscript[ψ, 1])|V|(Subscript[ψ, 2])⟩ - (Subscript[ϵ, 2])⟨(Subscript[ψ, 1])|(Subscript[ψ, 1])⟩

ϵ4=-ℰ*Integrate[ψ2[r,t] r t ψ1[r,t] r^2,{t,-1,1},{r,0,∞}] - ϵ2*Integrate[((ψ1[r,t])^2)  r^2,{t,-1,1},{r,0,∞}]