One-Electron Atomic Wave Functions

Remove["Global`*"]

Formulas follow Sakurai & Napolitano, 3e, Appendix Section B.6

Set the assumptions and form the wave function with the normalization factor

$Assumptions={Z>0,a0>0};

ψ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)!)];

Build some different wave functions and check the normalization integrals

ψ1s=ψnlm[r,θ,ϕ]/.{n->1,l->0,m->0};
Integrate[r^2 Sin[θ]Abs[ψ1s]^2,{r,0,∞},{θ,0,Pi},{ϕ,0,2Pi}]==1

ψ2s=ψnlm[r,θ,ϕ]/.{n->2,l->0,m->0};
Integrate[r^2 Sin[θ]Abs[ψ2s]^2,{r,0,∞},{θ,0,Pi},{ϕ,0,2Pi}]==1

ψ2p0=ψnlm[r,θ,ϕ]/.{n->2,l->1,m->0};
Integrate[r^2 Sin[θ]Abs[ψ2p0]^2,{r,0,∞},{θ,0,Pi},{ϕ,0,2Pi}]==1

ψ2p1p=ψnlm[r,θ,ϕ]/.{n->2,l->1,m->1};
Integrate[r^2 Sin[θ]Abs[ψ2p1p]^2,{r,0,∞},{θ,0,Pi},{ϕ,0,2Pi}]==1

ψ2p1m=ψnlm[r,θ,ϕ]/.{n->2,l->1,m->-1};
Integrate[r^2 Sin[θ]Abs[ψ2p1m]^2,{r,0,∞},{θ,0,Pi},{ϕ,0,2Pi}]==1

ψ3s=ψnlm[r,θ,ϕ]/.{n->3,l->0,m->0};
Integrate[r^2 Sin[θ]Abs[ψ3s]^2,{r,0,∞},{θ,0,Pi},{ϕ,0,2Pi}]==1

Stick with hydrogen and scale distance by the Bohr radius

Z=1;
a0=1;

Plot the radial wave functions

Plot[Nnl Rnl[r]/.{n->1,l->0},{r,0,6},
LabelStyle->15,PlotStyle->{Thick,Black},PlotRange->{0,2},PlotLegends->"1s"]

Plot[
{Nnl Rnl[r]/.{n->2,l->0},
Nnl Rnl[r]/.{n->2,l->1}},{r,0,10},
PlotRange->{-0.2,0.8},
LabelStyle->15,PlotStyle->Thick,
PlotLegends->{"2s","2p"}]

Plot[
{Nnl Rnl[r]/.{n->3,l->0},
Nnl Rnl[r]/.{n->3,l->1},
Nnl Rnl[r]/.{n->3,l->2}},{r,0,20},
PlotRange->{-0.1,0.4},
LabelStyle->15,PlotStyle->Thick,PlotLegends->{"3s","3p","3d"}]

These plots agree with Figure 7.7 in Krane, Modern Physics, Third Edition

Now make the 3D plots.

coordrep={
r->ToSphericalCoordinates[{x,y,z}][[1]],
θ->ToSphericalCoordinates[{x,y,z}][[2]],
ϕ->ToSphericalCoordinates[{x,y,z}][[3]]};

I use PerformanceGoal -> Speed just to get it to run faster

lim=5;
ContourPlot3D[
r^2Abs[ψ1s]^2/.coordrep,
{x,-lim,lim},{y,-lim,lim},{z,-lim,lim},
PerformanceGoal->"Speed"]

Using PerformanceGoal -> Quality gives a nicer plot, but takes MUCH longer for the last two (below)

lim=5;
ContourPlot3D[
r^2Abs[ψ1s]^2/.coordrep,
{x,-lim,lim},{y,-lim,lim},{z,-lim,lim},
PerformanceGoal->"Quality"]

Finish off the last two with PerformanceGoal -> Speed

lim=10;
ContourPlot3D[
r^2 Sin[θ]Abs[ψ2p1p]^2/.coordrep,
{x,-lim,lim},{y,-lim,lim},{z,-lim,lim},
PerformanceGoal->"Speed"]

lim=10;
ContourPlot3D[
r^2 Sin[θ]Abs[ψ2p0]^2/.coordrep,
{x,-lim,lim},{y,-lim,lim},{z,-lim,lim},
PerformanceGoal->"Speed"]