Variational solution to ground state hydrogen wave function

Remove["Global`*"]

$Assumptions={a>0,e<0,ℏ>0,m>0};

u=Exp[-r/a];

HavgNum=Integrate[
r^2u ((-ℏ^2/(2m))(1/r^2)D[r^2D[u,r],r]-(e^2/r) u),{r,0,∞}];
HavgDen=Integrate[r^2 u^2,{r,0,∞}];Havg=HavgNum/HavgDen

aSolve=a/.Solve[D[Havg,a]==0,a][[1]]

Havg/.a->aSolve

Havg/.a->aSolve/.m->ℏ^2/(a0 e^2)

Havg/.a->aSolve/.ℏ->e^2/(α c)

Variational solution infinite potential well

Remove["Global`*"]

$Assumptions={a>0,m>0,ℏ>0,norm>0,β>1/2};

ψExact=(Sqrt[1 /a])Cos[Pi x/(2a)];
Integrate[ψExact^2,{x,-a,a}]==1
EExact=Integrate[-ℏ^2/(2m)ψExact D[ψExact,{x,2}],{x,-a,a}]

ψApprox1=norm(a^2-x^2);
HavgNum=Integrate[-ℏ^2/(2m)ψApprox1 D[ψApprox1,{x,2}],{x,-a,a}];
HavgDen=Integrate[ψApprox1^2,{x,-a,a}];Havg1=HavgNum/HavgDen

N[(Havg1-EExact)/EExact]

ψApprox2=norm(a^2-x^2)^β;
HavgNum=Integrate[-ℏ^2/(2m)ψApprox2 D[ψApprox2,{x,2}],{x,-a,a}];
HavgDen=Integrate[ψApprox2^2,{x,-a,a}];Havg2=HavgNum/HavgDen

HavgScale=Havg2/(ℏ^2/(a^2 m));
Plot[HavgScale,{β,1,1.5}]

βMin=β/.FindRoot[D[HavgScale,β],{β,1.0}]

Havg2β=Havg2/.β->βMin

N[(Havg2β-EExact)/EExact]

ψApprox1n=Simplify[ψApprox1/Sqrt[Integrate[ψApprox1^2,{x,-a,a}]]];
ψApprox2n=Simplify[ψApprox2/Sqrt[Integrate[ψApprox2^2,{x,-a,a}]]];

aVal=1;
Plot[
{ψExact/.a->aVal,
ψApprox1n/.a->aVal,
ψApprox2n/.a->aVal/.β->βMin},
{x,-aVal,aVal},PlotLegends->{Exact,Quadratic,Minimized}]

Plot[
{((ψApprox1n/.a->aVal)-(ψExact/.a->aVal))/ψExact/.a->aVal,
((ψApprox2n/.a->aVal/.β->βMin)-(ψExact/.a->aVal))/ψExact/.a->aVal},
{x,-aVal,aVal}]

N[Integrate[ψExact ψExact,{x,-a,a}]]
N[Integrate[ψExact ψApprox1n,{x,-a,a}]]
N[Integrate[ψExact ψApprox2n/.β->βMin,{x,-a,a}]]

Variational approach to hydrogen polarizability

Remove["Global`*"]

$Assumptions={a0>0,m>0,ℏ>0,e<0,Ef>0};

First do the regular hydrogen atom to make sure we know what we’re doing

ψnlm[r_,θ_,ϕ_]=Nnl Rnl[r]Ylm[θ,ϕ];

Ylm[θ_,ϕ_]=SphericalHarmonicY[l,m,θ,ϕ];
Rnl[r_]=(2 Z r/(n a0))^l Exp[-Z r/(n a0)]Hypergeometric1F1[-n+l+1,2l+2,2Z r/(n a0)];
Nnl=1/(2l+1)! Sqrt[(2Z/(n a0))^3(n+l)!/(2n(n-l-1)!)];

Z=1;
ψ=ψnlm[r,θ,ϕ]/.{n->1,l->0,m->0}

Integrate[ψ^2 r^2Sin[θ],{r,0,∞},{θ,0,Pi},{ϕ,0,2Pi}]==1

z=r Cos[θ];
E0=Integrate[
ψ(
(-ℏ^2/(2m))((1/r^2)D[r^2D[ψ,r],r]+1/(r^2Sin[θ])D[Sin[θ]D[ψ,θ],θ])
-(e^2/r)ψ-e Ef z ψ
)r^2Sin[θ],
{r,0,∞},{θ,0,Pi},{ϕ,0,2Pi}]/.m->ℏ^2/(a0 e^2)

Now use a trial wave function that “sloshes” linearly to one side. Note that we expect δ<0 because the electron will be pushed in the direction opposite to the electric field (which is in the +z direction).

ψTrial=ψ(1+δ z/a0);
HavgNum=Integrate[
ψTrial(
(-ℏ^2/(2m))((1/r^2)D[r^2D[ψTrial,r],r]+1/(r^2Sin[θ])D[Sin[θ]D[ψTrial,θ],θ])
-(e^2/r)ψTrial-e Ef z ψTrial
)r^2Sin[θ],
{r,0,∞},{θ,0,Pi},{ϕ,0,2Pi}];
HavgDen=Integrate[
ψTrial ψTrial r^2Sin[θ],
{r,0,∞},{θ,0,Pi},{ϕ,0,2Pi}];
Havg=HavgNum/HavgDen/.m->ℏ^2/(a0 e^2);

Make a simple check

(Havg/.δ->0)==E0

Find the minimum with respect to the free parameter δ (which must be negative) and substitute into Havg

δMin=δ/.Solve[D[Havg,δ]==0,δ][[1]]

HavgMin=Series[Havg/.δ->δMin,{Ef,0,3}]