Quantum Mechanics of the hydrogen atom

Load Packages :

<<PhysicalConstants`
<<VectorAnalysis` (*for solving the differential equations in spherical coordinates*)

Set default coordinate system to spherical coordinates adapted to the problem

Laplacian[f[Xx,Yy,Zz]]

SetCoordinates[Spherical[r,θ,ϕ]];

Laplacian[f[r,θ,ϕ]]

Hamilton Operator:

H = ((p^2)/(2μ ))+V(r)=- ((ℏ^2)/(2μ))Δ+V(r)=-((ℏ^2)/(2 μ))(((∂^2)/(∂(r^2)))+(2/r)(∂/(∂r))+(1/(r^2))(Λ^2))+V(r)

with

reduced mass μ= ((mp me)/(me + mp))
(Λ^2)=(1/((sin[θ])^2))((∂^2)/(∂(ϕ^2)))+(1/(sin[θ]))(∂/(∂θ))sin[θ](∂/(∂θ))

(Λ^2)... Legendre Operator corresponds to the orbital angular momentum operator

(l^2)=-(ℏ^2)(Λ^2)=-(ℏ^2)((1/((sin[θ])^2))((∂^2)/(∂(ϕ^2)))+(1/(sin[θ]))(∂/(∂θ))sin[θ](∂/(∂θ)))

We have to find solutions to the (stationary) Schrödinger equation

H Ψ(r,θ,ϕ) = E Ψ(r,θ,ϕ)

The eigenvalues E correspond to the energy of the corresponding eigenstates, the wavefunctions of the system.

HHamiltonian[ψ_]:=-((ℏ^2)/(2 μ)) Laplacian[ψ]+V[r]ψ

Coulomb Potential V(r) (spherical symmetric 1/r potential)

V[r_,Z_:1]:=-((Z (Q^2))/(4π (Subscript[ϵ, 0])r))

Z...charge of the nucleus (here: Z=1, but can be adapted to include ions)
Q... electron charge

HHamiltonian[ψ[r,θ,ϕ]]

Separationsansatz (stationary problem) :

Since the Hamilton operator can be split up into a part depending only on r and a second part only on θ and ϕ we can try the separationsansatz

Ψ[r, θ, ϕ] -> R[r] Y[θ, ϕ]

The Hamiltonian is then

H R[r] Y[θ,ϕ] = E R[r] Y[θ,ϕ]

-((ℏ^2)/(2 μ))(Y(((∂^2)R)/(∂(r^2)))+((2Y)/r)((∂R)/(∂r))+(R/(r^2))(Λ^2)Y)+V(r) R Y = E R Y

Multiplication with ((r^2)/(R Y)) leads to

-((ℏ^2)/(2 μ))(((r^2)/R)(((∂^2)R)/(∂(r^2)))+((2r)/R)((∂R)/(∂r))+(1/Y)(Λ^2)Y)+V(r) (r^2)= E (r^2)

or

-((ℏ^2)/(2 μ))(((r^2)/R)(((∂^2)R)/(∂(r^2)))+((2r)/R)((∂R)/(∂r)))+(V(r)-E) (r^2)= λ=((ℏ^2)/(2 μY))(Λ^2)Y

Ψ[r,θ,ϕ]=R[r]Y[θ,ϕ];

HHamiltonian[Ψ[r,θ,ϕ]]-En Ψ[r,θ,ϕ]

Angular Part :

We have to solve the differential equation

λ=((ℏ^2)/(2 μY))(Λ^2)Y

AngularPart=(HHamiltonian[Y[θ,ϕ]]-V[r]Y[θ,ϕ])((2(r^2)μ)/(ℏ^2))

DSolve[AngularPart==λ Y[θ,ϕ],Y,{θ,ϕ}]

Unfortunately, DSolve cannot solve this type of partial differential equation. But we can split the equation once more into two parts (or look up the solution in any textbook on quantum mechanics).

Y[θ_,ϕ_]=T[θ]Φ[ϕ];

(AngularPart-λ Y[θ,ϕ])==0//Simplify

((Φ[ϕ] (Cot[θ] T'[θ]+T''[θ])+T[θ] (λ Φ[ϕ]+((Csc[θ])^2) Φ''[ϕ]))/(Φ[ϕ]T[θ]((Csc[θ])^2)))==0//Simplify

or

λ sin(((θ))^2)+ ((sin(θ))/(T(θ)))(cos(θ)((∂T(θ))/(∂θ))+sin(θ)(((∂^2)T(θ))/(∂(θ^2))))=-(((∂^2)Φ(ϕ))/(∂(ϕ^2)))=α

ϕ - part

DSolve[-((Φ''[ϕ])/(Φ[ϕ]))==α,Φ[ϕ],ϕ]

Uniqueness of the solutions on the sphere require that the solution is periodic in 2 π -> α =(m^2) with m =0,±1,±2,...
and by proper choice of normalization constants we find

Φ[ϕ_,m_]=(1/(Sqrt[2π]))Exp[I m ϕ]

Plugging α=(m^2)into the equation above yields

(λ ((Sin[θ])^2)+((Sin[θ] (Cos[θ] T'[θ]+Sin[θ] T''[θ]))/(T[θ]))-(m^2))T[θ]==0//Simplify

DSolve[((m^2)-λ ((Sin[θ])^2)) T[θ]==Sin[θ] (Cos[θ] T'[θ]+Sin[θ] T''[θ]),T[θ],θ]

Legendre Polynomials LegendreP[l,m,Cos[θ]]

Plot[Evaluate[Flatten[Table[Table[Tooltip[LegendreP[l,m,Cos[θ]],"l="<>ToString[l]<>", m="<>ToString[m]],{m,-l,l}],{l,0,1.5,.5}]]],{θ,-π,π},PlotRange->{-4,4},Frame->True]

-> integer l gives condition λ = l (l + 1) to obtain polynomial

(1/2) (-1+(Sqrt[1+4 λ]))/.λ->l(l+1)//Simplify[#,l>0]&

l ... angular momentum quantum number

Flatten[Table[Table[LegendreP[l,m,x],{m,-l,l}],{l,0,3}]]
Flatten[Table[Table[LegendreQ[l,m,x],{m,-l,l}],{l,0,2}]]

Plot[Evaluate[Flatten[Table[Table[LegendreP[l,m,x],{m,-l,l}],{l,0,3}]]],{x,-1,1},PlotRange->All]
Plot[Evaluate[Flatten[Table[Table[LegendreQ[l,m,x],{m,-l,l}],{l,0,3}]]],{x,-1,1},PlotRange->All]

Normalizable solutions: C . LegendreP

T[θ_,l_,m_]=C LegendreP[l,m,Cos[θ]]

Spherical Harmonics :

The solution of the angular part of the Schrödinger equation is given in terms of spherical harmonics:

Y[θ, ϕ] = T[θ, l, m] Φ[ϕ, m]∝ (e^(-i m ϕ))((Subscript[P, l])^m)(cos θ)

(Subscript[Y, lm]) are spherical harmonic function (of degree l and order m) with l=0,1,2,... and m=0,+±1,±2,...
l... angular momentum quantum number
m...magnetic quantum number

Calculate the spherical harmonics :

Y[l_Integer,m_Integer,θ_,ϕ_]:=SphericalHarmonicY[l,m,θ,ϕ]

Table[Table[Y[l,m,θ,ϕ],{m,-l,l}],{l,0,3}]

2D Plots in xz - plane

Modulus of the functions :

Table[Table[PolarPlot[Evaluate[Abs[Y[l,m,θ,0]]],{θ,0,2π},PlotRange->{{-0.65,0.65},{-0.65,0.65}},ImageSize->2*72,FrameLabel->{"x","z"},Frame->True,FrameTicks->None],{m,-l,l}],{l,0,2}]

Orthogonality and Normalization :

∫(Subscript[Y, lm]^( * ))(Subscript[Y, l'm'])\[DifferentialD]Ω=(Subscript[δ, ll'])(Subscript[δ, m m'])

Integrate[Y[1,1,t,phi]Conjugate[Y[1,1,t,phi]]Sin[t],{t,0,π},{phi,0,2π}]
Integrate[Y[1,1,t,phi]Conjugate[Y[1,0,t,phi]]Sin[t],{t,0,π},{phi,0,2π}]
Integrate[Y[2,1,t,phi]Conjugate[Y[2,1,t,phi]]Sin[t],{t,0,π},{phi,0,2π}]

Radial Part (bound solutions E<0):

((ℏ^2)/(2 μ))l(l+1) from the solution of the angular part:

-((ℏ^2)/(2 μ))(((r^2)/R)(((∂^2)R)/(∂(r^2)))+((2r)/R)((∂R)/(∂r)))+(V(r)-E) (r^2)+((ℏ^2)/(2 μ))l(l+1)==0

Multiplication with (R/(r^2))

-((ℏ^2)/(2 μ))((((∂^2)R)/(∂(r^2)))+(2/r)((∂R)/(∂r)))+(V(r)-E) R-(((ℏ^2) R)/(2 μ (r^2)))l(l+1)==0

and inserting V(r) = -((Z (Q^2))/(4 π (Subscript[ϵ, 0])r))

-((ℏ^2)/(2 μ))((((∂^2)R)/(∂(r^2)))+(2/r)((∂R)/(∂r)))+((((ℏ^2) )/(2 μ (r^2)))l(l+1)-((Z (Q^2))/(4 π (Subscript[ϵ, 0])r))-E) R==0

effective Potential: V(r) = (((ℏ^2) )/(2 μ (r^2)))l(l+1)-((Z (Q^2))/(4 π (Subscript[ϵ, 0])r))

R[r,0,0]=Rr[r];
Ψ[r,0,0]=Rr[r];

HHamiltonian[Rr[r,0,0]]

HHamiltonian[R[r,0,0]]-En Ψ[r,0,0]+((l(l+1) Rr[r])/(r^2))

The differential equation for the radial part is finally

(k^2) R(r)+((R(r) (Subscript[l, orb]))/(r^2))+((R(r) (Subscript[l, orb]^2))/(r^2))-((2 R(r))/(a r))-(2/r)((∂R(r))/(∂r))-(((∂^2)R(r))/((∂^2)r))]

with 
En =- (((ℏ^2)(k^2))/(2 μ))(relating a wavevector to the energy),
a= (((ℏ^2)4π (Subscript[ϵ, 0]))/(μ (Q^2))) (Bohr radius)and 
Z=1 (for the hydrogen)

bohrradius=(((ℏ^2)4π (Subscript[ϵ, 0]))/(μ (Q^2)));

RadialDiffEquation=((HHamiltonian[Rr[r]]-En Rr[r]))2 μ/(ℏ^2)+((l(l+1) Rr[r])/(r^2))/.{En->-(ℏ^2)(k^2)/(2 μ), (Q^2)-> (ℏ^2)4π (Subscript[ϵ, 0])/(a  μ)}//Simplify//Expand

radialsolution=Simplify[DSolve[RadialDiffEquation==0,Rr[r],r]]

HypergeometricU... confluent hypergeometric function U(a,b,z)
LaguerreL...generalized Laguerre polynomial (Subscript[L, n]^a)(x)

Solution at r = 0 :

We obtained two linearly independent solutions of the differential equations. Now we have to investigate which one is physical.

sol1=(E^(-k r)) (r^l)HypergeometricU[1-(( 1)/( a k))+l,2 (1+l),2 k r];
sol2=(E^(-k r)) (r^l)LaguerreL[-1+(( 1)/(a k))-l,1+2 l,2 k r];

TaylorRadial1=Series[sol1,{r,0,1}]//FullSimplify[#,Assumptions->{k>0,r>0,l==1}]&
TaylorRadial2=Series[sol2,{r,0,1}]//FullSimplify

Solution  Part 1 :

The first solution containing the confluent hypergeometric function goes to infinity at the origin, thus can be excluded.

Plot[Evaluate[Table[sol1//.{a->1,k->2},{l,0,3}]],{r,0,1},PlotRange->{0,10},FrameLabel->{"r [(Subscript[a, 0])]","R[r]"},Frame->True]

Solution  Part 2 :

The second solution is finite at the origin.

Plot[Evaluate[Table[sol2//.{a->1,k->2},{l,0,3}]],{r,0,2},PlotRange->{-0.5,0.5},FrameLabel->{"r [(Subscript[a, 0])]","R[r]"},Frame->True]

However, the wave functions are apparently not normalizable since they are not falling off exponentially as r -> ∞.

Plot[Evaluate[Table[sol2//.{a->1,k->2},{l,0,3}]],{r,0,10},PlotRange->All,FrameLabel->{"r [(Subscript[a, 0])]","R[r]"},Frame->True]

The trick is to know that the Laguerre function LaguerreL[a, b, c] found is polynomial whenever a is a positive integer.

This leads to the condition

wavenumber=Solve[-1+(( 1)/(a k))-l==nr,k][[1]]

sol2/.wavenumber/.a->1

Plot[Evaluate[Table[sol2/.wavenumber/.{a->1,nr->4},{l,0,3}]],{r,0,50},PlotRange->All,FrameLabel->{"r [(Subscript[a, 0])]","R[r]"},Frame->True]

Energy Quantum number: n = nr+l+1

The solution  which is finite at the origin and normalizable is finally given by :

R[r_,n_,l_]=Cr sol2/.wavenumber//.{nr->n-l-1}

where C is the normalization constant which has still to be calculated.

Energy eigenvalues

En[n_,l_,a_]=-(((ℏ^2)(k^2))/(2 μ))/.wavenumber/.{nr->n-l-1}

En[n,l,a]/.{ a->bohrradius}

Plug in numbers for n=1, l=0 yields the energy of the 1s bound state. This is the energy one needs to kick out the electron (ionization energy).

En[1,0,bohrradius]/ElectronCharge  Coulomb/.{ℏ->PlanckConstant/(2π)/Joule/Second,(Subscript[ϵ, 0])->VacuumPermittivity Meter Volt/Ampere/Second,Q->ElectronCharge/Coulomb,μ->(ElectronMass*ProtonMass)/(ElectronMass+ProtonMass)/Kilogram}

Energy of the ground state of the hydrogen atom: (Subscript[E, 1])=-13.6 eV

Normalization of the wavefunction - determine C

∫(Subscript[((SuperStar[Ψ])), nlm ])(Subscript[((SuperStar[Ψ])), nlm ])(d^3)r = 1

∫(Subscript[R, nl])(r)(Subscript[R, nl])(r)(r^2)dr ∫(SuperStar[(Subscript[Y, lm])])(θ,ϕ)(Subscript[Y, lm])(θ,ϕ)sin[θ]dθ dϕ= 1

∫(SuperStar[(Subscript[Y, lm])])(θ,ϕ)(Subscript[Y, lm])(θ,ϕ)sin[θ]dθ dϕ =1

(∫_0^∞)((R[r,1,0])^2)(r^2)\[DifferentialD]r

Solve[Integrate[Evaluate[((R[r,1,0])^2)](r^2),{r,0,∞},Assumptions->{a>0}]==1,Cr]

Normalize for different values of l and n

Assign function CNorm[n,l] for each value of n and l as far as we need it.

CNorm[1,0,a_]=Cr/.Solve[Integrate[Evaluate[((R[r,1,0])^2)](r^2),{r,0,∞},Assumptions->{a>0}]==1,Cr][[2]]

Calculate normalization constants for n=1,2,3,4,5:

Table[Table[
CNorm[i,j,a_]=Cr/.Solve[Integrate[Evaluate[((R[r,i,j])^2)](r^2),{r,0,∞},Assumptions->{a>0}]==1,Cr][[2]],
{j,0,i-1}],
{i,1,6}]

Rnl[r_,n_,l_,a_:1]:=CNorm[n,l,a] (E^(-(r/(a n)))) (r^l) LaguerreL[-1+n-l,1+2 l,((2 r )/(a n))]
Rnlun[r_,n_,l_,a_:1]:= (E^(-(r/(a n)))) (r^l) LaguerreL[-1+n-l,1+2 l,((2 r )/(a n))]

Ψ[r_,θ_,ϕ_,n_,l_,m_,a_:1]:=Rnl[r,n,l,a]Y[l,m,θ,ϕ]
Ψun[r_,θ_,ϕ_,n_,l_,m_,a_:1]:=Rnlun[r,n,l,a]Y[l,m,θ,ϕ]

Ψ[r,θ,ϕ,n,l,m,a]
Ψ[r,θ,ϕ,2,0,0]

Energies of the eigenstates

Solve the equation H Ψ = E Ψ with the eigenfunctions found :

Solve[HHamiltonian[Evaluate[Ψ[r,θ,ϕ,1,0,0,bohrradius]]]==E1 Ψ[r,θ,ϕ,1,0,0,bohrradius],E1]//FullSimplify
Solve[HHamiltonian[Evaluate[Ψ[r,θ,ϕ,2,0,0,bohrradius]]]==E2 Ψ[r,θ,ϕ,2,0,0,bohrradius],E2]//FullSimplify
Solve[HHamiltonian[Evaluate[Ψ[r,θ,ϕ,2,1,0,bohrradius]]]==E21 Ψ[r,θ,ϕ,2,1,0,bohrradius],E21]//FullSimplify

Radial part of the wave function R

l = 0 :

Plot[Evaluate[Table[Rnl[r,n,0],{n,1,4}]],{r,0,20},PlotRange->All,FrameLabel->{"r[(Subscript[a, 0])]","(Subscript[R, nl])(r)"},Frame->True]

l = 1 :

Plot[Evaluate[Table[Rnl[r,n,1],{n,2,4}]],{r,0,30},PlotRange->All,FrameLabel->{"r[(Subscript[a, 0])]","(Subscript[R, nl])(r)"},Frame->True]

l = 2 :

Plot[Evaluate[Table[Rnl[r,n,2],{n,3,4}]],{r,0,40},PlotRange->All,FrameLabel->{"r[(Subscript[a, 0])]","(Subscript[R, nl])(r)"},Frame->True]

Probability of finding the electron at a distance between r and r+dr: 4π (r^2)|R((|)^2)

l = 0 :

Plot[Evaluate[Table[4π (r^2)((Rnl[r,n,0])^2),{n,1,4}]],{r,0,20},PlotRange->All,FrameLabel->{"r[(Subscript[a, 0])]","(Subscript[R, nl])(r)"},Frame->True]

Maximum at r=(Subscript[a, 0])

Plot[Evaluate[Table[4π (r^2)((Rnl[r,n,0])^2),{n,1,4}]],{r,0,2},PlotRange->All,FrameLabel->{"r[(Subscript[a, 0])]","(Subscript[R, nl])(r)"},Frame->True,GridLines->{{1},None}]

l = 1 :

Plot[Evaluate[Table[4π (r^2)((Rnl[r,n,0])^2),{n,2,4}]],{r,0,30},PlotRange->All,FrameLabel->{"r[(Subscript[a, 0])]","(Subscript[R, nl])(r)"},Frame->True]

l = 2 :

Plot[Evaluate[Table[4π (r^2)((Rnl[r,n,0])^2),{n,3,4}]],{r,0,40},PlotRange->All,FrameLabel->{"r[(Subscript[a, 0])]","(Subscript[R, nl])(r)"},Frame->True]

Electron Density at fixed distance r

ElectronDensity[r_,θ_,ϕ_,n_,l_,m_,a_:1]:=Conjugate[Ψ[r,θ,ϕ,n,l,m,a]]Ψ[r,θ,ϕ,n,l,m,a]
ElectronDensity[Ψ_]:=Conjugate[Ψ] Ψ

coordinateline=Graphics3D[{Thick,Line[{{-300,0,0},{300,0,0}}],Line[{{0,-300,0},{0,300,0}}],Line[{{0,0,-300},{0,0,300}}]}];
myplotrange={-0.06,0.06};

n = 1 (r=1(Subscript[a, 0]))

SphericalPlot3D[Evaluate[ElectronDensity[1,θ,ϕ,1,0,0]],{θ,0,π},{ϕ,0,2π},PlotRange->myplotrange,PlotStyle->Directive[Orange,Specularity[White,10]],Boxed->False,Axes->False];
Show[%,coordinateline]

n=2, l = 0 (r=1(Subscript[a, 0]))

SphericalPlot3D[Evaluate[
ElectronDensity[1,θ,ϕ,2,0,0]],{θ,0,π},{ϕ,0,2π},PlotRange->myplotrange,PlotStyle->Directive[Orange,Specularity[White,10]],Boxed->False,Axes->False]

Higher Density at the origin -> take r=0.05 (Subscript[a, 0]):

SphericalPlot3D[Evaluate[
ElectronDensity[0.05,θ,ϕ,2,0,0]],{θ,0,π},{ϕ,0,2π},PlotRange->myplotrange,PlotStyle->Directive[Orange,Specularity[White,10]],Boxed->False,Axes->False]

n=2, l = 1

Table[SphericalPlot3D[Evaluate[
ElectronDensity[1,θ,ϕ,2,1,m]],{θ,0,π},{ϕ,0,2π},PlotRange->All,PlotStyle->Directive[Orange,Specularity[White,10]],Boxed->False,Axes->False],{m,-1,1}]

n = 3, l = 0

SphericalPlot3D[Evaluate[
ElectronDensity[0.01,θ,ϕ,3,0,0]],{θ,0,π},{ϕ,0,2π},PlotRange->All,PlotStyle->Directive[Orange,Specularity[White,10]],Boxed->False,Axes->False]

n = 3, l = 1

Table[SphericalPlot3D[Evaluate[
ElectronDensity[1,θ,ϕ,3,1,m]],{θ,0,π},{ϕ,0,2π},PlotRange->All,PlotStyle->Directive[Orange,Specularity[White,10]],Boxed->False,Axes->False],{m,-1,1}]

n = 3,l = 2

Table[SphericalPlot3D[Evaluate[
ElectronDensity[5,θ,ϕ,3,2,m]],{θ,0,π},{ϕ,0,2π},PlotRange->All,PlotStyle->Directive[Orange,Specularity[White,10]],Boxed->False,Axes->False],{m,-2,2}]

n = 4, l = 0

SphericalPlot3D[Evaluate[
ElectronDensity[0.01,θ,ϕ,4,0,0]],{θ,0,π},{ϕ,0,2π},PlotRange->All,PlotStyle->Directive[Orange,Specularity[White,10]],Boxed->False,Axes->False]

n = 4, l = 1

Table[SphericalPlot3D[Evaluate[
ElectronDensity[2,θ,ϕ,4,1,m]],{θ,0,π},{ϕ,0,2π},PlotRange->All,PlotStyle->Directive[Orange,Specularity[White,10]],Boxed->False,Axes->False],{m,-1,1}]

n = 4, l = 3

Table[SphericalPlot3D[Evaluate[
100ElectronDensity[12,θ,ϕ,4,3,m]],{θ,0,π},{ϕ,0,2π},PlotRange->All,PlotStyle->Directive[Orange,Specularity[White,10]],Boxed->False,Axes->False],{m,-3,3}]

(Factor 100 for getting better plot, otherwise values too small)

n = 5, l = 4

Table[SphericalPlot3D[Evaluate[
1000ElectronDensity[12,θ,ϕ,5,4,m]],{θ,0,π},{ϕ,0,2π},PlotRange->All,PlotStyle->Directive[Orange,Specularity[White,10]],Boxed->False,Axes->False],{m,-3,3}]

Representation by a density plot of |Ψ((|)^2)

Cut in the X - Z Plane:

DensityPlot[Evaluate[ElectronDensity[Ψ[(Sqrt[(z^2)+(x^2)]),ArcCos[z/(Sqrt[(z^2)+(x^2)])],0,5,4,1]]],{x,-40,40},{z,-40,40},(*ColorFunction->Function[{x,y,z},Hue[z]],*)PlotPoints->90,PlotRange->All,AxesLabel->{"x [(Subscript[a, 0])]","z [(Subscript[a, 0])]"}]

Hybridization :

For given n  all states have same energy independent of their angular momentum  quantum number l and magnetic quantum number m. Superpositions can be formed, which are also valid solutions of the Schrödinger Equation. For example we can construct functions with symmetry axis in x or y direction.

s2[r_,θ_,ϕ_]=Ψ[r,θ,ϕ,2,0,0];
p2x[r_,θ_,ϕ_]=(Ψ[r,θ,ϕ,2,1,+1]+Ψ[r,θ,ϕ,2,1,-1])/((Sqrt[2])I);
p2y[r_,θ_,ϕ_]=(Ψ[r,θ,ϕ,2,1,+1]-Ψ[r,θ,ϕ,2,1,-1])/((Sqrt[2]));
p2z[r_,θ_,ϕ_]=Ψ[r,θ,ϕ,2,1,0];

2sp Orbital:

hybridcoordinates=Graphics3D[{Thick,Line[{{-0.005,0,0},{0.005,0,0}}],Line[{{0,-0.005,0},{0,0.005,0}}],Line[{{0,0,-0.005},{0,0,0.005}}]}];
SphericalPlot3D[Evaluate[ElectronDensity[s2[1,θ,ϕ]+p2x[1,θ,ϕ]]],{θ,0,π},{ϕ,0,2π},PlotPoints->30,Boxed->False,Axes->False,PlotStyle->Directive[Orange,Specularity[White,10]],PlotRange->All];
Show[%,%%]

In the XY - Plane :

2s, 2 px and 2 py orbitals :

Plot[ElectronDensity[s2[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]]/.{x->0},{y,-10,10},PlotRange->All,Axes->False,Frame->True,FrameLabel->{"y [(Subscript[a, 0])]","|(Subscript[Ψ, 20])((|)^2)"}];
Plot[Evaluate[ElectronDensity[p2x[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]/.{x->0}]],{y,-10,10},PlotRange->All,Axes->False,Frame->True,FrameLabel->{"y [(Subscript[a, 0])]","|(Subscript[Ψ, 21  x])((|)^2)"}];
Plot[Evaluate[ElectronDensity[p2y[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]/.{y->0}]],{x,-10,10},PlotRange->All,Axes->False,Frame->True,FrameLabel->{"x [(Subscript[a, 0])]","|(Subscript[Ψ, 21  y])((|)^2)"}];
GraphicsGrid[{{%%%,%%,%}},ImageSize->700]

DensityPlot[Evaluate[ElectronDensity[s2[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]]],{x,-10,10},{y,-10,10},FrameLabel->{"x [(Subscript[a, 0])]","y [(Subscript[a, 0])]"},ColorFunction->"CherryTones",PlotPoints->50,PlotRange->{0,0.005}];
DensityPlot[Evaluate[ElectronDensity[p2x[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]]],{x,-10,10},{y,-10,10},FrameLabel->{"x [(Subscript[a, 0])]","y [(Subscript[a, 0])]"},ColorFunction->"CherryTones",PlotPoints->30,PlotRange->All];
DensityPlot[Evaluate[ElectronDensity[p2y[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]]],{x,-10,10},{y,-10,10},ColorFunction->"CherryTones",PlotRange->All,PlotPoints->30,FrameLabel->{"x [(Subscript[a, 0])]","y [(Subscript[a, 0])]"}];
GraphicsGrid[{{%%%,%%,%}},ImageSize->700]

2s+2px

Plot[ElectronDensity[s2[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]+p2x[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]]/.{x->0},{y,-10,10},PlotRange->All,Axes->False,Frame->True,FrameLabel->{"y [(Subscript[a, 0])]","|(Subscript[Ψ, 20])+(Subscript[Ψ, 21  x])((|)^2)"}]

DensityPlot[Evaluate[ElectronDensity[s2[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]+ p2y[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]]],{x,-6,6},{y,-6,6},ColorFunction->"CherryTones",PlotPoints->50,PlotRange->All,FrameLabel->{"x [(Subscript[a, 0])]","y [(Subscript[a, 0])]","|(Subscript[Ψ, 20])+(Subscript[Ψ, 21  x])((|)^2)"}]

2(sp^3) Orbital:

SphericalPlot3D[Evaluate[ElectronDensity[s2[1,θ,ϕ]+p2x[1,θ,ϕ]+p2y[1,θ,ϕ]+p2z[1,θ,ϕ]]],{θ,0,π},{ϕ,0,2π},PlotPoints->30,PlotRange->All,Boxed->False,Axes->False,PlotStyle->Directive[Orange,Specularity[White,10]]];
Show[%,hybridcoordinates,ViewPoint->{-2,2,1}]

DensityPlot[Evaluate[ElectronDensity[s2[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]+ p2x[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]+ p2y[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]+ p2z[(Sqrt[(y^2)+(x^2)]),π/2,ArcTan[x,y]]]],{x,-7,7},{y,-7,7},ColorFunction->"CherryTones",PlotPoints->50,PlotRange->All,FrameLabel->{"x [(Subscript[a, 0])]","y [(Subscript[a, 0])]","|(Subscript[Ψ, 20])+(Subscript[Ψ, 21  x])+(Subscript[Ψ, 21  y])+(Subscript[Ψ, 21  z])((|)^2)"}]

Beyond the Hydrogen Atom: 1/(r^2)-Potential

V2[r_]:=-((  (Q^2))/(4π (Subscript[ϵ, 0])r))-((A (Q^2))/(4π (Subscript[ϵ, 0])(r^2)))

HHamiltonianMod[ψ_]:=-((ℏ^2)/(2 μ)) Laplacian[ψ]+V2[r]ψ

RadialDiffEquation2=((HHamiltonianMod[R2[r]]-En R2[r]))2 μ/(ℏ^2)+(((Subscript[l, orb])((Subscript[l, orb])+1) R2[r])/(r^2))/.{En->-(ℏ^2)(k^2)/(2 μ), (Q^2)-> (ℏ^2)4π (Subscript[ϵ, 0])/(a  μ),Z->1}//Simplify//Expand

radialsolution2=Simplify[DSolve[RadialDiffEquation2==0,R2[r],r]]

-(1/2)+(1/(a k))-(1/2) (Sqrt[1-((8 A)/a)+4 (Subscript[l, orb])+4 (Subscript[l, orb]^2)]) shall be integer:

wavenumber2=Solve[-(1/2)+(1/(a k))-(1/2) (Sqrt[1-((8 A)/a)+4 (Subscript[l, orb])+4 (Subscript[l, orb]^2)])==nr,k][[1]]//Simplify

R2[r_,n_,l_]=C LaguerreL[-(1/2)+(1/(a k))-(1/2) (Sqrt[1-((8 A)/a)+4 (Subscript[l, orb])+4 (Subscript[l, orb]^2)]),(Sqrt[1-((8 A)/a)+4 (Subscript[l, orb])+4 (Subscript[l, orb]^2)]),2 k r]/.wavenumber2//.{nr->n-(Subscript[l, orb])-1}/.{(Subscript[l, orb])->l}//Simplify

En2[n_,l_,a_,A_]=-(((ℏ^2)(k^2))/(2 μ))/.wavenumber2/.{nr->n-(Subscript[l, orb])-1}/.{(Subscript[l, orb])->l}//Simplify

Energy is now l - dependent!

Compare Energy Levels

rsqrenergy=Table[Table[En2[n,l,bohrradius,4 (10^(-12))],{l,0,n-1}],{n,1,5}]/ElectronCharge*Coulomb/.{ℏ->PlanckConstant/(2π)/Joule/Second,μ->ElectronMass/Kilogram,Q->ElectronCharge/Coulomb,(Subscript[ϵ, 0])->VacuumPermittivity Meter Volt/ Ampere/Second};
renergy=Table[Table[En[n,l,bohrradius],{l,0,n-1}],{n,1,5}]/ElectronCharge*Coulomb/.{ℏ->PlanckConstant/(2π)/Joule/Second,μ->ElectronMass/Kilogram,Q->ElectronCharge/Coulomb,(Subscript[ϵ, 0])->VacuumPermittivity Meter Volt/ Ampere/Second};

Graphics[{Text["1/r",{2,0.5}],Text["1/r+A/(r^2)",{9,0.5}],Text["E [eV]",{-4,0.5}],Text[Style[#,Larger],{-4,#-0.1}]&/@Flatten[renergy]}];
Graphics[{Red,Sequence@@(Line[{{0,#},{5,#}}]&/@Flatten[renergy])}];
Graphics[{Green,Sequence@@(Line[{{7,#},{12,#}}]&/@Flatten[rsqrenergy])}];
Show[%,%%,%%%]

-> Degeneracy of Energy Levels is lifted.